Hamlet Medina Ruiz

Michel Kieffer

Béatrice Pesquet-Popescu

TCP and Network Coding Equilibrium and Dynamic Properties

Keywords: congestion control, network coding, queue management, stability

This paper analyzes the impact on the TCP-Reno congestion control mechanism of a network coding (NC) layer inserted in the TCP/IP stack. A multi-source multi-link model is considered to study the equilibrium and dynamic properties of the TCP-NC protocol with RED as active queue management mechanism. The existence and uniqueness of some network equilibrium is demonstrated and characterized in terms of average throughput, loss rate, and queue length. Global stability is proved in absence of forward delay, and the effects of the NC redundancy factor and of the delay on the local stability of TCP-NC-RED is studied around the equilibrium point. Results provided by the proposed model are compared to those obtained by simulation for N sources sharing a single link. TCP-NC-RED becomes unstable when delay or capacity increases, as TCP-Reno does. Its stability region is characterized as a function of the redundancy factor.

I. INTRODUCTION

In [START_REF] Sundararajan | Network Coding Meets TCP[END_REF], a new TCP-friendly protocol is proposed that interfaces network coding (NC) with TCP by introducing a new coding layer in the TCP/IP stack between TCP and IP. The idea is to benefit from NC to improve TCP throughput. In this layer TCP segments are encoded at the sender and decoded at the receiver. In [START_REF] Sundararajan | Network Coding Meets TCP: Theory and Implementation[END_REF] the authors additionally describe several important practical aspects of incorporating NC below TCP-Reno [START_REF] Jacobson | Congestion avoidance and control[END_REF], such as the TCP-compatible sliding window code, and new rules for acknowledging bytes to the TCP layer. Redundant packets are transmitted by the sender to mask non-congestion losses, e.g., losses on the wireless part of the network, from TCP congestion control. Simulation results in [START_REF] Sundararajan | Network Coding Meets TCP: Theory and Implementation[END_REF], [START_REF] Sundararajan | Network Coding Meets TCP[END_REF] show an improvement in TCP throughput for a constant value of the redundancy factor ρ. The purpose of ρ is to match the rate at which data are obtained at the receiver to the sending rate of TCP. If ρ is too small, the losses are not masked from TCP and the NC decoding probability decreases. On the other hand if ρ is too large, the congestion may increase, reducing the throughput. Thus, an adaptation mechanism for ρ has been proposed in [START_REF] Medina | Redundancy adaptation scheme for network coding with TCP[END_REF].

This paper characterizes the effect of NC on the equilibrium and local stability of TCP-Reno. The equilibrium or average throughput has been characterized in [START_REF] Kim | Modeling network coded tcp throughput: a simple model and its validation[END_REF], where a framework 1 This work has been partly supported by DIM-LSC SWAN.

is proposed to predict TCP-NC throughput when Vegas is used as TCP congestion control. Here, we consider the tools from control and optimization theory used to characterize TCP [START_REF] Kelly | Rate Control for Communication Networks: Shadow Prices, Proportional Fairness and Stability[END_REF], [START_REF] Low | A Duality Model of TCP and Queue Management Algorithms[END_REF], [START_REF] Low | Understanding TCP Vegas: A Duality Model[END_REF], [START_REF] Tang | Equilibrium of Heterogeneous Congestion Control: Optimality and Stability[END_REF] to study TCP-NC as a distributed primal-dual algorithm that maximizes some aggregate utility function.

In what follows, TCP-NC is briefly described in Section II, details about the network model and properties of the TCP-NC protocol are given in Sections III and IV, respectively. Our model is evaluated in Section V and some conclusions and futures research directions are drawn in Section VI.

II. PRINCIPLES OF THE TCP-NC PROTOCOL

The TCP-NC protocol introduces a NC layer between the TCP and IP layers in the TCP/IP protocol stack [START_REF] Sundararajan | Network Coding Meets TCP[END_REF]. The sender stores packets generated by TCP in a coding buffer. For every arrival from TCP, it transmits ρ random linear combinations of the most w recently arrived packets in the coding buffer. To convey information about the linear combination, an NC header is added to the coded packets, which contains information about the coefficients used to mix the packets, the first and last byte SEQ of each packet involved in the mix, and the first byte that has not been acknowledged. The original packets remain in the coding buffer until an appropriate TCP ACK arrives from the receiver side.

On the receiver side, upon receiving a packet containing a new linear combination, the decoder places it in the decoding buffer, appends the corresponding coefficients to the decoding matrix (which contains the coefficients used to mix already received packets), and performs Gaussian elimination. This process helps determining whether the received packet contains new information. Then the receiver sends a TCP cumulative ACK to the sender requesting the first unseen packet in order, like in conventional TCP. The Gaussian elimination may result in a new packet being decoded. In this case the decoder delivers this packet to the receiver TCP. Any ACK generated by the receiver TCP is suppressed and not sent to the sender. These ACKs may be used for managing the decoding buffer. An important point is that the introduction of the new NC layer does not require any change in the basic features of TCP. For more detail, see [START_REF] Sundararajan | Network Coding Meets TCP: Theory and Implementation[END_REF].

III. DYNAMIC MODEL

In this section we develop a model of TCP-NC/RED to study the equilibrium and dynamic properties of the TCP-NC protocol. We start with a non-linear model, study its equilibrium properties in absence of forward delay using optimization theory [START_REF] Kelly | Rate Control for Communication Networks: Shadow Prices, Proportional Fairness and Stability[END_REF], [START_REF] Low | Internet Congestion Control[END_REF], [START_REF] Tang | Equilibrium of Heterogeneous Congestion Control: Optimality and Stability[END_REF], and then linearize the model around its equilibrium to study the effect of the delay and of the redundancy factor ρ on the stability of the system.

A. Network Model

The model introduced in [START_REF] Low | Internet Congestion Control[END_REF] is used to describe the network. A network is represented by a set of |L| = L links, indexed by = 1, . . . , L, with finite capacities c . They are shared by a set of |S| = S sources indexed by i. Each source i uses a set L i ⊆ L of links to convey information through the network to some receiver. The set L i defines an L × S routing matrix R, with R i = 1 if ∈ L i and R i = 0 otherwise. Following the notation in [START_REF] Low | Internet Congestion Control[END_REF], each source i has an associated TCP congestion window size w i (t), a NC redundancy factor ρ i , and an aggregate congestion measure or aggregate price q i (t). Here the redundancy factor ρ i is introduced to account for NC layer. Each link has a congestion measure or price p (t) at time t.

Define the round trip time (RTT) τ i (t) for source i as

τ i (t) = d i + R i b (t) c , (1)
where d i is the round-trip propagation delay assumed to be constant and b (t) is the backlog (queue length) at link . The rate for the i-th source is

x i (t) = w i (t) τ i (t) . (2)
The aggregate input rate at link is the sum of the delayed source rates

y (t) = i R i x i (t -τ f i ), (3)
where τ f i is the forward delay from source i to link . The active queue management (AQM) [START_REF] Low | Internet Congestion Control[END_REF] models describe how the congestion measure p (t) at link is updated. The congestion control problem involves a continuous interaction between TCP algorithms implemented in the sources, and AQM algorithms implemented in the links. The aggregate congestion measure of user i (chosen here as the end-to-end loss probability) is

q i (t) = l R i p (t -τ b i ). (4)
It is the sum of delayed link loss probabilities, where τ b i is the backward delay from link to source i. The RTT time is the sum of the forward and backward delays

τ i (t) = τ f i (t) + τ b i (t). (5)
As in [START_REF] Low | A Duality Model of TCP and Queue Management Algorithms[END_REF], we assume that at each link the congestion measure is modeled through a price mapping function m that maps the queue length b (t) at link to a price p (t). Source i can observe only its own rate x i (t) and its aggregate congestion measure q i (t). Link observes only its own local congestion p (t) and the aggregate flow rate y (t).

TCP algorithms adjust the i-th source rate using some function F i based on q i (t) and ρ i (to account for the presence of the NC layer). The dynamics of the source rate is given by

ẋi (t) = F i (x i (t), q i (t), ρ i ) , (6)
for all sources i. AQM algorithms adjust link prices according to some function G based on the aggregate link rate on link . In the network side, the evolution of the congestion of link is

ṗ (t) = G (p (t), y (t)) . (7)
F i depends on the considered TCP algorithm (Reno, Vegas, etc.) and G i depends on the AQM mechanism (FIFO, REM, RED, etc.) [START_REF] Low | Internet Congestion Control[END_REF], [START_REF] Low | Understanding TCP Vegas: A Duality Model[END_REF]. The interaction between the control at sources and the control at links results in a feedback system, see Figure 1. This network model is used in Section III-B to derive the dynamics of TCP-NC combined with RED-AQM.

B. TCP-NC/RED-AQM Model

To describe the dynamics of TCP-NC, we focus on the congestion avoidance phase of TCP-Reno, in which a large continuous flow set up by TCP spends most of the time. We derive the fluid model approximations from the packet level description of the TCP-NC protocol and study its dynamics. At time t, ρ i x i (t) is the rate at which packets are sent to the receiver. A fraction of these packets is acknowledged at a rate E[R ACKs ] and each ACK increments the size w i (t) of the TCP congestion window by 1/w i (t). Negative acknowledgments, at an average rate of E[R NACKs ], return back to the sender, each one reducing w i (t) by βw i (t). The constant β is the rate at which x i (t) is decreased each time there are losses (typically β = 1/2). The dynamic of the source rate during the congestion avoidance phase of TCP-Reno with NC may be described as follows

ẋi (t) = 1 w i (t) E[R ACKS ] -βw i (t)E[R NACKS ] (8)
Accounting for the effects of the NC layer,

E[R NACKS ] is E[R NACKS ] = x i (t -τ i (t))q i (t)- (ρ i -1)(1 -q i (t))x i (t -τ i (t))- (ρi-1)xi(t-τi(t))-1 k=0 [k -(ρ i -1)x i (t -τ i (t))] × ρx i (t -τ i (t)) k (q i (t)) k (1 -q i (t)) ρixi(t-τi(t))-k (9)
The first term in [START_REF] Shakkottai | Network Optimization and Control[END_REF] is the rate at which packets are lost; this term appears also in TCP-Reno without NC. The second and third terms in [START_REF] Shakkottai | Network Optimization and Control[END_REF] are related to NC. The second is the rate at which NC compensate losses, and the third expresses the fact that x i (t) packets/s are enough to perform network decoding at the receiver. The rate at which ACKs are received is

E[R ACKS ] = ρ i x i (t -τ i (t)) -E[R NACKS ]. In congestion avoidance phase, x i (t)q i (t) ≥ (1 -q i (t))(ρ i -1)x i (t)
, and ( 9) can be approximated as

E[R NACKS ] ≈ x i (t-τ i (t))q i (t)-(ρ i -1)(1-q i (t))x i (t-τ i (t)) (10
) Combining ( 8) and ( 10), the evolution of source rates under the control of TCP-NC (6) may be rewritten as

ẋi (t) = x i (t -τ i (t)) (ρ i -1 + (1 -q i (t)) ρ i ) 1 w i (t) - βx i (t -τ i (t)) (1 -(1 -q i (t)) ρ i ) w i (t). (11)
To study the evolution of the queue length in each link , we consider random early detection (RED) as AQM mechanism. With RED-AQM, b (t) evolves as follows [START_REF] Low | Internet Congestion Control[END_REF] [START_REF] Tang | Equilibrium of Heterogeneous Congestion Control: Optimality and Stability[END_REF] where (a) + z = a if z > 0, and (a) + z = max(0, a) if z = 0. RED averages the instantaneous queue length b (t) by an exponentially weighted average r (t), given by

ḃ (t) = (y (t) -c ) + b (t)
ṙ (t) = -α c (r (t) -b (t)) (13)
for some constant 0 < α < 1. Given r (t), the local congestion measure p (t) at link is given by

p (t) = m (r (t)) =          0, r (t) ≤ b , γ (r (t) -b ), b < r (t) < b , η r (t) -(1 -2p ), b < r (t) < 2b , 1, r t (t) ≥ 2b , (14)

IV. PROPERTIES OF TCP-NC, CONSTANT REDUNDANCY

In this section we characterize some equilibrium and dynamic properties of TCP-NC. To study the equilibrium of a network characteristic (average throughput, loss rate, delay, etc.) under the TCP-NC-AQM control, we will use the fact that the equilibrium of ( 6) and ( 7) can be easily characterized by studying its underlying optimization problem, where we maximize the aggregate utility of the users subject to link capacity constrains [START_REF] Kelly | Rate Control for Communication Networks: Shadow Prices, Proportional Fairness and Stability[END_REF]. To prove the global stability of the algorithm we assume the absence of forward delay, τ f i = 0. The local stability is studied by linearizing the system around the equilibrium as a function of the round-trip time and the redundancy factor.

A. Equilibrium and utility maximization framework

Let (x * i , w * i , q * i , . . . , p * ) denotes the equilibrium point of user i and link . It satisfies

0 = F i (x * i , q * i , ρ i ) 0 = G (p * , y * ) (15)
We consider the utility maximization framework introduced in [5], [START_REF] Low | Understanding TCP Vegas: A Duality Model[END_REF], where the equilibrium source rate vector x * and the equilibrium price vector p * of ( 6)-( 7) are the primal and dual optimal solutions to an equivalent utility maximization problem. A fixed point of ( 6) defines an implicit relationship between x * i and q * i , expressed as

q * i = f i (x * i , ρ i ) (16)
where f i is a positive, strictly monotone decreasing function [START_REF] Low | Internet Congestion Control[END_REF]. This assumption allows to introduce the duality between the optimal solution of a utility maximization problem and the equilibrium of ( 6) and ( 7) by defining U i , the utility function of user i as

∂Ui(x * i ) ∂x * i = f i (x * i ).
Using the assumptions on f i , one deduces that U i has a positive decreasing derivative and is therefore monotone increasing and strictly concave. Now, consider the problem of maximizing the aggregate utility formulated in [START_REF] Kelly | Rate Control for Communication Networks: Shadow Prices, Proportional Fairness and Stability[END_REF] max

x≥0 i U i (x i ) subject to Rx ≤ c (17)
where the constraint says that at each link the aggregate rate y does not exceed the capacity. An optimal rate vector x * exists since the objective function U i is continuous and the feasible solution set is compact [START_REF] Low | Internet Congestion Control[END_REF]. It is unique since U i is strictly concave.

The key to understand the equilibrium of ( 6)-( 7) is to regard the vectors x and p as primal and dual variables, respectively, and (F, G) = (F i , G , i ∈ S, ∈ L) as distributed primal-dual algorithms to solve the primal problem (17) and its Lagrangian dual min

p≥0 i max xi≥0 (U i (x i ) -x i q i ) + p c . ( 18
)
Hence the dual variable is a precise measure of congestion in the network. The dual problem has an optimal solution since the primal problem is feasible. The equilibrium points (x * , p * ) of ( 6) and ( 7) are interpreted as the optimal solution of the primal-dual problem, and that (F, G) iterates on both the primal and dual variables to solve both problems. The rate vector x * is optimal if and only if there exists p * ≥ 0 such that

∂U i (x * i ) ∂x * i = q * i = R i p * (19)
and y * = i R i x * i ≤ c , with equality if p l ≥ 0. The last condition is the complementary slackness conditions. This condition indicates that all bottlenecks are fully utilized at equilibrium or that the AQM should match the input rate to link capacity to maximize the utilization at every bottleneck link. Any AQM algorithm that stabilizes queues possesses this property, and generates a Lagrange multiplier p * that solves the dual problem, as is the case for classical TCP without NC [START_REF] Kelly | Rate Control for Communication Networks: Shadow Prices, Proportional Fairness and Stability[END_REF], [START_REF] Low | Internet Congestion Control[END_REF].

In summary, the utility U i of source i is completely determined by the source dynamics F i , i.e., by the TCP algorithm. The link dynamics G i guarantees complementary slackness condition, and the prices p * play the role of the Lagrange multipliers of the dual problem. The strict concavity of U i guarantees the existence and uniqueness of the optimal rates x * of (17) and the global convergence of the algorithms, and p * is unique when R has full row rank. Thus, ( 6) and ( 7) solve the utility maximization problem given by (17) in a distributed way. Moreover, various TCP-AQM protocols can be modeled as different distributed primal-dual algorithms (F, G) to solve the primal optimization problem (17) and its dual (18).

1) Global stability in absence of forward delay, τ f i = 0: Neglecting the forward delay, from the equilibrium of ( 11) obtained for ẋi = 0, one deduces the average aggregate price seen by user i

q * i = β (ρ i -1) (x * i τ * i ) 2 + 1 ρ i (1 + β(x * i τ * i ) 2 ) = f (x * i ). (20)
The utility function of TCP-NC using ( 19) and ( 20) is

U i (x i ) = 1 -(ρ i -1) √ βτ i ρ i tan -1 βτ i x i + ρ i -1 ρ i x i .
(21) This utility function differs from the one found for TCP-Reno in [START_REF] Kelly | Rate Control for Communication Networks: Shadow Prices, Proportional Fairness and Stability[END_REF] in the amplitude of the inverse tangent and in the addition of a linear term. Both are due to the effect of the NC layer. Nevertheless, the utility for TCP-NC-Reno is still concave and yields a unique equilibrium point. The average throughput is found by considering again the equilibrium of [START_REF] Sundararajan | Network Coding Meets TCP[END_REF], ẋi (t) = 0, leading to

x * i = ρ i -1 + (1 -q * i )ρ i β (1 -(1 -q * i ) ρ i ) 1 τ * i . ( 22
)
To study the global stability of TCP-NC, consider the following function used in [START_REF] Shakkottai | Network Optimization and Control[END_REF] V

(p) = (c -y * ) p + i qi q * i x * i -(U i ) -1 (σ) dσ.
(23) One can show that ( 23) is a Lyapunov function, even when considering NC, V (p) ≥ 0, V (p) ≤ 0, ∀p ≥ 0. By the Lyapunov theorem one deduces that the system is globally stable in the absence of forward delay [START_REF] Khalil | Nonlinear Systems[END_REF].

B. Dynamics of TCP-NC-RED: linearized model

So far we studied the equilibrium point of the TCP-NC model. In this section, the TCP-NC-RED equations are linearized to study the effect of the delay and of the redundancy factor around the equilibrium point. Assuming R has full rank, there is a unique equilibrium (w * i , q * i , τ * i ). Expressing (11) as a function of the congestion window, one obtains

ẇi (t) = w i (t -τ * i ) (ρ i -1 + (1 -q i (t)) ρ i ) τ i (t -τ * i ) 1 w i (t) - βw i (t -τ * i ) (1 -(1 -q i (t)) ρ i ) w i (t) τ i (t -τ * i ) , (24)
where b * denotes the equilibrium queue lengths and

τ * i = d i + R i b *
c the equilibrium RTT. Let δw i (t) = w i (t)w * i . By linearizing (24) around the equilibrium, we obtain the linearized window dynamics

δw i (t) = - 2β (1 -(1 -q * i ) ρ i ) w * i τ * i δw i (t)- ρ 2 i (1 -(1 -q * i ) ρ i ) τ * i δq i (t). ( 25
)
In a similar way, with δp (t) = p (t) -p * and by linearizing (14) around the equilibrium, we obtain the linearized price dynamics

δp (t) = -α c δp (t) + α c γ δb (t). ( 26
)
For the purpose of linearizing the link dynamics of RED, the non-bottleneck links can be ignored. The linearized queue length dynamics around the equilibrium RTT is given in [START_REF] Low | Internet Congestion Control[END_REF]

δb (t) = i R i δw i t -τ f * i τ * i - i k R i R ki w * i (τ * i ) 2 c k δb k t -τ f * i , (27)
Expressing ( 25), (26), and (27) in the Laplace domain, we get

δw(s) = -(sI + D 1 ) -1 D 2 R T b (s)δp(s) (28) δp(s) = (sI + D 3 ) -1 D 4 δb(s) (29) δb(s) = (sI + R f (s)D 5 R T D 6 ) -1 R f (s)D 7 δw(s) (30)
where the D i are diagonal matrices, defined as

D 1 = diag - 2β(1-(1-q * i )ρi)w * i τ * i D 2 = diag ρ 2 i (1-(1-q * i )ρi)τ * i D 3 = diag (α l c l ) D 4 = diag (α l c l γ l ) D 5 = diag w * i (τ * i ) 2 D 6 = diag 1 c l D 7 = diag 1 τ * i
and where R f (s) and R b (s) are the forward and backward delay matrices, with

[R f (s)] i = e -τ f i s if ∈ L i and [R f (s)] i = 0 else, and [R b (s)] i = e -τ b i s if ∈ L i , and [R b (s)] i = 0 alse.
Using (28)-(30), one can derive the loop gain L(s) = δq(s)/δw(s) as

L(s) = R f (s)D 7 (sI + D 1 ) -1 D 2 R T b (s)(sI + D 3 ) -1 × D 4 (sI + R f (s)D 5 R T D 6 ) -1 . (31)
which may be used to characterize stability.

V. VERIFICATION AND STABILITY REGION

In this section we present the properties of the linearized model (31) and the effect of the delay and of the redundancy on the stability region of TCP-NC.

Consider N identical sources sharing a single link |L| = 1 with capacity c, RED as AQM, and equal redundancy factors ρ i = ρ. In this case, dropping the subscripts from (28), the transfer function between the congestion window and the loss congestion measure is

δw(s) = ρ 2 p * z(ρ) e -τ b * s τ * s + 2βz(ρ)p * w * δp(s), (32)
where

z(ρ) = 1-(1-p * )ρ p *
. The equilibrium congestion window is given by w * = τ * c N , where τ * is the equilibrium RTT. The equilibrium price p *

p * = β(ρ -1)(w * ) 2 + 1 ρ (1 + β(w * ) 2 ) ,
is obtained from (20). From ( 28)-( 29), the transfer function between the congestion measure and the congestion window is

δp(s) = αcγ s + αc N e -τ f * s τ * s + e -τ f * s δw(s) (33)
Using ( 30), (32), and (33) the close-loop transfer function of the system is

L N C (s) = αcγ s + αc 1 τ * s + e -τ f * s N ρ 2 e -τ * s z(ρ)p * (τ * s + 2βz(ρ)p * w * ) (34)
The first factor of (34) is due to the queue averaging (13), the second one describes the relation between the congestion window and the buffer size, and the third one illustrates the effect of TCP-NC. Considering a similar scenario without NC, the transfer function for the standard TCP-Reno is obtained from (34) by taking ρ = 1

L T CP (s) = αcγ s + αc 1 τ * s + e -τ f * s N e -τ * s p * (τ * s + 2βp * w * ) (35)
The closed-loop system is stable if and only if the closeloop function (34) does not encircle the point (-1, 0) as s traverses the closed D contour in the complex plane. The pure delay term in (34) adds significant phase at frequencies of the order of 1/τ * and higher, and the loop gain at the crossover frequency will determine stability. To characterize the stability region, we examine the Nyquist plot of The previous simulation scenario is implemented in OPNET using the FTP persistent sessions. The TCP-Reno-NC protocol and RED with explicit congestion notification (ECN) marking are used, as well as the same parameter values for N , c, d, and packet size. For each (N, c) pair, the queue and TCP congestion window trajectories are examined to determine the critical delay d S (N, c) at which the system transits from stability to instability. This delay is measured at 5 ms accuracy. The fundamental frequency of oscillation f S (N, c) is obtained from the FFT of the queue trajectory.

Figure 2 shows the results of our model predictions and of the simulation results obtained from OPNET. Figure 2(a) shows the critical delay d M (N, c) computed using (34), and the critical delay d S (N, c) obtained from OPNET simulations. Each data point corresponds to a particular value of (N, c). The solid lines is the geometric locus where all the points should be if the prediction from the model would perfectly agree with the simulation results. 3(a) shows the size of an individual window and the average window size for the 20 sources, both as a function of time. Oscillations are present in the individual window, and more important the average size shows a deterministic limit cycle, that is also present in the queue trace, Figure 3(b). In this case, the protocol is in unstable regime. This result can be predicted from the model (34) by plotting its Nyquist diagram for the considered scenario. Figure 3(c) shows that the Nyquist diagram encircles the point -1 + j0, leading to an unstable closed-loop system.

Figure 4 shows the stability region for the TCP-NC and TCP-Reno protocols using ρ i = ρ = 1.2 for all users. For each N the critical delay d M (N, c) is plotted versus the capacity c of the common link. Each curve separates the stable region, below the curve, from the unstable one, above it. The TCP-NC curve is dotted, and the TCP-Reno one is solid. TCP-Reno and TCP-NC become unstable when the delay or capacity become large, and the stability region increases with the number of users sharing the link, as with the standard TCP-Reno without NC [START_REF] Low | Internet Congestion Control[END_REF]. The reduction in the stability region of TCP-NC with respect to the standard TCP-Reno is a result of the redundancy factor. Thus, increasing the redundancy factor increases instability.

Figure 5 shows the stability region for TCP-Reno and TCP-NC for redundancy factors ρ i = 1.05, 1.1, 1.15, 1.2 when N = 60 users share the network. Similar to the last experiment, for each pair (ρ i , c) the critical delay d M (N, c) is plotted versus the capacity c of the common link. Taking the stability region of TCP without NC as reference, one can observe a reduction of the size of the stability region for TCP-NC as ρ i increases.

VI. CONCLUSIONS AND FUTURE WORK

In this paper we have used tools from control theory to describe the dynamics of the TCP-NC protocol. We showed that TCP-NC solves a concave utility maximization problem, which yields a unique equilibrium point. We proposed a multilink multi-source model that can be used to study the stability of a general TCP-NC-RED network. The theoretical model we developed and simulation results have shown that TCP-NC becomes unstable when the network scales up in delay and capacity, as the standard TCP-Reno does, but additionally the stability region of TCP-NC reduces, compared to standard TCP-Reno, when the redundancy factor increases. In future work we want to derive stability conditions for the original TCP-NC, and develop a new model to study global and local stability of TCP-NC with an adaptive redundancy factor.

Fig. 1 .

 1 Fig. 1. Interactions between the TCP and the AQM mechanisms

  where b , b , p , are the parameters of the RED algorithm. Morover, γ = p /(b -b ) and η = (1 -p )/b .

  L N C (jω) considering N = 20, 30, . . . , 60 sources, with link capacity c = 8, 9, . . . , 15 packets/ms, and equal propagation delays d = 50, 55, 60, . . . , 100 ms. A constant packet size of 1000 bytes is assumed. For each (N, c) pair, we determine the delay d M (N, c) at which the Nyquist plot intersects the real axis closest to -1. This is the delay at which the system with parameters (N, c) transits from stability to instability. For this delay the critical frequency f M (N, c) at which the phase of L N C (jω) equals -π is also computed. Note that the computation of L N C (jω) requires the equilibrium RTT τ * , which is equal to the sum the of propagation delay d M (N, c) and the equilibrium queuing delay. The queuing delays is computed from the equilibrium model. For all evaluations, the following RED parameters have been used: b = 540 packets, b = 40 packets, α = 10 -4 , and p = 0.1.

  Figure 2(b) represents the corresponding critical frequencies f S (N, c) and f M (N, c). It can be observed from Figure 2 that our predictions are acceptable, since the resolution for the delay is 5 ms.

Fig. 2 .

 2 Fig. 2. Model vs OPNET simulation

Figure

  Figures 3(a) and3(b) show the OPNET simulation results corresponding to the trajectories of the sizes of the congestion window and of the queue when N = 20 sources share a link of capacity c = 9 packets/ms, with a round trip propagation delay d = 50 ms. Figure3(a)shows the size of an individual window and the average window size for the 20 sources, both as a function of time. Oscillations are present in the individual window, and more important the average size shows a deterministic limit cycle, that is also present in the queue trace, Figure3(b). In this case, the protocol is in unstable regime. This result can be predicted from the model (34) by plotting its Nyquist diagram for the considered scenario. Figure3(c)shows that the Nyquist diagram encircles the point -1 + j0, leading to an unstable closed-loop system.Figure4shows the stability region for the TCP-NC and TCP-Reno protocols using ρ i = ρ = 1.2 for all users. For each N the critical delay d M (N, c) is plotted versus the capacity c of the common link. Each curve separates the stable

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Congestion window trajectory, queue trajectory, and Nyquist stability diagram for N = 20, c = 9 packets/ms and d = 50 ms


Christophe Sibertin-Blanc

Olivier Therond
email: olivier.therond@toulouse.inra.fr

Claude Monteil
email: monteil@ensat.fr

Pierre Mazzega
email: pierre.mazzega@lmtg.obs-mip.fr

Formal Modeling of Social-Ecological Systems

The success of Integrated Assessment and Modeling of social-ecological systems requires a framework allowing the members of such a process to share and gather their respective knowledge about the system under consideration and to get confidence into the reliability of the software that implements the system's model they have produced. To this end, this paper presents an Entity-Process meta-model of SESs and outlines its use.

Introduction

Social-ecological systems (SES) are considered as complex adaptive systems that are characterised by self-organisation and distributed control. They encompass multiple actors with diverse and contrasting management interests and objectives, acting at different spatial and temporal levels [START_REF] Reed | Stakeholder participation for environmental management: A literature review[END_REF][START_REF] Pahl-Wostl | Transitions towards adaptive management of water facing climate and global change[END_REF][START_REF] Giampietro | Complexity and Scales: The Challenge for Integrated Assessment[END_REF]. In such systems, social, economic and ecological processes, in interaction at many organization levels, give rises to emergent structures and functions at the whole SES level.

During the last decades integrated assessment (IA) has been playing an increasing role to better understand how complex systems function. [START_REF] Rotmans | Integrated assessment: A growing child on its way to maturity[END_REF] define integrated assessment as an interdisciplinary and participatory process to combine, interpret and communicate knowledge from diverse scientific disciplines. IA approaches to assess system responses to scenarios of change have been commonly used to design sustainable management and development strategies. Mainstream in IA is based on modeling and simulation. This Integrated Assessment and Modeling (IAM) approach aims at describing quantitatively the causal relationships and the interactions between the various components of the investigated system. According to the problem at hand, IAM approach assembles the relevant information from a wide range of scientific disciplines and put them into a "policy oriented context" in order to analyze complex system behaviour [START_REF] Toth | Integrated environmental assessment methods: Evolution and applications[END_REF]. The potential of model-based methods is well-established for handling management and policy problems [START_REF] Schreinemachers | An agent-based simulation model of humaneenvironment interactions in agricultural systems[END_REF][START_REF] Mcintosh | Tools to think with? Towards understanding the use of computer-based support tools in policy relevant research[END_REF][START_REF] Jakeman | Ten iterative steps in development and evaluation of environmental models[END_REF][START_REF] Oxley | Integrated modelling and decision-support tools: a Mediterranean example[END_REF].

When applying IAM to social-ecological systems, one of the main challenges lies in the integration of many different disciplinary viewpoints and forms of knowledge. One fruitful way to reach this objective consists in building up a well-funded and shared ontology of the SES, i.e. a formal specification of the concepts and relationships among these concepts [START_REF] Beck | Ontology-based simulation in agricultural systems modeling[END_REF].

Another challenge is the confidence of the participants of the modeling process in the simulation results. This confidence is problematic if, once they have built the model together, they lose the control on the following of the process, that is the implementation by the computer scientists of the model as a piece of software. The use of a meta-model allows to produce a well-structured model whose implementation suffers from no technical distortion.

This paper presents an ontology expressed as a meta-model, intended to provide scientists and stakeholders from various competencies and backgrounds with a generic and formal framework to coherently arrange their respective knowledge about SESs. It guides collective modeling processes aiming to assess SESs' scenarios by using simulation models. Examples of models and corresponding specific diagrams (with a well-defined semantics) used for their graphical representation are presented in another paper [START_REF] Therond | A Meta-Model to support Integrated Modeling of Social-Ecological Systems[END_REF].

An Entity-Process Meta-model for SESs Modelling

Working in an interdisciplinary environment requires a shared conceptual framework for a collective interpretation of the concepts put forward by one another. Accordingly, our meta-model uses terms understandable and usable both by stakeholders and by a large scientific community, rather than referring to domain specific terminologies. The primary function of the meta-model is to provide a formal framework facilitating the coherent description of SESs as formal models. The meta-model identifies the types of the constitutive elements that are considered in the design of SESs' models and how they are related, as shown in Fig. 1 represented as an UML class diagram. According to this meta-model, the model of a SES includes a set of entities and a set of processes, together with relations between elements of these sets. The structure of a system is composed of entities and relations between entities, both being perceivable and more or less persistent. The dynamic character of SESs is handled through processes that involve one or more entities of the system, and make the state and also the structure of the system to evolve. Finally, we account for interactions between a SES and its environment.

System structure: entities and relations

The meta-model distinguishes three categories of entities: material resources, cognitive resources and actors. Instances of entities are the considered individual resources and actors that may appear or disappear over time, or experience changes of their state. Entities are characterized by properties whose values constitute the state of an entity instance. They are endowed with operations that process these property values.

There are two types of relation between entities. A structural relation permanently associates entities related by their very nature (e.g. field -soil type). A non-structural relation puts in association entities as the result of some action (e.g. farmer ([wner of an acquired] field). A link is an instance of a relation; it ties instances of the related entities. Actors are human agents, be they an individual, a population of similar individuals or a group such as an organization, association, committee, enterprise, etc. Each actor executes at least one activity.

Resources are objects of uses, or objects for uses. Material resources are physical resources spatially and temporally distributed (e.g. a water body, a field plot, a crop). Cognitive resources stay in the minds of human beings: they are information, believes or expectations about facts, procedures, values or aims, that humans use or consider in the activities they undertake, in designing strategies or in the formation of goals. A cognitive resource is for a single actor (e.g. the cropping system of a farmer) or shared by several actors (e.g. the weather forecasts). Many cognitive resources have a material counterpart -for example a legal norm (cognitive resource) has an official publication support (material resource) -which is most often not relevant to consider. It should be stressed that any knowledge about a material resource used by an actor or a process implies the existence of an own cognitive resource; such a cognitive resource will be considered in the model only if required (e.g. if its proprieties do not correspond to those of the material resource). Every resource must be either impacted or used by at least one process (this statement is the main integrity constraint of the meta-model, see 3.1 below).

At each time the structure of a system is defined as the sets of its entity instances and of its structural links (the instance of the structural relations). At each time the system's state is defined as the states of these entities (i.e. the values of their properties) together with the existing non-structural links.

System dynamics: processes and behaviors

The processes generate the observed phenomena and make the structure or the state of the system to change. Each execution of a process impacts entities by changing their property values, or by creating or destroying entity instances or links. We distinguish three types of processes: activities, socio-economic processes and ecological processes.

An Ecological process corresponds to an enactment of a biophysical law. It determines the evolution of the states of material resources (and possibly of actors).

A Socio-economic process generates phenomena resulting from the human beings' economic or social activities. Only the effects of a socio-economic process is described in the model, considering that it is not necessary to represent the actors' activities at the origin of this phenomenon, that is how it is generated.

An Activity is a process executed by an actor intending to achieve some goals. A goal is linked to the activity aiming at its achievement.

Each actor is associated with one behavior that describes how he orchestrates his various activities, solves the conflict that could appear between them, and decides to execute (or not) at a given time an action rather than another one.

The description of a process is broken down into actions that can themselves be decomposed into sub-actions. The dynamics of a process describes the scheduling of actions executions as a system of equations, an automaton, an algorithm, or a set of (event-condition-action) rules, according to the most appropriate representation for the considered process. A process dynamics uses entities (resources or actors) whose states condition its development.

Processes interact by means of entities. Indeed, it is through changes in the states of the impacted entities that a process affects the course of the processes that depend on these entities.

Fig. 2:

The interactions between a SES and its environment (outside the dashed-dotted line box).

Interactions between the system and its environment

Social-ecological systems are open systems that interact with their environment, that is the entities and processes of the reality that are outside the system while having an influence on it. The meta-model considers these interactions through the setting of resources at the interface of the system. We call external process a process whose execution is not affected by the state of the modeled system and which changes the state of some entities that are, from this fact, qualified as interface entities (see Fig. 2). An interface entity may be impacted by both external and internal processes and it is similar in nature to other entities. Classically the dynamics of an external process is not represented but rather its impacts. These impacts are commonly framed as a (spatio-) temporal series, representing for each time step of the simulation, the new values (or the variations) of interface entities' properties (e.g. climate data).

We call external actor an actor who performs only external activities. Outputs of the system to the environment are considered only whether they constitute indicators to be considered.

The Modeling Process

Any interdisciplinary modeling of a SES requires, in addition to a conceptual framework such the one presented in the previous section, a methodological framework to structure and support the modeling process. The structure of the model to be built guides the structure of the modeling process by defining the objects it must produce. We here mainly address the main issues related to a tractable representation of models conforming with our meta-model, in order to illustrate its practical usability.

The meta-model described in section 2 is a tool for supporting integrated modeling of a given reference system. By "reference system" we mean the collective representation that observers build of the knowledge they share about the reality they investigate regarding a given question (see Fig. 3). The first step of the modeling process consists in an instantiation of the meta-model into a conceptual model that identifies relevant entities (actors, material resources and cognitive resources), relations and processes when considering the question the model is intended to answer.

Practically, several conceptual models are produced in the course of a progressive and iterative modeling process, along with the improvement of the understanding of the question at hand. They differ by their levels of granularity or of generality, their focus on the core of the system, or their boundary. The conceptual model to be simulated must be internally consistent in terms of (spatial, temporal, sociological, physical, etc.) granularity and ontological level of the represented entities and processes.

The structural part of the conceptual model is in turn to be instantiated into a concrete model of the reference system. Indeed, the entities and relations of the conceptual model are just the types (or classes) of the reference system's entities and relations to be considered. For example, if the conceptual model includes a "Water Agency" actor, the concrete model will include a single instance of this type whose properties values describe the particularities of the Water Agency at hand; a cognitive resource "Weather forecast" will be instantiated as many different instances (one for each time periods of the simulation); a material resource "Field plot" will be instantiated into the number of concrete plots necessary to cover the agricultural territory of the reference system. The same holds for the relations between entities: a relation "owns" between an actor "Farmer" and a resource "Field plot" in the conceptual model will give rise in the concrete model to concrete owns links between each Farmer instance and the field plot instances he posses.

Finally, the concrete model will be implemented by the development of a computer program, the simulation model.

The conceptual model of a SES includes too many elements to be represented at one single glance. Thus, its representation must be broken down into several diagrams, each one focused on an aspect and using specific notations. All these diagrams must be documented with the explanations necessary for their unambiguous understanding by all the participants of the modeling process. Indeed, the comprehensive view of the system's representation that is embedded into the model is necessary for the correct interpretation of the simulation results, and thus for the achievement of the modeling process's goals. All the diagrams regarding the system's structure and dynamics must be carefully aligned to ensure the global coherence of the system's conceptual model. Accordingly, this alignment is obtained through back and forth steps between the various diagrams.

Fig. 3:

The different models and their relations to the reference system.

Actor-resource diagrams and integrity constraints

The actor-resource diagram shows the actors and the resources of the conceptual model together with their relations, that is the structural part of the model.

It may be represented by means of the widely used UML notation for class diagrams [START_REF] Omg | Object Management Group, UML 2.0 Specification[END_REF], but it could be represented using the Entity-Relationship formalism as well (Chen et Pin-Shan, 1976). To distinguish the actors, material resources, cognitive resources and interface entities, their surrounding boxes may have different colors or shape (see Plate 1 (a) for the conventions we use).

It is likely that the readability constraints of the actor-resource diagram do not agree with a layout on a single page. Helpfully the diagram can be split into several subdiagrams, leading to the identification of sub-domains of the system's model.

An actor-resource diagram must be accompanied with the statement of its integrity constraints. Integrity constraints are restrictions on the instantiation of entities and relations that ensure the possibility to interpret the concrete model as a description of a coherent and feasible world. Cardinalities of the relations express some of such constraints, e.g. in the UML class diagram of Fig. 1, each process has to impact at least one entity. Other constraints deserve to be explicitly stated, which bear on the values of properties of an entity or of linked entities, or on the existence of links or entities. For example, in a conceptual model including actors "Farmer", the value of the property "income" of any farmer must be the sum of the incomes of the fields crop he works.

It is often the case that a system' model requires to refer to different theories that represent a same item of the reality by different entities [START_REF] Quesnel | VLE: a multimodelling and simulation environment[END_REF]. For example, the representation of water flow and plant growth processes often needs to represent the investigated area by resources that differ with regard to their properties or spatial delineation. Whether the area is represented through different resources by the two processes, the state of resources used and impacted by the plant growth process must take into account the effect of the water flow process and conversely. Integrity constraints allow specifying how the relations between the states of the two types of resources have to be handled. The UML proposes to formally express these constraints using the Object Constraints Language (OMG, 2010) but any means for unambiguous statement may be used. (e) Actor A performs an activity that impact the resource R. (f)

Actor A uses the resource R to execute at least one of his activities; the relation can be labeled or not. (g)

Actor A performs an activity that simultaneously impacts resources R1 and R2, their respective states being inter-dependent.

(h) Actor A performs an activity that changes the realization of a process being dependent on the state of resource R1 and impacting resource R2.

Plate 1: Graphical conventions for entities are given in panel (a). In panels (b) to (h) resources can be indifferently material or cognitive.

Process diagrams

The process diagram displays the processes and the entities they impact, the actors performing the activities and the most relevant entities used by their dynamics. For readability, it may be broken down into several sub-diagrams. To get a synthetic representation, Plate 1 suggests graphical conventions for basic patterns.

To be illustrative, Fig. 4 shows the process diagram of the quantitative water management in a French association of irrigators. In such associations the water withdrawals are generally carried out through collective equipment (pump) and water is distributed to irrigators' fields by water pipes. Association often withdraws water in several water resources like reservoir hills and rivers. Water withdrawals in river are generally limited by a water quota attributed to the association. The association defines its management rules to distribute available water resources among association's members. The president of the association is generally one of the irrigators and often plays the role of coordinator and regulator. The process diagram of Fig. 4 presents is aimed at simulating impacts of management rules scenarios on different indicators regarding water withdrawals (quantity and dynamics), nitrogen pollution and farmer income. The considered time horizon of scenarios is a short term one (about 15 years). In the model presented in Fig. 4 two types of actors perform activities: irrigators and the association president (AP). "Cooperative group" and "Irrigator association" are external actors. They are represented to ensure that modelers share understanding of the nature and origin of interface resources. Irrigators execute two activities: "allocate crop acreage" and "irrigate". The AP checks the compliance of irrigators' behaviors with the association rules and if necessary puts fine to irrigators to regulate their behavior. The resource "Water level" is generated by the AP's activity "measure". Both the AP and irrigators use "Association rules" as cognitive resources. The definition of association is provided by external process scenarios. The relations between the cognitive resource "Water level" and the material resources "River water" and "Hill reservoir water" express integrity constraints between these entities: "Water level" directly corresponds to "River water" and "Hill reservoir water". To execute the activity "allocate crop acreage" irrigators use "Market price", "Expected CAP premium" as resources. To perform "irrigate" they use the cognitive resource "Weather forecast".

Only one internal ecological process is considered in this model: "soil water and plant dynamics" (SW&P dynamics). It allows simulating crop growth, water and nitrogen cycle and crop yield in each irrigated field plot. The graphical features of plate 1 allow representing that the activity "irrigate" performs a water transfer from rivers or hill reservoirs to soil-crop resources (i.e. plot field).

Only one internal socio-economic process is considered in this model: irrigation technology dissemination. It allows representing evolution of agricultural irrigation practices over the temporal extent of scenarios without considering by who and how it is performed.

Dynamics and Behavior diagrams

A dynamics diagram describes the dynamics of one process, so that there is a dynamics diagram for each internal process. It describes the actions of the process and which ones must be performed at each time step of the simulation. As there are many ways to formalize the dynamics of a process, no standard representation is proposed. What about the multiplicity of the entities instances of the concrete model impacted by the execution of a process at each time step of the simulation? The better is to consider that each process occurs only once and applies to every (concrete) instance of the impacted entities. As an example in the model shown in Fig. 4, every occurrence of the Soil Water & Plant dynamics process applies to each instance of the Soil crop resource. The correct time and spatial resolutions of the model should be the finest one making sense for all the processes.

The processes interact through the entities of the model's structure. To ensure that they act within the same structure, there must be no ambiguity about the identities and the properties of the entities involved in their respective executions. They must also specify the entities' operations they need for their processing, in order to complete the actor-resource diagram. A dynamics diagram is to be precise enough so that it can be implemented as a computer routine without ambiguity.

An actor's behavior diagram describes how an actor having several activities to execute selects, among the actions of these activities enabled by the current system's state, the ones to be effectively executed. The various activities of an actor should be described using the same formalism, which can also be used for the description of the actor's behavior diagram.

Scenarios and indicators

Scenarios intend to study the effect upon the behavior of the system of various hypothesizes. They are expressed as variations of the concrete model. They can encompass changes into (i) the dynamics of internal processes, notably their parameters (ii) the impact of external processes upon interface entities, or (iii) the initial states of entities and links of the concrete model. In our example on water management (see Fig. 4), the scenarios of rainfall over several decades are introduced in the form of spatiotemporal series on the hydrological.

The integrated assessment of scenarios is achieved by analyzing the value of indicators derived from the final value or from the evolution, in the course of the simulation, of the state of entities of interest for scientists and stakeholders.

Documentation of the concrete model and scenarios

The system's concrete model describes the initial structure and state of the system to be simulated that is the instances of the entities and their state, together with their links. It includes too many elements to be displayed as diagrams. A comprehensive documentation of the concrete model must describe: 1. the concrete instances of the model's entities, including the values of their properties and their links, together with the sources of the data used to define these instances; if the instances of some entities are very numerous, they are not described individually but by a probabilistic law that distributes the values of their properties; 2. the sources of the data that are used to feed the external processes and how they are processed, if necessary, to fit the structure of the interface entities; 3. an implementable specification of the dynamics and the actions of each process; 4. the algorithms allowing to calculate the values of indicators.

For each scenario, the changes have to be documented in the same way.

Simulation model

The simulation model is an implementation of the concrete model (or one of the scenarios) as a software3 . Classically in IAM approaches the simulation model looks like a black box for most of the participants to the modeling process [START_REF] Resnick | Beyond Black Boxes: Bringing Transparency and Aesthetics Back to Scientific Investigation[END_REF]. They do not know to what extent the simulation model is a faithful implementation of the concrete one. This lack of comprehensive view of the simulation model often leads to little or even no confidence in the simulation results.

Using the entity-process meta-model, the structure of the simulation model can mirror the one of the concrete model and so prevent the occurrence of distortions between the two models. Using an object-oriented language such as Java, each entity of the conceptual model is implemented as an object class provided with the corresponding attributes and methods. The relations are implemented as entities' attributes, according to the navigability 4 of the relations that can be deduced from the processes' actions. These classes will be instantiated as objects according to the concrete model. Each process will be implemented as a class including class attributes and a method that implements its dynamics. Each actor class will be in addition provided with a method that implements its behavior.

The simulation model can be developed using a simulation platform providing useful services such as Cormas (Bousquet et al. 1998), Mimosa (Müller, 2010), Record (Bergez et al., submitted), Gamma (Taillandier et al., 2010), Repast (North et al. 2006) or many others -see the simulation platforms reviewed in [START_REF] Schreinemachers | An agent-based simulation model of humaneenvironment interactions in agricultural systems[END_REF].

Thanks to the structure of the behavior of the model, as a set of processes that interact in an asynchronous way by means of entities, the engine that makes the simulation to progress is quite simple. Assuming that the time scales of all the processes are nested, at each time step it has to carry out the following: 1. for every internal ecological and socio-economic process, activate its dynamics to select the actions to be performed; 2. for each actor, activate its behavior to select the actions of each of his activities that are enabled, and then, according to his behavior, select the action(s) to be performed; 3. perform the update of the state of the interface entities (as a result of the effect of external processes); 4. perform the actions selected by each internal bio-physical process, by each socioeconomic process and then by each actor's behavior; in case of conflicts between some actions (e.g. the execution of one action disabled a selected other one), the conflict must be solved either at random or by following well-specified rules.

A simulation corresponds to a numerical experiment performed by applying the simulation engine to the simulation model of a given scenario. Each numerical experiment allows assessing indicators. Comparisons of indicator values of a given scenario to the ones of a reference, either the current situation or a reference scenario, allow the integrated assessment of scenarios.

Conclusion

The concepts proposed by the Entity-Process meta-model are elementary, while the proper modeling of complex SESs requires recourse to much more elaborated concepts such as "territory", "public resource" or "policy". We believe this meta-model can be specialized for specific modeling purpose (regarding the type of the investigated reference SES and investigated question) by providing rigorous definitions of higher level concepts, built upon concepts of this meta-model. This is nothing but the "profile" feature of the UML.

The meta-model presented in the paper can be used for two different purposes: to support a priori the Integrated Assessment and Modeling process of a SES or to make explicit a posteriori the conceptual model of an already implemented model or of a quite informal model that does not satisfactorily manage the complexity of the system at hand. In both cases, the formal nature of the meta-model allows a quite easy comparison of models and of their components [START_REF] Bezivin | On the Unification Power of Models[END_REF].

Fig. 1 :

 1 Fig. 1: A meta-model for the formal modeling of social-ecological systems represented as a UML class diagram. Cardinalities use the following conventions: *: any number of elements (including none) ; +: at least one element ; n: exactly n elements.

  socio-economic process being dependent on the state of the resource it impacts.

Fig. 4 :

 4 Fig. 4: A process diagram of water management at the scale of an irrigator association in France. SW&P means "Soil Water and Plant". See Plate 1 for the graphical conventions. Resources positioned on the dashed-dotted boundary line are interface resources.

According to the terminology of logic, the concrete model would be considered as a theory, and the SES and the simulation model as two models (or interpretations) of this theory in the worlds of perceivable reality and software respectively.

The navigability of a relation between entities specifies which of the related entities needs to access to the other.

Acknowledgements. This work is funded by the Thematic Network for Advanced Research «Sciences & Technologies for Aeronautics and Space» (RTRA STAE http://fondation-stae.net/) in Toulouse, France, under the MAELIA project.


Fabien Caleyron

Alain Combescure

Vincent Faucher

S Potapov

SPH modeling of fluid-solid interaction for dynamic failure analysis of fluid-filled thin shells

  L'archive ouverte pluridisciplinaire


Sushma Patil

Eitan Altman
email: eitan.altman@sophia.inria.fr

Manjesh K Hanawal
email: shanawal@sophia.inria.fr

Julio Rojas-Mora

Modeling and Simulation of Mobility of Crowds

Keywords: Mobility, flocking, statistical analysis

Mobility models studied in the networking community usually assume independence between the movement of individuals. While this may well model sparse networks, there are many scenarios that might not follow this assumption. In contrast, within other communities, such as road traffic engineering, biology and computer graphics, models of mobility usually take into account the dependence of the mobility pattern of an individual with respect to that of its neighbors. Our goal in this paper is to study how this dependence impacts the performance measure from the networking point of view. In particular, we implement a bio-inspired model for mobility of crowds and, by simulation, we study how mobility influences the performance measures of a distributed network. We perform statistical analysis on the samples obtained through simulations. In particular, we study the distribution of the message delivery time and show that it is light tailed, with exponential tail distribution.

I. INTRODUCTION

Mobility of nodes in a wireless network have a direct influence on the performance of various protocols used for message delivery; performance metrics like delay in message delivery, energy consumed, probability of successful message transmission, etc., are functions of mobility pattern of nodes. It is desirable to know the performance of a protocol by testing it on real world traces of the node mobility patterns collected through experiments. However, often many a times traces of the desired scenario are not available as the wireless network may not yet be deployed. In such scenario, one relies on the data generated from synthetic models to study network performance. Synthetic models that describe a mobility scenario under consideration, in a realistic fashion, can be complex. The models available in literature tradeoff between ease of implementation/analysis and being more realistic. Depending on which synthetic model best captures node mobility, the choice of the synthetic model is made for evaluation of the performance of a protocol.

Quite often, simple mobility models like Random Waypoint or its variants [START_REF] Hyytia | Spatial Node Distribution of the Random Waypoint Mobility Model with Applications[END_REF] are assumed as they are relatively easy to implement and analyze [START_REF] Bettstetter | The Node Distribution of the Random Waypoint Mobility Model for Wireless Ad Hoc Networks[END_REF]. These models assume that movement of each mobile node (MN) is independent of others, and also their movement does not depend on their past locations or speeds. These assumptions are justifiable in a very sparse network where MNs are spread far from each other and hence their movement can be assumed to be uncorrelated. These models, though easy to analyze, lead to unrealistic scenarios, like sharp turns or sudden stops. To avoid such cases models are proposed that correlate a node's mobility on its past movements to varied degrees of generality, but with increasing complexity of tractability. In Gauss-Markov models [START_REF] Liang | Predictive distance-based mobility management for PCS networks[END_REF], current location of a node is made to depend on its past location through a controlled parameter. For very thorough surveys on mobility models see [START_REF] Camp | A Survey of Mobility Models for Ad Hoc Network Research[END_REF] and [START_REF] Ranjan | Handbook of Mobile Ad Hoc Networks for Mobility Models[END_REF].

In many scenarios MNs may not move independent of each other, but in a group, or following a particular reference node in the group. For example, a group of people visiting a museum or touring a city following an instructor. To cater to this kind of dependent mobility, several models are proposed by the network community. In the Reference Point Group Mobility (RPGM) model all of the nodes follow a path traveled by a logical center [START_REF] Hong | A group mobility model for ad hoc wireless networks[END_REF]. The logical center may be made to follow a predefined path or can follow random waypoint mobility. The other MNs move according to some random waypoint in the vicinity of the logical center. Several other useful group mobility models can be derived as variants of RPGM: column group mobility, nomadic group mobility, pursue mobility [START_REF] Camp | A Survey of Mobility Models for Ad Hoc Network Research[END_REF], etc.

We shall be interested in crowd movements where MNs are humans carrying wireless terminals. Here the MNs may not follow any particular MN, nor their movement is independent of each other. In such crowd movement, they avoid colliding with each other and maintain a safe distance from their neighbors. For example, people moving in a market or a busy commercial area. In this paper, we aim to study such crowd mobility where mobility of a MN depends on its neighbors. In biology, such kind of mobility is extensively studied to understand movement of flock of birds, school of fish or heard of animals. Rules governing such mobility are extremely interesting if one wishes to recreate patterns of such flock movements through artificial life. Indeed, scientists in the field of computer graphics synthesize mobility of such crowds to create beautiful patterns [1]. To synthesize such patterns they implement rules that govern how mobility of an entity in the group depends on its neighbors.

Craig Reynolds [START_REF] Reynolds | Herd, and Schools: A Distributed Behavioral Model[END_REF] proposed a systematic way to synthesize aggregate motion of crowds using computer graphics. He referred to each entity in the crowd as 'boid'. The synthesized movement of the boids is based on his finding that the complex auto-organization of the group into complex macroscopic patterns is determined by three simple microscopic rules that each individual in the group follows: (i) Boids try to fly towards the center of mass of neighboring boids. (ii) Boids try to keep a small distance away from other objects (including other boids). (iii) Boids try to match velocity with near boids. If each boid is implemented as an independent actor, which navigates according to its local perception of the dynamic environment obeying the above rules, the interaction of the boids results in the aggregate action that looks very synchronized as if centrally controlled.

The primary aim of [START_REF] Reynolds | Herd, and Schools: A Distributed Behavioral Model[END_REF] is to recreate beautiful patterns observed in nature. Variants of the rules proposed by Reynolds are used to study local interaction and evolutionary patterns in biology [1]. In one of such studies, Ariel Dolan [2] implemented creatures that guard their territory against intruders, and study evolutions of populations. He referred to such lifelike creature as 'Floy', which interact with their neighbors as follows: (i) stay close to your fellows but not too close. (ii) if you see an intruder, move towards it and attack. The first rule of Dolan is similar to the first two rules of Reynolds, but they differ in how the Floys track its neighbors: In Dolan's Floy model, each Floy tracks two of its randomly chosen neighbors and tries to be close to them, whereas in Reynolds' model, each boid tries to be in the center of its neighbors that are within a certain radius. Also, in Dolan's model there is no matching of velocities among the Floys.

In Reynolds' boids model the three rules make the mobility of boids highly dependent on each other; the boids stay close to each other, and result in patterns that give a feeling that all the boids are guided centrally. However, relaxed rules of Dolan's Floys model makes Floys to avoid too much togetherness, which lead them to move in an ensemble that spreads over a wide region. Each Floy may stray away from its group and wander somewhat randomly in a given territory around its group members. But it will not go too far away, as soon it will tend to move towards its group due to the first rule. This behavior is similar to that of a crowd movement where the nodes move from one location to other. Obviously, MNs avoid colliding with each other and do not stray far from its neighboring fellows while it moves in a given territory. Thus, we can use rules similar to that in the Floy model to simulate and study crowd movements in ad hoc networks.

Our aim in this paper is to study performance in our crowd mobility model, from the networking point of view. We consider the delay tolerant network scenario. Each MN in the crowd can have a message that another MN in the crowd is interested. The MN which have this message (source MN) can spread the message in the crowd by giving it to any other MN that comes close to it. The message is thus relayed among the crowd and can finally reach the MN that is interested in this message (destination MN). We shall study the time taken for a message sent by a source MN, using the other MNs as relays, to reach the destination MN. In particular, we shall be interested in identifying the tail distribution of the message transfer delay and its statistics, i.e., the probability that the message takes larger than a threshold time to reach destination. Such analysis is of primary importance as a message may become irrelevant if it does not reach the destination within a stipulated time. We also study how the transmission distance influences the message delivery time. This is equivalent to analyzing energy-delay trade off in wireless networks [START_REF] Small | Resource and Performance Tradeoffs in Delay-Tolerant Wireless Networks[END_REF], as larger transmission range means more power transmission with lesser delay and vice versa.

The paper is organized as follows: In Section II, we discuss the simulator we built to study crowd mobility based on Craig Reynolds' boid model. In Section III we discuss how we adopt the boid model to study network crowd movement and explain the simulation settings. In Section IV we discuss the statistical results of the message delivery time. In Section V we discuss the effect of transmission range on the message delivery time. Finally, in Section VI we discuss the possible extensions and summary of our work.

II. CROWD MOBILITY SIMULATOR

In this section we discuss the simulator we built to study crowd mobility based on Reynolds' boid model. We refer to each boid as a node, henceforth.

Below we discuss how each node updates its location and velocity according to Reynolds' three rules discussed in the introduction. In the 2-D version of the simulator each node begins with a fixed velocity v 0 at an arbitrary location in a rectangular region. In applying rule 1 (cohesion), each node looks for other nodes within a radius, say R 1 , and finds the position that is the center of mass with respect to the location of these nodes. In applying rule 2 (separation), each nodes looks for other nodes within a radius, say R 2 , and finds the position that is the center of mass with respect to the the locations of these nodes. In applying rule 3 (alignment), each node looks for relative velocity of the other nodes within a distance, say R 3 , with respect to itself, and finds the average velocity. The new velocity is obtained by combining the values obtained from the three rules with weights γ 1 , γ 2 , and γ 3 , respectively, and adding to the current velocity. For example, let A 1 = (a 1 , b 1 ), A 2 = (a 2 , b 2 ), and A 3 = (a 3 , b 3 ) be the values obtained from rule 1,2, and 3, respectively, for a node which has current velocity V n = (V x n , V y n ) and is at position

X n = (x n , y n ).
Here the first component in each vector corresponds to the x-axis, and the second to the y-axis. Then, the node moves to the new position X n+1 = (x n+1 , y n+1 ) with velocity V n+1 = (V x n+1 , V y n+1 ), given by

V n+1 = V n + γ 1 A 1 -γ 2 A 2 + γ 3 A 3 and X n+1 = X n + ∆ * V n+1 ,
where ∆ is a constant that governs how fast the nodes update their position. Note that the weights γ 1 , γ 2 , and γ 3 govern the relative importance assigned to each rule. By adjusting these parameters one gets varied degree of dependency among the nodes' movement. For example, if γ 1 is large while γ 2 and γ 3 are small, then nodes stay very close to each other, but seem to change directions and come close to each other often. If γ 2 is large while γ 1 and γ 3 are small, then nodes are equally spaced from each other but move more randomly and change direction often. If γ 3 is large while γ 1 and γ 2 are small, then all the nodes align themselves in a particular direction and continue to move in the same direction till they hit the boundaries. By varying these weights one gets varied degree of randomness in the movements of the nodes. We put an upper limit on the velocity of each node. If the velocity of a node exceeds this limit, denoted by v max , then its velocity is set to v max .

We implemented Reynolds' three rules in Java with an interface to control the weights. A snapshot of our implementation is shown in figures 1 and 2. In the simulator, the number of nodes (N ), and the parameters R 1 , R 2 , and R 3 can be set. If the values of R 1 , R 2 , and R 3 are high, i.e., nodes take into account movement of other nodes from a larger area in deciding their new velocity, and the movement becomes more dependent. When the values of these parameters are decreased, the movement tends to become more random. Thus, by changing various parameters in the simulator, the movement can be changed from being completely random to perfectly synchronized.

In building the aforementioned simulator, we used the Dolan's Java templates which he developed to implement his Floy model. We implemented Reynolds' boid with the help of these templates and built in various functionalities to control movements of boids that help us study crowd mobility in ad hoc networks as discussed in the following sections. The simulator and its source code are available online at [17].

III. AN ADAPTATION OF BOIDS MODEL TO NETWORKING

In crowd movement, people try to avoid colliding with each other by maintaining safe distance, while also not straying far away from its neighbors. This scenario can be simulated by appropriately setting the weights γ 1 , γ 2 , and γ 3 and the parameters R 1 , R 2 , and R 3 in the simulator. Note that in crowd mobility, velocity matching is not an important criterion. However, avoiding collisions and not straying away far from neighbors are important conditions, with the former having more priority. These criteria correspond to setting γ 3 = 0, and γ 2 > γ 1 in the simulator.

We adopt the boid model to study crowd movement in wireless ad hoc networks. In particular, we will be interested in scenarios like delay tolerant networks (DTNs) where endto-end connectivity between source and destination is not available. Everybody in the crowd is assumed to have a transmitter and receiver, and can transmit, receive and relay messages. We henceforth refer to each boid-like entity or each person in the crowd as a "mobile node" (MN). The MNs are battery powered and can not transmit beyond a certain transmission range. Let T denote the transmission range of each MN. They can transmit or receive messages when within distance T from other MNs. If the other MNs already have this message, then they simply ignore the message.

In this preliminary study, we consider epidemic routing scenario in a delay tolerant network [START_REF] Altman | Optimal monotone forwarding policies in delay tolerant mobile ad-hoc networks[END_REF], i.e., a MN having the message transmits it to all the MNs that are within the transmission range T . We leave other protocols like, two hop routing, direct contact, spray and wait [START_REF] Spyropoulos | Spray and wait: an efficient routing scheme for intermittently connected mobile networks[END_REF] for future study. Also, we assume that the transmissions are error free and instantaneous. Further, all the nodes are identical and participate equally in relaying the message. In the simulator it is possible to set a certain fraction of MNs not to participate in message spreading. However, in this paper we restrict to the case where all the nodes participate in message spreading.

In epidemic routing the copy of the message is available with a large number of MNs and it is likely that it will reach the destination with high probability. But the time required for the message to reach the destination, or the delay in delivery of the message, can be still very high. For example, small MN density, small transmission range, can result in longer time for two MNs to be within transmission distance of each other hence leading to large time delays. In DTNs, the message delivery is of primary interest than the time taken it to be delivered at the destination, but the sooner the message gets delivered the better: in some case the message may lose its relevance if it fails to reach within a stipulated time. Thus, it is of interest to know probability of message taking too long to be delivered at the destination. Or, more generally, to know the tail distribution of the delay.

IV. EXPERIMENTAL SETUP AND STATISTICAL ANALYSIS

A. Experimental Setup

In this section we discuss simulation setup and data collection for empirical study of the delay distribution. In the simulator we set the parameters as follows: R 1 = R 3 = 100, R 2 = 25, γ 1 = 0.001, γ 2 = 1.0, γ 3 = 0, v 0 = 2, ∆ = 0.095, T = 25, v max = 20. The choice of these parameters is based on visual monitoring, such that the simulated mobility looks close to that of a crowd mobility. The above parameters are set as default values in the simulator. The simulation is run with 50 MNs that start at random locations with M N 1 as the source and M N n as the destination. In each run, we set a warmup period of 30 seconds before collecting any data. Once this warmup period lapses, we start a timer at the destination and the note the time when it receives the message for the first time. We made 500 runs to collect samples of the random variable D.

In this paper, we present statistical analysis with the parameter set to default values. However, we note that the same observations continues to hold for other set of parameters as well. In [18] we present numerical results for a range of parameters and some preliminary observation on sensitivity analysis. In the following subsections, we explain the statistical analysis for one set of samples collected from the simulator. These samples are available online at [17].

B. Statistical Analysis

In this section we analyze the delay distribution using statistical tools available in MATLAB. We obtain the empirical distribution using the samples generated through the simulator. It is observed that the exponent of the tail distribution is linear, and the slope is obtained using linear regression. Finally, we derive the theoretical distribution that best represents the tail of the delay distribution.

Let F D (•) denote the cumulative distribution function (CDF) of D. We will characterize the delay distribution by evaluating the empirical distribution (ECDF) of D which we denote as FD (•). We also calculated upper and lower bounds (see Figure 3) that corresponds to 95% confidence intervals of the ECDF [START_REF] Greenwood | The natural duration of cancer[END_REF]. We next evaluate the exponent of the tail distribution, i.e., rate of change or slope of:

y(d) = -log(1 -FD (d)).
This function is plotted in Figure 4, as the ECDF, with its upper and lower confidence interval bounds. As seen from Figure 4, y(d) is an almost linear function for large values of D. In this subsection we evaluate the slope of y(d) using linear regression.

Let d = (d 1 , d 2 , • • • , d m ) and f = (f 1 , f 2 , • • • , f m )
denote the time sequence and the corresponding values of the cumulative probability. Define

y i = -log(1- f i ) for each i = 1, 2, • • • m and denote y = (y 1 , y 2 , • • • , y m ).
As our interest is in the behavior of the tail distribution, we take into account only the last 40% of sample of the function (d, y) to evaluate the slope. Let M denote the size of this truncated vector (d, y), and re-index its samples from 1 to M. For the samples we have, the sample size after truncation is M = 201.

Before we proceed to apply regression analysis on the samples, let us first verify that there is significant linear relationship between d and y using regression the t-test. The correlation coefficient for the pair (d, y) is R = 0.9960, and the corresponding test statistics for the data is:

t = r M -2 1 -R 2 = 156.7324,
whereas the critical value is 1.972 at significance level 0.05.

Clearly, at significance level 0.05, there is a significant linear relation between d and y as t > 1.9720.

We next evaluate the slope and intercept of the linear function that best represents the linear relation between d and y. These values are given, respectively, as follows:

m = M d i y i -( d i )( y i ) M d 2 i -( d i ) 2 = 0.0001443 and: b = y i -m d i M = -2.602.
The best fit linear function in shown in figure 4 in a continuous line. Note that, in the region of interest, its tail, the linear function lies within the upper and lower confidence interval bounds. Thus the logarithm of the tail of the ECDF can be approximated by the linear line with slope and intercept as computed above, and the confidence interval for this fit

Data quantiles

Fig. 5. q-q plot is 95 %. Also, the adjusted R 2 value as a measure of the coefficient of determination for this fit is 0.9919. Hence, with a high confidence value, we can approximate y(d) by a linear function, which in turn implies that tail of the empirical distribution can be approximated by an exponential function, i.e.:

Pr{D ≥ d} ≈ exp{-md} f or d > d * ,
where d * = d N -M is the threshold beyond which we apply the regression to fit y(d) to the data, with m as the slope of the best fit linear function. Note that in arriving at the exponential approximation for the tail distribution we only took into account the slope of the linear regression and not the intercept. In the next section we study the goodness of fit of the exponential distribution to the tail of the empirical distribution by studying quantile-quantile plots and performing a Kolmogorov-Smirnov test.

D. Quantile-Quantile Plot

The exponential approximation of the empirical distribution is good when the samples are larger than the threshold d * . Let D denote the distribution of the delay conditioned that it is larger than d * . Then the CDF of D is given by

F D (d) = 1 -exp{-m(d -d * )}, if d ≥ d * 0, otherwise (1)
The q-quantile of this CDF for a given p ∈ [0 1] can be evaluated as

F -1 D (p) = d * - log(1 -p) m ,
where F -1 D (p) is the inverse of the CDF of the tail. We compare this theoretical quantile function with the quantile of the ECDF obtained from the samples collected. For a fair comparison, we discard all the samples that are smaller that d * in computing the ECDF. The q-q plot is shown in figure 5. Note that the quantile functions closely follows the 45 degree line, with deviations only at the edges. This is an expected behavior as we are approximating an exponential function with a linear fit. Thus, we overstimate the quantiles at the lower end of the curve, while overstimating them at the upper end. Nevertheless, the correlation coefficient between the two quantile functions is 0.9938. Hence, the exponential approximation of the tail distribution is a good fit.

We next perform the Kolomogrov-Smirnov (KS) test on the data samples. The test statistic for this test is

D = max d |F D(d) -F D (d)| = 0.0978,
where F D denotes the ECDF of the samples that takes value larger than d * in the data. At a significance level of 5%, the hypothesis that the samples are drawn from the exponential distribution with λ = m is not rejected by the KS test. Thus, it can be believed that the tail is exponentially distributed.

V. DELAY-ENERGY TRADEOFF

In this section we study the effect of transmission range T on the mean message delivery time. If nodes can transmit at higher power, then they can transmit the message successfully over a larger distance. However, since the nodes are battery powered, this will require the nodes to recharge often. In Figure 6 we plot the mean message delivery time versus the transmission range T . In generating the plots we set the parameters to the same values given in Section III. Note that if T < 25, then the nodes come rarely within the transmission range of each other as they will try to maintain a separation of 25 units (R 2 = 25). Thus, we generate the plot for T ≥ 25. We also obtain the linear regression of the data points and plot it in Figure 6. From the figure we note that there is a linear relation between the mean message delivery time and the transmission range. Thus, to decrease the mean transmission delay by a certain factor, the MS need to increase the transmission range by the same factor. Mobility modeling is crucial to analyze the performance of wireless networks. To understand the performance of a protocol, simplifying assumptions, like independence of node mobility, are made in wireless networks, which is not realistic. In this paper, we proposed a model to study crowd mobility that captures dependency among the movements of individuals in a group. We implemented a simulator in which the degree of dependency of movements among the individuals can be varied by controlling the parameters. We established through systematic statistical analysis that the distribution of message delivery time is light tailed with exponential distribution.

In this preliminary study, we have not take into account errors in transmissions and minimum contact duration for successful transmission. In future, we like to bring in these aspects into simulations. Also, in the current model, when two nodes come close to each other they sometimes make a complete U-turn, which is not realistic. We would like to make such turns smooth, and study how this impacts the results.

Fig. 1 .

 1 Fig. 1. Snap shot of Crowd Mobility simulator: Blue node has message

  Fig. 3. empirical CDF of delay

Fig. 4 .

 4 Fig. 4. Exponent of tail distribution

Fig. 6 .

 6 Fig. 6. Transmission range vs mean delay