**Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach**

Ciro D'Apice1, Rosanna Manzo1 and Benedetto Piccoli2

<sup>1</sup>*Department of Electronic and Information Engineering, University of Salerno, Fisciano (SA)* <sup>2</sup>*Department of Mathematical Sciences, Rutgers University, Camden, New Jersey* <sup>1</sup>*Italy* <sup>2</sup>*USA*

#### **1. Introduction**

420 Telecommunications Networks – Current Status and Future Trends

Bandler, J. W., Srinivasan, T.V., Charalambous, (1972) Minimax Optimization of networks by Grazor Search, IEEE Trans. Microwave Theory Tech., MTT-20, 596-604. Bandler, J. W., Charalambous, C. (1972), Practical least pth optimization of networks, IEEE

Brayton, R.K., S.W. Director, G.D. Hachtel, and L.Vidigal, (1979), A New Algorithm for

Banner, R., Orda A., (2007), Multipath Routing Algorithms for Congestion Minimization, IEEE/ACM Transactions on Networking, Vol. 15, No. 2, pp. 413-424 Chen, J., Chan, S.-H. Gary, Li V. O. K., (2004) Multipath Routing for Video Delivery Over

Lin, X., Shroff, N. B., (2006)"Utility Maximization for Communication Networks With

Güven, T., La, R. J., Shayman, M. A., Bhattacharjee, B., (2008), A Unified Framework for

Jaffe J. M., "Bottlneck Flow Control", IEEE Transactions on Communications, Vol 29, No 7,

Tsai, D., Liau, T. C., Tsai, Wei K., (2006) Least Square Approach to Multipath Maxmin Rate

Tsai, Wei K., Kim, Y., (1999) Re-Eximining Maxmin Protocols: A Fundamental Study on

Delli Priscoli, F., (2010) A Fully Cognitive Approach for Future Internet, Future Internet ISSN 1999-5903 available at www.mdpi.com/journal/futureinternet.

Malozemov, V. N., (1974) Introduction to minimax, John Wiley & Sons. Cidon, I., Rom R., Shavitt Y., (1999) Analysis of Multipath Routing, IEEE/ACM Transactions

Statistical Circuit Design Based on Quasi-Newton Methods and Function Splitting, IEEE Trans. Circuits and Systems, Vol. CAS-26, pp. 784-794. Demyanov, V. F.,

Bandwidth-Limited Networks, IEEE Journal on Selected Areas in Communications,

Multipath Routing", IEEE Transactions on Automatic Control, Vol. 51, No. 5, pp.

Multipath Routing for Unicast and Multicast Traffic, IEEE/ACM Transactions on

Allocation, 14th IEEE International Conference on Networks, 2006 (ICON '06), Vol.

Convergence, Complexity, Variations and Performance, 18th Annual Joint Conference of the IEEE Computer and Communications Societies (INFOCOM '99),

Trans. Microwave Theory Tech., MTT-20, 834-840.

ON Networking, Vol. 7, No. 6, pp. 885-896

Networking, Vol. 16, No. 5, pp. 1038-1051

Vol. 22, No. 10, pp. 1920-1932

July 1981, pp. 954-962

Vol. 2, March 1999, pp. 811-818

766-781

1, pp. 1-6.

There are various approaches to telecommunication and data networks (see for example Alderson et al. (2007), Baccelli et al. (2006), Baccelli et al. (2001), Kelly et al. (1998), Tanenbaum (1999), Willinger et al. (1998)). A first model for data networks, similar to that used for car traffic, has been proposed in D'Apice et al. (2006), where two algorithms for dynamics at nodes were considered and existence of solutions to Cauchy Problems was proved. Then in D'Apice et al. (2008), following the approach of Garavello et al. (2005) for road networks (see also Coclite et al. (2005); Daganzo (1997); Garavello et al. (2006); Holden et al. (1995); Lighthill et al. (1955); Newell (1980); Richards (1956)), sources and destinations have been introduced, thus taking care of the packets paths inside the network.

In this Chapter we deal with the fluid-dynamic model for data networks together with optimization problems, reporting some results obtained in Cascone et al. (2010); D'Apice et al. (2006; 2008; 2010).

A telecommunication network consists in a finite collection of transmission lines, modelled by closed intervals of **R** connected by nodes (routers, hubs, switches, etc.). Taking the Internet network as model, we assume that:


Since each lost packet is sent again until it reaches next node, looking at macroscopic level, it is assumed that the packets number is conserved. This leads to a conservation law for the packets density *ρ* on each line:

$$
\rho\_t + f\left(\rho\right)\_x = 0.\tag{1}
$$

The flux *f*(*ρ*) is given by *v*(*ρ*)· *ρ* where *v* is the average speed of packets among nodes, derived considering the amount of packets that may be lost.

Due to the nonlinear relation among cost functionals, the optimization of velocity and travel

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 423

The analytical treatment of a complex network is very hard due to the high nonlinearity of the dynamics and discontinuities of the I/O maps. For these reasons, a decentralized strategy has

Step 1. The optimal controls for asymptotic costs in the case of a single node with constant

Step 2. For a complex network, the (locally) optimal parameters at every node are used. Thus,

The optimization problem for nodes of 2 × 2 type, i.e. with two entering and two exiting lines, and traffic distribution coefficient *α* and priority parameter *p* as control parameters, constant

Then a test telecommunication network, consisting of 24 nodes, each one of 2 × 2 type has been studied. Three different choices have been tested for the traffic distribution coefficients and priority parameters: (locally) optimal, static random and dynamic random. The first choice is given by Step 1. By static random parameters, we mean a random choice done at the beginning of the simulation and then kept constant. Finally, dynamic random coefficients are

The results present some interesting features: the performances of the optimal coefficients are definitely superior with respect to the other two. Then, how the dynamic random choice, which sometimes is equal in performance to the optimal ones, may be not feasible

The Chapter is organized as follows. Section 2 reports the model for data networks. Then, in Section 3, we consider possible choices of the traffic distribution functions, and how to compute the traffic distribution matrix from the latter functions and the traffic-type function. We describe two routing algorithms, giving explicit unique solutions to RPs. In Section 4, we discuss the validity of the assumption on the loss probability function, the velocity and flux. The subsequent Section 5 is devoted to the analysis of the optimal control problem introducing the cost functionals. Simulations for three different choices of parameters (optimal, static and dynamic random) in the case of a complex network are presented. The paper ends with

A telecommunication network is a finite collection of transmission lines connected together by nodes, some of which are sources and destinations. Formally we introduce the following

**Definition 1.** *A telecommunication network is given by a* 7*-tuple* (*N*, I*,* F*,* J *,* S*,* D*,* R) *where*

**Nodes** J *is a collection of subsets of* {±1, ..., ±*N*} *representing nodes. If j* ∈ *J* ∈ J *, then the transmission line I*|*j*<sup>|</sup> *is crossing at J as incoming line (i.e. at point bi) if j* <sup>&</sup>gt; <sup>0</sup> *and as outgoing line*

*<sup>i</sup>* ] �→ **R***, i* = 1, ..., *N;*

**Cardinality** *N is the cardinality of the network, i.e. the number of lines in the network;* **Lines** I *is the collection of lines, modelled by intervals Ii* = [*ai*, *bi*] ⊆ **R***, i* = 1, ..., *N;*

the optimal control is determined at each node independently.

initial data and asymptotic functionals has been completely solved.

chosen randomly at every instant of time for every node.

for modelling and robustness reasons has been discussed.

**Fluxes** <sup>F</sup> *is the collection of flux functions fi* : [0, *<sup>ρ</sup>*max

time can give different control parameters.

been adapted as follows:

conclusions in Section 6.

**2. Basic definitions**

definition:

initial data is computed.

The key point of the model is the loss probability, used to define the flux function. Indeed the choice of a non reasonable loss probability function could invalidate the model. To achieve the goal of the validation of the model assumptions, the loss probability function has been compared with the behaviour of the packet loss derived from known models used in literature to infer network performance and the shape of the velocity and flux functions has been discussed. All the comparisons confirm the validity of the assumptions underlying the fluid-dynamic model (see D'Apice et al. (2010)).

To describe the evolution of networks in which many lines intersect, Riemann Problems (RPs) at junctions were solved in D'Apice et al. (2006) proposing two different routing algorithms:


One of the drawback of (RA2) is that it does not take into account the global path of packets, therefore leading to possible cycling to bypass congested nodes. These cyclings are avoided if we consider that the packets originated from a source and with an assigned destination have paths inside the network.

Taking this in mind the model was refined in D'Apice et al. (2008). On each transmission line a vector *π* describing the traffic types, i.e. the percentages of packets going from a source to a destination, has been introduced. Assuming that packets velocity is independent from the source and the destination, the evolution of *π* follows a semilinear equation

$$
\pi\_t + v(\rho)\pi\_x = 0,\tag{2}
$$

hence inside transmission lines the evolution of *π* is influenced by the average speed of packets.

Different distribution traffic functions describing different routing strategies have been analysed:


In particular two ways according to which the traffic at a junction is splitted towards the outgoing lines have been defined. Starting from the distribution traffic function, and using the vector *π*, the traffic distribution matrix, which describes the percentage of packets from an incoming line that are addressed to an outgoing one, has been assigned. Then, methods to solve RPs according to the routing algorithms (RA1) and (RA2) have been proposed. Optimizations results have been obtained for the model consisting of the conservation law (1). In particular priority parameters and traffic distribution coefficients have been considered as controls and two functionals to measure the efficiency of the network have been defined in Cascone et al. (2010):


Due to the nonlinear relation among cost functionals, the optimization of velocity and travel time can give different control parameters.

The analytical treatment of a complex network is very hard due to the high nonlinearity of the dynamics and discontinuities of the I/O maps. For these reasons, a decentralized strategy has been adapted as follows:


The optimization problem for nodes of 2 × 2 type, i.e. with two entering and two exiting lines, and traffic distribution coefficient *α* and priority parameter *p* as control parameters, constant initial data and asymptotic functionals has been completely solved.

Then a test telecommunication network, consisting of 24 nodes, each one of 2 × 2 type has been studied. Three different choices have been tested for the traffic distribution coefficients and priority parameters: (locally) optimal, static random and dynamic random. The first choice is given by Step 1. By static random parameters, we mean a random choice done at the beginning of the simulation and then kept constant. Finally, dynamic random coefficients are chosen randomly at every instant of time for every node.

The results present some interesting features: the performances of the optimal coefficients are definitely superior with respect to the other two. Then, how the dynamic random choice, which sometimes is equal in performance to the optimal ones, may be not feasible for modelling and robustness reasons has been discussed.

The Chapter is organized as follows. Section 2 reports the model for data networks. Then, in Section 3, we consider possible choices of the traffic distribution functions, and how to compute the traffic distribution matrix from the latter functions and the traffic-type function. We describe two routing algorithms, giving explicit unique solutions to RPs. In Section 4, we discuss the validity of the assumption on the loss probability function, the velocity and flux. The subsequent Section 5 is devoted to the analysis of the optimal control problem introducing the cost functionals. Simulations for three different choices of parameters (optimal, static and dynamic random) in the case of a complex network are presented. The paper ends with conclusions in Section 6.

#### **2. Basic definitions**

2 Telecommunications Networks

The key point of the model is the loss probability, used to define the flux function. Indeed the choice of a non reasonable loss probability function could invalidate the model. To achieve the goal of the validation of the model assumptions, the loss probability function has been compared with the behaviour of the packet loss derived from known models used in literature to infer network performance and the shape of the velocity and flux functions has been discussed. All the comparisons confirm the validity of the assumptions underlying the

To describe the evolution of networks in which many lines intersect, Riemann Problems (RPs) at junctions were solved in D'Apice et al. (2006) proposing two different routing algorithms: (RA1) Packets from incoming lines are sent to outgoing ones according to their final destination (without taking into account possible high loads of outgoing lines); (RA2) Packets are sent to outgoing lines in order to maximize the flux through the node.

One of the drawback of (RA2) is that it does not take into account the global path of packets, therefore leading to possible cycling to bypass congested nodes. These cyclings are avoided if we consider that the packets originated from a source and with an assigned destination have

Taking this in mind the model was refined in D'Apice et al. (2008). On each transmission line a vector *π* describing the traffic types, i.e. the percentages of packets going from a source to a destination, has been introduced. Assuming that packets velocity is independent from the

hence inside transmission lines the evolution of *π* is influenced by the average speed of

Different distribution traffic functions describing different routing strategies have been

• at a junction the traffic started at source *s* and with *d* as final destination, coming from the

• at a junction the traffic started at source *s* and with *d* as final destination, coming from the

In particular two ways according to which the traffic at a junction is splitted towards the outgoing lines have been defined. Starting from the distribution traffic function, and using the vector *π*, the traffic distribution matrix, which describes the percentage of packets from an incoming line that are addressed to an outgoing one, has been assigned. Then, methods to solve RPs according to the routing algorithms (RA1) and (RA2) have been proposed. Optimizations results have been obtained for the model consisting of the conservation law (1). In particular priority parameters and traffic distribution coefficients have been considered as controls and two functionals to measure the efficiency of the network have been defined in

transmission line *i*, is routed on every outgoing lines or on some of them.

*π<sup>t</sup>* + *v*(*ρ*)*π<sup>x</sup>* = 0, (2)

source and the destination, the evolution of *π* follows a semilinear equation

transmission line *i*, is routed on an assigned line *j*;

1) The velocity of packets travelling through the network. 2) The travel time taken by packets from source to destination.

fluid-dynamic model (see D'Apice et al. (2010)).

paths inside the network.

packets.

analysed:

Cascone et al. (2010):

A telecommunication network is a finite collection of transmission lines connected together by nodes, some of which are sources and destinations. Formally we introduce the following definition:

**Definition 1.** *A telecommunication network is given by a* 7*-tuple* (*N*, I*,* F*,* J *,* S*,* D*,* R) *where*

**Cardinality** *N is the cardinality of the network, i.e. the number of lines in the network;*

**Lines** I *is the collection of lines, modelled by intervals Ii* = [*ai*, *bi*] ⊆ **R***, i* = 1, ..., *N;*

**Fluxes** <sup>F</sup> *is the collection of flux functions fi* : [0, *<sup>ρ</sup>*max *<sup>i</sup>* ] �→ **R***, i* = 1, ..., *N;*

**Nodes** J *is a collection of subsets of* {±1, ..., ±*N*} *representing nodes. If j* ∈ *J* ∈ J *, then the transmission line I*|*j*<sup>|</sup> *is crossing at J as incoming line (i.e. at point bi) if j* <sup>&</sup>gt; <sup>0</sup> *and as outgoing line*

given by the averaged density times the average velocity. A possible choice of *p* is the

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 425

0, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> *<sup>σ</sup>*, *<sup>ρ</sup>max* (*ρ*−*σ*)

*<sup>v</sup>*¯, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> *<sup>σ</sup>*,

*<sup>v</sup>*¯*ρ*, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> *<sup>σ</sup>*, *<sup>v</sup>*¯*σ*(*ρmax*−*ρ*)

Σ Ρ*max*

To simplify the treatment of the corresponding conservation laws, we will assume the

(*F*) Setting *ρmax* = 1, on each line the flux *f* : [0, 1] → *R* is concave, *f*(0) = *f*(1) = 0 and there

<sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*, <sup>1</sup>

where we use the assumption (*F*). Therefore, the network load evolution is described by a

*<sup>i</sup>* ]. Moreover, inside each line *Ii* we define a traffic-type function *πi*, which measures the portion of the whole density coming from each source and travelling towards each destination:

*<sup>ρ</sup>*, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> <sup>1</sup>

Ρ

2 ,

*f* (*ρ*) = *ρ* (1 − *ρ*), *ρ* ∈ [0, 1] , (8)

*∂tρ<sup>i</sup>* + *∂<sup>x</sup> fi* (*ρi*) = 0, (9)

*v <sup>2</sup>*�Σ � Ρ Ρ

*v*¯ *<sup>σ</sup>*(*ρmax*−*ρ*)

*<sup>ρ</sup>* (*ρmax*−*σ*), *<sup>σ</sup>* <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> *<sup>ρ</sup>*max, (4)

*<sup>ρ</sup>* (*ρmax*−*σ*), *<sup>σ</sup>* <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> *<sup>ρ</sup>*max, (5)

*<sup>ρ</sup>max*−*<sup>σ</sup>* , *<sup>σ</sup>* <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> *<sup>ρ</sup>*max. (6)

*<sup>f</sup>* <sup>Ρ</sup>

Σ Ρ*max*

Ρ

<sup>2</sup> , *v*¯ = 1.

*<sup>v</sup>* <sup>Ρ</sup> *<sup>v</sup><sup>2</sup>*Σ�Ρ

<sup>2</sup> <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> 1, (7)

*p* (*ρ*) =

*v* (*ρ*) =

*f* (*ρ*) =

Ρ

*v*

Fig. 1. Loss probability, average velocity and flux behaviours for *ρmax* = 1, *σ* = <sup>1</sup>

*f* (*ρ*) =

On each transmission line *Ii* we consider the evolution equation

*<sup>v</sup>* <sup>Ρ</sup>

Σ

Notice that the "tent" function

and the parabolic flux

satisfy the assumption (F).

**2.2 Dynamics on the network**

exists a unique maximum point *σ* ∈]0, 1[.

finite set of functions *<sup>ρ</sup><sup>i</sup>* : [0, <sup>+</sup>∞[ <sup>×</sup> *Ii* �→ [0, *<sup>ρ</sup>*max

*<sup>2</sup>* Ρ�Σ Ρ

*<sup>p</sup>* <sup>Ρ</sup>

following:

from which

following:

and

*(i.e. at point ai) if j* < 0*. For each junction J* ∈ J *, we indicate by* Inc(*J*) *the set of incoming lines, that are Ii's such that i* ∈ *J, while by* Out(*J*) *the set of outgoing lines, that are Ii's such that* −*i* ∈ *J. We assume that each line is incoming for (at most) one node and outgoing for (at most) one node;*


#### **2.1 Dynamics on lines**

Following D'Apice et al. (2008), we recall the model used to define the dynamics of packet densities along lines. We make the following hypothesis:


Since the packet transmission velocity on the line is assumed constant, it is possible to compute an average velocity function and thus an average flux function.

Let us focus on two consecutive nodes *Nk* and *Nk*+1, assume a static situation, i.e. *Rk* and *Rk*<sup>+</sup><sup>1</sup> are constant. Indicate by *δ* the distance between the nodes, Δ*tav* the packets average transmission time, *v*¯ = *<sup>δ</sup>* <sup>Δ</sup>*t*<sup>0</sup> the packet velocity without losses and *<sup>v</sup>* <sup>=</sup> *<sup>δ</sup>* <sup>Δ</sup>*tav* the average packets velocity. Then, we can compute:

$$\Delta t\_{a\upsilon} = \sum\_{n=1}^{M} n \Delta t\_0 (1 - p(R\_{k+1})) p^{n-1} (R\_{k+1})\_{\upsilon}$$

where *M* = [*T*/Δ*t*0] (here [·] indicates the floor function) represents the number of attempts of sending a packet and *T* is the length of a processing time slot. The hypothesis (H4) corresponds to assume Δ*t*<sup>0</sup> << *T* or, equivalently, *M* ∼ +∞. Making the identification, *M* = +∞, we get:

$$
\Delta t\_{av} = \frac{\Delta t\_0}{1 - p(R\_{k+1})}.
$$

$$
v = \overline{v}(1 - p(R\_{k+1})).\tag{3}
$$

and

Let us call now *ρ* the averaged density and *ρmax* its maximum. We can interpret the probability loss function *p* as a function of *ρ* and, using (3), determine the corresponding flux function, given by the averaged density times the average velocity. A possible choice of *p* is the following:

$$p\left(\rho\right) = \begin{cases} 0, & 0 \le \rho \le \sigma, \\ \frac{\rho\_{\max}\left(\rho - \sigma\right)}{\rho\left(\rho\_{\max} - \sigma\right)}, \sigma \le \rho \le \rho\_{\max} \end{cases} \tag{4}$$

from which

$$v\left(\rho\right) = \begin{cases} \overline{v}, & 0 \le \rho \le \sigma, \\ \overline{v} \frac{\sigma(\rho\_{\max} - \rho)}{\rho \left(\rho\_{\max} - \sigma\right)}, \sigma \le \rho \le \rho\_{\max} \end{cases} \tag{5}$$

and

4 Telecommunications Networks

**Destinations** D *is the subset of* {1, ..., *N*} *representing lines leading to traffic destinations, Thus,*

**Traffic distribution functions** R *is a finite collection of functions (also multivalued) rJ* : Inc(*J*) × S×D→ Out(*J*)*. For every J, rJ*(*i*,*s*, *d*) *indicates the outgoing direction of traffic that started at*

Following D'Apice et al. (2008), we recall the model used to define the dynamics of packet

(H1) Lines are composed of consecutive processors *Nk*, which receive and send packets. The

(H2) There are two time-scales: Δ*t*0, the physical travel time of a single packet from node to node (assumed to be independent of the node for simplicity); *T*, the processing time,

(H3) Each processor *Nk* tries to send all packets *Rk* at the same time. Packets are lost according to a loss probability function *p* : [0, *Rmax*] → [0, 1], computed at *Rk*+1, and lost packets are

(H4) The number of packets not transmitted for a whole processing time slot is negligible. Since the packet transmission velocity on the line is assumed constant, it is possible to

Let us focus on two consecutive nodes *Nk* and *Nk*+1, assume a static situation, i.e. *Rk* and *Rk*<sup>+</sup><sup>1</sup> are constant. Indicate by *δ* the distance between the nodes, Δ*tav* the packets average

<sup>Δ</sup>*t*<sup>0</sup> the packet velocity without losses and *<sup>v</sup>* <sup>=</sup> *<sup>δ</sup>*

where *M* = [*T*/Δ*t*0] (here [·] indicates the floor function) represents the number of attempts of sending a packet and *T* is the length of a processing time slot. The hypothesis (H4) corresponds to assume Δ*t*<sup>0</sup> << *T* or, equivalently, *M* ∼ +∞. Making the identification,

Let us call now *ρ* the averaged density and *ρmax* its maximum. We can interpret the probability loss function *p* as a function of *ρ* and, using (3), determine the corresponding flux function,

1 − *p*(*Rk*<sup>+</sup>1)

,

*v* = *v*¯(1 − *p*(*Rk*<sup>+</sup>1)). (3)

<sup>Δ</sup>*tav* <sup>=</sup> <sup>Δ</sup>*t*<sup>0</sup>

*<sup>n</sup>*Δ*t*0(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*(*Rk*<sup>+</sup>1))*pn*−1(*Rk*<sup>+</sup>1),

<sup>Δ</sup>*tav* the average packets

during which each processor tries to operate the transmission of a given packet;

*and only if j is not outgoing for any node. We assume that* S �= ∅*;*

densities along lines. We make the following hypothesis:

packets number at *Nk* is indicated by *Rk* ∈ [0, *Rmax*];

Δ*tav* =

sent again for a time slot of length *T*;

**2.1 Dynamics on lines**

transmission time, *v*¯ = *<sup>δ</sup>*

*M* = +∞, we get:

and

velocity. Then, we can compute:

*j* ∈ D *if and only if j is not incoming for any node. We assume that* D �= ∅*;*

*source s has d as final destination and reached J from the incoming road i.*

compute an average velocity function and thus an average flux function.

*M* ∑ *n*=1

*(i.e. at point ai) if j* < 0*. For each junction J* ∈ J *, we indicate by* Inc(*J*) *the set of incoming lines, that are Ii's such that i* ∈ *J, while by* Out(*J*) *the set of outgoing lines, that are Ii's such that* −*i* ∈ *J. We assume that each line is incoming for (at most) one node and outgoing for (at most) one node;* **Sources** S *is the subset of* {1, ..., *N*} *representing lines starting from traffic sources. Thus, j* ∈ S *if*

$$f\left(\rho\right) = \begin{cases} \bar{v}\rho, & 0 \le \rho \le \sigma, \\ \frac{\sigma\tau(\rho\_{\max}-\rho)}{\rho\_{\max}-\sigma}, \sigma \le \rho \le \rho\_{\max}. \end{cases} \tag{6}$$

Fig. 1. Loss probability, average velocity and flux behaviours for *ρmax* = 1, *σ* = <sup>1</sup> <sup>2</sup> , *v*¯ = 1.

To simplify the treatment of the corresponding conservation laws, we will assume the following:

(*F*) Setting *ρmax* = 1, on each line the flux *f* : [0, 1] → *R* is concave, *f*(0) = *f*(1) = 0 and there exists a unique maximum point *σ* ∈]0, 1[.

Notice that the "tent" function

$$f\left(\rho\right) = \begin{cases} \rho\_{\prime} & 0 \le \rho \le \frac{1}{2} \\ 1 - \rho\_{\prime} \frac{1}{2} \le \rho \le 1 \end{cases} \tag{7}$$

and the parabolic flux

$$f\left(\rho\right) = \rho\left(1-\rho\right), \rho \in \left[0,1\right],\tag{8}$$

satisfy the assumption (F).

#### **2.2 Dynamics on the network**

On each transmission line *Ii* we consider the evolution equation

$$
\partial\_l \rho\_{\mathbf{i}} + \partial\_{\mathbf{x}} f\_{\mathbf{i}} \left( \rho\_{\mathbf{i}} \right) = 0,\tag{9}
$$

where we use the assumption (*F*). Therefore, the network load evolution is described by a finite set of functions *<sup>ρ</sup><sup>i</sup>* : [0, <sup>+</sup>∞[ <sup>×</sup> *Ii* �→ [0, *<sup>ρ</sup>*max *<sup>i</sup>* ].

Moreover, inside each line *Ii* we define a traffic-type function *πi*, which measures the portion of the whole density coming from each source and travelling towards each destination:

**3. Riemann Solvers at junctions**

line as follows:

Finally denote with

solution of the RP at the junction.

1) *rJ* : Inc(*J*) ×S×D→ Out(*J*);

*rJ*(*i*,*s*, *d*) in two different ways:

**2b)** *rJ* : Inc(*J*) ×S×D→ [0, 1]

*i*,*s*,*d*,*j*

*rj*(*i*,*s*, *d*) ⊆ Out(*J*);

*rJ*(*i*,*s*, *<sup>d</sup>*)=(*αi*,*s*,*d*,*n*+<sup>1</sup>

with 0 ≤ *α*

outgoing lines.

*α i*,*s*,*d*,*j <sup>J</sup>* .

**2a)** *rJ* : Inc(*J*) ×S×D → Out(*J*),

and

the outgoing ones and by (*ρ*1,0, ..., *ρn*+*m*,0) the initial datum.

*γ*max *<sup>i</sup>* =

*γ*max *<sup>j</sup>* =

Consider a junction *J* of *n* × *m* type. We denote with *ρi*(*t*, *x*), *i* = 1, ..., *n* and *ρj*(*t*, *x*), *j* = *n* + 1, ..., *n* + *m* the traffic densities, respectively, on the incoming transmission lines and on

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 427

Define the maximum flux that can be obtained by a single wave solution on each transmission

*<sup>i</sup>* ], *i* = 1, ..., *n*,

and with *γ*ˆ*inc* = (*f*(*ρ*ˆ*i*), ..., *f*(*ρ*ˆ*n*)), *γ*ˆ *out* = (*f*(*ρ*ˆ*n*+1), ..., *f*(*ρ*ˆ*n*+*m*)) where *ρ*ˆ = (*ρ*ˆ1, ..., *ρ*ˆ*n*+*m*) is the

If *rJ* is of type 1), then each packet has a deterministic route, it means that, at the junction *J*, the traffic that started at source *s* and has *d* as final destination, coming from the transmission

Instead if *rJ* is of type 2), at the junction *J*, the traffic with source *s* and destination *d* coming from a line *Ii* is routed on every line *Ij* ∈ Out(*J*) or on some lines *Ij* ∈ Out(*J*). We can define

> *n*+*m* ∑ *j*=*n*+1 *α i*,*s*,*d*,*j <sup>J</sup>* = 1.

In case 2a) we have to specify in which way the traffic at junction *J* is splitted towards the

The definition 2b) means that, at the junction *J*, the traffic with source *s* and destination *d* coming from line *Ii* is routed on the outgoing line *Ij*, *j* = *n* + 1, ..., *n* + *m* with probability

*<sup>j</sup>* ], *j* = *n* + 1, ..., *n* + *m*,

*<sup>f</sup>*(*σ*), if *<sup>ρ</sup>i*,0 <sup>∈</sup> ]*σ*, 1] , *<sup>i</sup>* <sup>=</sup> 1, ..., *<sup>n</sup>*, (11)

*<sup>f</sup>*(*ρj*,0), if *<sup>ρ</sup>j*,0 <sup>∈</sup> ]*σ*, 1] , *<sup>j</sup>* <sup>=</sup> *<sup>n</sup>* <sup>+</sup> 1, ..., *<sup>n</sup>* <sup>+</sup> *<sup>m</sup>*. (12)

*<sup>f</sup>*(*ρi*,0), if *<sup>ρ</sup>i*,0 <sup>∈</sup> [0, *<sup>σ</sup>*],

*<sup>f</sup>*(*σ*), if *<sup>ρ</sup>j*,0 <sup>∈</sup> [0, *<sup>σ</sup>*],

Ω*<sup>i</sup>* = [0, *γ*max

Ω*<sup>j</sup>* = [0, *γ*max

Now, we discuss some possible choices for the traffic distribution function:

Out(*J*) ,

Let us analyze how the distribution matrix *A* is constructed using *π* and *rJ*.

*<sup>J</sup>* )

*<sup>J</sup>* ≤ 1, *j* ∈ {*n* + 1, ..., *n* + *m*},

2) *rJ* : Inc(*J*) ×S×D → Out(*J*), i.e. *rJ* is a multifunction.

line *Ii*, is routed on an assigned line *Ij* (*rJ*(*i*,*s*, *d*) = *j*).

*<sup>J</sup>* , ..., *<sup>α</sup>i*,*s*,*d*,*n*+*<sup>m</sup>*

**Definition 2.** *A traffic-type function on a line Ii is a function*

$$
\pi\_i \colon [0, \infty[ \times [a\_i, b\_i] \times \mathcal{S} \times \mathcal{D} \mapsto [0, 1] ]
$$

*such that, for every t* ∈ [0, ∞[ *and x* ∈ [*ai*, *bi*]

$$\sum\_{s \in \mathcal{S}, d \in \mathcal{D}} \pi\_i(t, \mathbf{x}, \mathbf{s}, d) = 1.$$

In other words, *πi*(*t*, *x*,*s*, *d*) specifies the density fraction *ρi*(*t*, *x*) that started from source *s* and is moving towards the final destination *d*.

Assuming, on the discrete model, that a FIFO policy is used at nodes, it is natural that the averaged velocity, obtained in the limit procedure, is independent from the original sources of packets and their final destinations. In other words, we make the following hypothesis:

(H5) On each line *Ii*, the average velocity of packets depends only on the value of the density *ρ<sup>i</sup>* and not on the values of the traffic-type function *πi*.

As a consequence of hypothesis (H5), we deduce the semilinear equation

$$
\partial\_t \pi\_i(t, \mathbf{x}, \mathbf{s}, d) + \partial\_{\mathbf{x}} \pi\_i(t, \mathbf{x}, \mathbf{s}, d) \cdot v\_i(\rho\_i(t, \mathbf{x})) = 0. \tag{10}
$$

This equation is coupled with equation (9) on each line *Ii*. More precisely, equation (10) depends on the solution of (9), while in turn at junctions the values of *π<sup>i</sup>* will determine the traffic distribution on outgoing lines as explained below.

For simplicity and without loss of generality, we assume from now on that the fluxes *fi* are all the same and we indicate them with *f* . Thus, the model for a single transmission line, consists in the system of equations:

$$\begin{cases} \rho\_t + f\left(\rho\right)\_x = 0, \\ \pi\_t + \pi\_\chi \cdot v(\rho) = 0. \end{cases}$$

To treat the evolution at junctions, let us introduce some notations. Fix a junction *J* with *n* incoming transmission lines, say *I*1, ..., *In*, and *m* outgoing transmission lines, say *In*+1, ..., *In*+*<sup>m</sup>* (junction of *n* × *m* type). The basic ingredient for the solution of Cauchy Problems by Wave Front Tracking method is the solution of Riemann Problems (RPs) (see Bressan (2000), Dafermos (1999), Serre (1999)).

We call RP for a junction the Cauchy Problem corresponding to an initial data *ρ*1,0, ..., *ρn*+*m*,0 ∈ [0, 1], and *πs*,*<sup>d</sup>* <sup>1</sup> , ...,*πs*,*<sup>d</sup> <sup>n</sup>*+*<sup>m</sup>* ∈ [0, 1] which are constant on each transmission line.

**Definition 3.** *A Riemann Solver (RS) for the junction J is a map that associates to Riemann data ρ*<sup>0</sup> = (*ρ*1,0,..., *ρn*+*m*,0) *and* Π<sup>0</sup> = (*π*1,0,..., *πn*+*m*,0) *at J the vectors ρ*ˆ = (*ρ*ˆ1,..., *ρ*ˆ*n*+*m*) *and* Πˆ = (*π*ˆ 1,..., *π*ˆ *<sup>n</sup>*+*m*) *so that the solution on an incoming transmission line Ii, i* = 1, . . . , *n, is given by the wave* (*ρi*,0, *ρ*ˆ*i*) *and on an outgoing one Ij, j* = *n* + 1, . . . , *n* + *m, is given by the waves* (*ρ*ˆ*j*, *ρj*,0) *and* (*π*ˆ*j*, *πj*,0)*. We require the following consistency condition:*

*(CC) RS*(*RS*(*ρ*0, Π0)) = *RS*(*ρ*0, Π0)*.*

Once a RS is defined and the solution of the RP is obtained, we can define admissible solutions at junctions.

#### **3. Riemann Solvers at junctions**

Consider a junction *J* of *n* × *m* type. We denote with *ρi*(*t*, *x*), *i* = 1, ..., *n* and *ρj*(*t*, *x*), *j* = *n* + 1, ..., *n* + *m* the traffic densities, respectively, on the incoming transmission lines and on the outgoing ones and by (*ρ*1,0, ..., *ρn*+*m*,0) the initial datum.

Define the maximum flux that can be obtained by a single wave solution on each transmission line as follows:

$$\gamma\_i^{\text{max}} = \begin{cases} f(\rho\_{i,0})\_\prime \text{ if } \rho\_{i,0} \in [0, \sigma]\_\prime \\ f(\sigma)\_\prime \text{ if } \rho\_{i,0} \in [\sigma, 1]\_\prime \end{cases} \\ i = 1, \ldots, n,\tag{11}$$

and

6 Telecommunications Networks

*π<sup>i</sup>* : [0, ∞[×[*ai*, *bi*] ×S×D �→[0, 1]

In other words, *πi*(*t*, *x*,*s*, *d*) specifies the density fraction *ρi*(*t*, *x*) that started from source *s* and

Assuming, on the discrete model, that a FIFO policy is used at nodes, it is natural that the averaged velocity, obtained in the limit procedure, is independent from the original sources of packets and their final destinations. In other words, we make the following hypothesis:

(H5) On each line *Ii*, the average velocity of packets depends only on the value of the density

This equation is coupled with equation (9) on each line *Ii*. More precisely, equation (10) depends on the solution of (9), while in turn at junctions the values of *π<sup>i</sup>* will determine the

For simplicity and without loss of generality, we assume from now on that the fluxes *fi* are all the same and we indicate them with *f* . Thus, the model for a single transmission line, consists

We call RP for a junction the Cauchy Problem corresponding to an initial data *ρ*1,0, ..., *ρn*+*m*,0 ∈

Once a RS is defined and the solution of the RP is obtained, we can define admissible solutions

*<sup>n</sup>*+*<sup>m</sup>* ∈ [0, 1] which are constant on each transmission line. **Definition 3.** *A Riemann Solver (RS) for the junction J is a map that associates to Riemann data ρ*<sup>0</sup> = (*ρ*1,0,..., *ρn*+*m*,0) *and* Π<sup>0</sup> = (*π*1,0,..., *πn*+*m*,0) *at J the vectors ρ*ˆ = (*ρ*ˆ1,..., *ρ*ˆ*n*+*m*) *and* Πˆ = (*π*ˆ 1,..., *π*ˆ *<sup>n</sup>*+*m*) *so that the solution on an incoming transmission line Ii, i* = 1, . . . , *n, is given by the wave* (*ρi*,0, *ρ*ˆ*i*) *and on an outgoing one Ij, j* = *n* + 1, . . . , *n* + *m, is given by the waves* (*ρ*ˆ*j*, *ρj*,0)

 *<sup>ρ</sup><sup>t</sup>* <sup>+</sup> *<sup>f</sup>* (*ρ*)*<sup>x</sup>* <sup>=</sup> 0, *π<sup>t</sup>* + *π<sup>x</sup>* · *v*(*ρ*) = 0. To treat the evolution at junctions, let us introduce some notations. Fix a junction *J* with *n* incoming transmission lines, say *I*1, ..., *In*, and *m* outgoing transmission lines, say *In*+1, ..., *In*+*<sup>m</sup>* (junction of *n* × *m* type). The basic ingredient for the solution of Cauchy Problems by Wave Front Tracking method is the solution of Riemann Problems (RPs) (see Bressan (2000),

*∂tπi*(*t*, *x*,*s*, *d*) + *∂xπi*(*t*, *x*,*s*, *d*) · *vi*(*ρi*(*t*, *x*)) = 0. (10)

*πi*(*t*, *x*,*s*, *d*) = 1.

∑*s*∈S,*d*∈D

**Definition 2.** *A traffic-type function on a line Ii is a function*

*ρ<sup>i</sup>* and not on the values of the traffic-type function *πi*.

traffic distribution on outgoing lines as explained below.

*and* (*π*ˆ*j*, *πj*,0)*. We require the following consistency condition:*

in the system of equations:

Dafermos (1999), Serre (1999)).

<sup>1</sup> , ...,*πs*,*<sup>d</sup>*

*(CC) RS*(*RS*(*ρ*0, Π0)) = *RS*(*ρ*0, Π0)*.*

[0, 1], and *πs*,*<sup>d</sup>*

at junctions.

As a consequence of hypothesis (H5), we deduce the semilinear equation

*such that, for every t* ∈ [0, ∞[ *and x* ∈ [*ai*, *bi*]

is moving towards the final destination *d*.

$$\gamma\_{j}^{\text{max}} = \begin{cases} \begin{array}{c} f(\sigma), \text{ if } \rho\_{j,0} \in [0, \sigma], \\ f(\rho\_{j,0}), \text{ if } \rho\_{j,0} \in [\sigma, 1] \end{array} \end{cases} \text{; } j = n+1, \ldots, n+m. \tag{12}$$

Finally denote with

$$\begin{aligned} \Omega\_i &= [0, \gamma\_i^{\max}]\_\prime i = 1, \dots, n, \\ \Omega\_j &= [0, \gamma\_j^{\max}]\_\prime j = n + 1, \dots, n + m\_\prime \end{aligned}$$

and with *γ*ˆ*inc* = (*f*(*ρ*ˆ*i*), ..., *f*(*ρ*ˆ*n*)), *γ*ˆ *out* = (*f*(*ρ*ˆ*n*+1), ..., *f*(*ρ*ˆ*n*+*m*)) where *ρ*ˆ = (*ρ*ˆ1, ..., *ρ*ˆ*n*+*m*) is the solution of the RP at the junction.

Now, we discuss some possible choices for the traffic distribution function:


If *rJ* is of type 1), then each packet has a deterministic route, it means that, at the junction *J*, the traffic that started at source *s* and has *d* as final destination, coming from the transmission line *Ii*, is routed on an assigned line *Ij* (*rJ*(*i*,*s*, *d*) = *j*).

Instead if *rJ* is of type 2), at the junction *J*, the traffic with source *s* and destination *d* coming from a line *Ii* is routed on every line *Ij* ∈ Out(*J*) or on some lines *Ij* ∈ Out(*J*). We can define *rJ*(*i*,*s*, *d*) in two different ways:

$$\begin{aligned} \mathbf{2a}) \ r\_{I} &: \text{Inc}(I) \times \mathcal{S} \times \mathcal{D} \hookrightarrow \text{Out}(I), \\ r\_{j}(i, s, d) &\subseteq \text{Out}(I); \\ \mathbf{2b}) \ r\_{I} &: \text{Inc}(I) \times \mathcal{S} \times \mathcal{D} \to [0, 1]^{\text{Out}(I)}, \\ r\_{I}(i, s, d) &= (a\_{I}^{i, s, d, n+1}, \dots, a\_{I}^{i, s, d, n+m}) \\ \text{with } 0 \le a\_{I}^{i, s, d, j} &\le 1, j \in \{n+1, \dots, n+m\}\_{\prime} \sum\_{j=n+1}^{n+m} a\_{I}^{i, s, d, j} = 1. \end{aligned}$$

In case 2a) we have to specify in which way the traffic at junction *J* is splitted towards the outgoing lines.

The definition 2b) means that, at the junction *J*, the traffic with source *s* and destination *d* coming from line *Ii* is routed on the outgoing line *Ij*, *j* = *n* + 1, ..., *n* + *m* with probability *α i*,*s*,*d*,*j <sup>J</sup>* .

Let us analyze how the distribution matrix *A* is constructed using *π* and *rJ*.

(B) respecting (A) the router chooses to send packets in order to maximize fluxes (i.e., the

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 429

(RA2) We assume that the number of packets through the junction is maximized both over

Ω*adm* = {*γ*ˆ : *γ*ˆ ∈ Ω<sup>1</sup> × ... × Ω*n*, ∃*A* ∈ A *t*.*c*.*Aγ*ˆ ∈ Ω*n*+<sup>1</sup> × ... × Ω*n*+*m*} ,

If the region Ω*adm* is convex than rules (A) and (B) amount to the Linear Programming

**Definition 5.** *Let τ* : [0, 1] → [0, 1] *be the map such that f*(*τ*(*ρ*)) = *f*(*ρ*) *for every ρ* ∈ [0, 1] *and*

We need some assumption on the matrix *A* (satisfied under generic conditions for *m* = *n*). Let {*e*1, ...,*en*} be the canonical basis of **<sup>R</sup>***<sup>n</sup>* and for every subset *<sup>V</sup>* <sup>⊂</sup> **<sup>R</sup>***<sup>n</sup>* indicate by *<sup>V</sup>*<sup>⊥</sup> its orthogonal. Define for every *i* = 1, ..., *n*, *Hi* = {*ei*}⊥, i.e. the coordinate hyperplane orthogonal to *ei* and for every *<sup>j</sup>* <sup>=</sup> *<sup>n</sup>* <sup>+</sup> 1, ..., *<sup>n</sup>* <sup>+</sup> *<sup>m</sup>* let *<sup>α</sup><sup>j</sup>* <sup>=</sup> {*αj*1, ..., *<sup>α</sup>jn*} ∈ **<sup>R</sup>***<sup>n</sup>* and define *Hj* <sup>=</sup> {*αj*}⊥. Let K be the set of indices *k* = (*k*1, ..., *kl*), 1 ≤ *l* ≤ *n* − 1, such that 0 ≤ *k*<sup>1</sup> < *k*<sup>2</sup> < ... < *kl* ≤ *n* + *m*

**Theorem 6.** *(Theorem 3.1 in Coclite et al. (2005) and 3.2 in Garavello et al. (2005)) Let* (*N*, I*,* F*,* J *,* S*,* D*,* R) *be an admissible network and J a junction of n* × *m type. Assume that the flux f* : [0, 1] → **R** *satisfies (F) and the matrix A satisfies condition (C). For every ρ*1,0, ..., *ρn*+*m*,0 ∈ [0, 1]*, and for*

<sup>1</sup> , ...,*πn*+*m*(0, ·,*s*, *<sup>d</sup>*) = *<sup>π</sup>s*,*<sup>d</sup>*

*ρ*1(0, ·) ≡ *ρ*1,0, ..., *ρn*+*m*(0, ·) ≡ *ρn*+*m*,0,

*<sup>τ</sup>*(*ρi*,0), 1

*Hh*. Letting **<sup>1</sup>** = (1, ..., 1) <sup>∈</sup> **<sup>R</sup>***n*, we assume

*<sup>n</sup>*+*<sup>m</sup>* ∈ [0, 1] *, there exist densities ρ*ˆ1, ..., *ρ*ˆ*n*+*<sup>m</sup> and a unique admissible centered weak*

, *if* 0 ≤ *ρi*,0 ≤ *σ*,

[*σ*, 1] , *if <sup>σ</sup>* <sup>≤</sup> *<sup>ρ</sup>i*,0 <sup>≤</sup> 1, *<sup>i</sup>* <sup>=</sup> 1, ..., *<sup>n</sup>*, (14)

*<sup>n</sup>*+*m*,(*s* ∈ S, *d* ∈ D).

(*γ*ˆ1 + *γ*ˆ2).

In case 2a) first we observe that the set A is convex. The admissible region given by

max *γ*ˆ∈Ω*adm*

number of packets which are processed).

is convex at least for the case of junctions of 2 × 2.

This problem has clearly a solution, which may not be unique. Let us consider the case 2b). We need some more notations.

> *l h*=1

*k* .

incoming and outgoing lines.

We have to distinguish case 2a) and 2b).

*τ*(*ρ*) �= *ρ for every ρ* ∈ [0, 1]\{*σ*}*.*

and for every *<sup>k</sup>* ∈ K set *Hk* =

In case 2b) the following result holds

*solution, ρ* = (*ρ*1, ..., *ρn*+*m*) *at J such that*

*<sup>π</sup>*1(0, ·*s*, *<sup>d</sup>*) = *<sup>π</sup>s*,*<sup>d</sup>*

{*ρi*,0} ∪

*ρ*ˆ*<sup>i</sup>* ∈

(C) for every *k* ∈ K, **1** ∈/*H*<sup>⊥</sup>

<sup>1</sup> , ...*πs*,*<sup>d</sup>*

*every πs*,*<sup>d</sup>*

*We have*

**3.1 Algorithm (RA1)**

problem:

**Definition 4.** *A distribution matrix is a matrix*

$$A \doteq \left\{ \alpha\_{j,i} \right\}\_{j=n+1,\ldots,n+m, i=1,\ldots,n} \in \mathbb{R}^{m \times m}$$

*such that*

$$0 < \alpha\_{j,i} < 1, \quad \sum\_{j=n+1}^{n+m} \alpha\_{j,i} = 1,$$

*for each i* = 1, ..., *n and j* = *n* + 1, ..., *n* + *m, where αj*,*<sup>i</sup> is the percentage of packets arriving from the i-th incoming transmission line that take the j-th outgoing transmission line.*

In case 1) we can define the matrix *A* in the following way. Fix a time *t* and assume that for all *i* ∈ Inc(*J*),*s* ∈ S and *d* ∈ D, *πi*(*t*, ·,*s*, *d*) admits a limit at the junction *J*, i.e left limit at *bi*. For *i* ∈ {1, ..., *n*}, *j* ∈ {*n* + 1, ..., *n* + *m*}, we set

$$\mathfrak{a}\_{j,\dot{i}} = \sum\_{\substack{\mathbf{s} \in \mathcal{S}, \mathbf{d} \in \mathcal{D}, \\ r\_{\mathbf{l}}(i, \mathbf{s}, \mathbf{d}) = \dot{j}}} \pi\_{\mathbf{i}}(t, \mathbf{b}\_{\mathbf{i}} - \mathbf{s}, \mathbf{d}).$$

The fluxes *fi*(*ρi*) to be consistent with the traffic-type functions must satisfy the following relation:

$$f\_{\vec{l}}(\rho\_{\vec{l}}(\cdot, a\_{\vec{l}} +)) = \sum\_{i=1}^{n} a\_{\vec{l},i} f\_{\vec{l}}(\rho\_{\vec{l}}(\cdot, b\_{\vec{l}} -))\_{\vec{l}}$$

for each *j* = *n* + 1, ..., *n* + *m*.

Let us analyze how to define the matrix *A* in the case 2a). We may assign *ϕ*(*i*,*s*, *d*) ∈ *rJ*(*i*,*s*, *d*) and set

$$\begin{aligned} \mathfrak{a}\_{j,i} &= \sum\_{\substack{s \in \mathcal{S}, d \in \mathcal{D}\_{\prime} \\ i \cdot \wp(i, s, d) = j}} \pi\_i(t, b\_i - , s, d), \\ \mathfrak{a}\_{j,i} &= 0, \text{ if } j \notin r\_f(i, s, d). \end{aligned}$$

However, it is more natural to assign a flexible strategy defining a set of admissible matrices *A* in the following way

$$\mathcal{A} = \left\{ \begin{aligned} A: & \exists a\_{\boldsymbol{I}}^{i,s,d,j} \in [0,1], \sum\_{j=n+1}^{n+m} a\_{\boldsymbol{I}}^{i,s,d,j} = 1, a\_{\boldsymbol{I}}^{i,s,d,j} = 0, \text{ if } j \notin r\_{\boldsymbol{I}}(i,s,d): \\ & \mathfrak{a}\_{\boldsymbol{j},i} = \sum\_{\substack{s \in \mathcal{S}, d \in \mathcal{D}, \\ j \in r\_{\boldsymbol{I}}(i,s,d)}} \pi\_{\boldsymbol{i}}(\mathfrak{t}, b\_{\boldsymbol{i}}-,s,d)a\_{\boldsymbol{I}}^{i,s,d,j} \end{aligned} \right\}.$$

Finally, we treat now the case 2b). In this case the matrix *A* is unique and is defined by

$$a\_{\mathbf{j},\mathbf{i}} = \sum\_{s \in \mathcal{S}, d \in \mathcal{D}} \pi\_{\mathbf{i}}(\mathbf{t}, \mathbf{b}\_{\mathbf{i}} - \mathbf{s}, \mathbf{d}) a\_{\mathbf{j}}^{\mathbf{i}, \mathbf{s}, \mathbf{d}, \mathbf{j}}.\tag{13}$$

We describe two different RSs at a junction that represent two different routing algorithms:

(RA1) We assume that

(A) the traffic from incoming transmission lines is distributed on outgoing transmission lines according to fixed coefficients;


#### **3.1 Algorithm (RA1)**

8 Telecommunications Networks

*<sup>j</sup>*=*n*+1,...,*n*+*m*,*i*=1,...,*<sup>n</sup>* <sup>∈</sup> **<sup>R</sup>***m*×*<sup>n</sup>*

*αj*,*<sup>i</sup>* = 1,

*πi*(*t*, *bi*−,*s*, *d*).

*αj*,*<sup>i</sup> fi*(*ρi*(·, *bi*−)),

*πi*(*t*, *bi*−,*s*, *d*),

*i*,*s*,*d*,*j*

*πi*(*t*, *bi*−,*s*, *d*)*α*

*πi*(*t*, *bi*−,*s*, *d*)*α*

*<sup>J</sup>* = 0, if *j* ∈/ *rJ*(*i*,*s*, *d*) :

⎫ ⎪⎪⎪⎬

⎪⎪⎪⎭ .

*<sup>J</sup>* . (13)

*i*,*s*,*d*,*j J*

*i*,*s*,*d*,*j*

*n*+*m* ∑ *j*=*n*+1

*for each i* = 1, ..., *n and j* = *n* + 1, ..., *n* + *m, where αj*,*<sup>i</sup> is the percentage of packets arriving from the*

In case 1) we can define the matrix *A* in the following way. Fix a time *t* and assume that for all *i* ∈ Inc(*J*),*s* ∈ S and *d* ∈ D, *πi*(*t*, ·,*s*, *d*) admits a limit at the junction *J*, i.e left limit at *bi*. For

The fluxes *fi*(*ρi*) to be consistent with the traffic-type functions must satisfy the following

*n* ∑ *i*=1

Let us analyze how to define the matrix *A* in the case 2a). We may assign *ϕ*(*i*,*s*, *d*) ∈ *rJ*(*i*,*s*, *d*)

**Definition 4.** *A distribution matrix is a matrix*

*i* ∈ {1, ..., *n*}, *j* ∈ {*n* + 1, ..., *n* + *m*}, we set

*such that*

relation:

and set

for each *j* = *n* + 1, ..., *n* + *m*.

*A* in the following way

A =

(RA1) We assume that

according to fixed coefficients;

⎧ ⎪⎪⎪⎨

*A* : ∃*α*

*i*,*s*,*d*,*j <sup>J</sup>* ∈ [0, 1],

⎪⎪⎪⎩

*A*=˙ � *αj*,*<sup>i</sup>* �

0 < *αj*,*<sup>i</sup>* < 1,

*αj*,*<sup>i</sup>* = ∑

*fj*(*ρj*(·, *aj*+)) =

*αj*,*<sup>i</sup>* = ∑

*n*+*m* ∑ *j*=*n*+1 *α i*,*s*,*d*,*j <sup>J</sup>* = 1, *α*

*<sup>α</sup>j*,*<sup>i</sup>* <sup>=</sup> <sup>∑</sup> *<sup>s</sup>*∈S,*d*∈D, *j*∈*rJ*(*i*,*s*,*d*)

*<sup>α</sup>j*,*<sup>i</sup>* <sup>=</sup> <sup>∑</sup>*<sup>s</sup>*∈S,*d*∈D

*s*∈S,*d*∈D, *i*:*ϕ*(*i*,*s*,*d*)=*j*

*αj*,*<sup>i</sup>* = 0, if *j* ∈/ *rJ*(*i*,*s*, *d*).

Finally, we treat now the case 2b). In this case the matrix *A* is unique and is defined by

We describe two different RSs at a junction that represent two different routing algorithms:

(A) the traffic from incoming transmission lines is distributed on outgoing transmission lines

However, it is more natural to assign a flexible strategy defining a set of admissible matrices

*s*∈S,*d*∈D, *rJ*(*i*,*s*,*d*)=*j*

*i-th incoming transmission line that take the j-th outgoing transmission line.*

We have to distinguish case 2a) and 2b).

In case 2a) first we observe that the set A is convex. The admissible region given by

$$\Omega\_{adm} = \{ \gamma : \gamma \in \Omega\_1 \times \dots \times \Omega\_{\mathfrak{n}}, \exists A \in \mathcal{A} \,\, t.c.A \,\gamma \in \Omega\_{\mathfrak{n}+1} \times \dots \times \Omega\_{\mathfrak{n}+\mathfrak{m}} \} \,\,.$$

is convex at least for the case of junctions of 2 × 2.

If the region Ω*adm* is convex than rules (A) and (B) amount to the Linear Programming problem:

$$\max\_{\hat{\gamma}\in\Omega\_{atm}} (\hat{\gamma}\_1 + \hat{\gamma}\_2).$$

This problem has clearly a solution, which may not be unique.

Let us consider the case 2b). We need some more notations.

**Definition 5.** *Let τ* : [0, 1] → [0, 1] *be the map such that f*(*τ*(*ρ*)) = *f*(*ρ*) *for every ρ* ∈ [0, 1] *and τ*(*ρ*) �= *ρ for every ρ* ∈ [0, 1]\{*σ*}*.*

We need some assumption on the matrix *A* (satisfied under generic conditions for *m* = *n*). Let {*e*1, ...,*en*} be the canonical basis of **<sup>R</sup>***<sup>n</sup>* and for every subset *<sup>V</sup>* <sup>⊂</sup> **<sup>R</sup>***<sup>n</sup>* indicate by *<sup>V</sup>*<sup>⊥</sup> its orthogonal. Define for every *i* = 1, ..., *n*, *Hi* = {*ei*}⊥, i.e. the coordinate hyperplane orthogonal to *ei* and for every *<sup>j</sup>* <sup>=</sup> *<sup>n</sup>* <sup>+</sup> 1, ..., *<sup>n</sup>* <sup>+</sup> *<sup>m</sup>* let *<sup>α</sup><sup>j</sup>* <sup>=</sup> {*αj*1, ..., *<sup>α</sup>jn*} ∈ **<sup>R</sup>***<sup>n</sup>* and define *Hj* <sup>=</sup> {*αj*}⊥. Let K be the set of indices *k* = (*k*1, ..., *kl*), 1 ≤ *l* ≤ *n* − 1, such that 0 ≤ *k*<sup>1</sup> < *k*<sup>2</sup> < ... < *kl* ≤ *n* + *m* and for every *<sup>k</sup>* ∈ K set *Hk* = *l h*=1 *Hh*. Letting **<sup>1</sup>** = (1, ..., 1) <sup>∈</sup> **<sup>R</sup>***n*, we assume

(C) for every *k* ∈ K, **1** ∈/*H*<sup>⊥</sup> *k* .

In case 2b) the following result holds

**Theorem 6.** *(Theorem 3.1 in Coclite et al. (2005) and 3.2 in Garavello et al. (2005)) Let* (*N*, I*,* F*,* J *,* S*,* D*,* R) *be an admissible network and J a junction of n* × *m type. Assume that the flux f* : [0, 1] → **R** *satisfies (F) and the matrix A satisfies condition (C). For every ρ*1,0, ..., *ρn*+*m*,0 ∈ [0, 1]*, and for every πs*,*<sup>d</sup>* <sup>1</sup> , ...*πs*,*<sup>d</sup> <sup>n</sup>*+*<sup>m</sup>* ∈ [0, 1] *, there exist densities ρ*ˆ1, ..., *ρ*ˆ*n*+*<sup>m</sup> and a unique admissible centered weak solution, ρ* = (*ρ*1, ..., *ρn*+*m*) *at J such that*

$$\begin{aligned} \rho\_1(0, \cdot) &\equiv \rho\_{1,0\prime} \dots \rho\_{n+m}(0, \cdot) \equiv \rho\_{n+m,0\prime} \\ \pi^1(0, \cdot s, d) &= \pi\_1^{s,d} \dots \pi^{n+m}(0, \cdot, s, d) = \pi\_{n+m\prime}^{s,d} (s \in \mathcal{S}, d \in \mathcal{D}). \end{aligned}$$

*We have*

$$\rho\_i \in \begin{cases} \{\rho\_{i,0}\} \cup \left[\tau(\rho\_{i,0}), 1\right], \text{ if } 0 \le \rho\_{i,0} \le \sigma\_\prime \\ \left[\sigma, 1\right], & \text{if } \sigma \le \rho\_{i,0} \le 1 \end{cases} \tag{14} \\ \text{ o } = 1, \dots, n, \tag{14}$$

*rq*

*<sup>j</sup>* , *j* = 3, 4. Let us determine *γ*ˆ*<sup>j</sup>* in the second case, using the traffic

Γ1

<sup>P</sup> *rq* P Q

Γ1 max

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 431

In the first case we set (*γ*ˆ1, *γ*ˆ2) = *P*, while in the second case we set (*γ*ˆ1, *γ*ˆ2) = *Q*, with

distribution parameter *α*. Since not all packets can go on the outgoing transmission lines, we let *C* be the amount that goes through. Then *αC* packets go on the outgoing line *I*<sup>3</sup> and (1 − *α*)*C* on the outgoing line *I*4. Consider the space (*γ*3, *γ*4) and define the following lines:

*<sup>r</sup><sup>α</sup>* : *<sup>γ</sup>*<sup>4</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>α</sup>*

*r*<sup>Γ</sup> : *γ*<sup>3</sup> + *γ*<sup>4</sup> = Γ.

and *G*<sup>1</sup> and *G*<sup>2</sup> its endpoints. Since in case 2a) we have an infinite number of matrices *A*, each of one determines a line *rα*, we choose the most "natural" line *r<sup>α</sup>* , i.e. the one nearest to the

(*γ*3, *<sup>γ</sup>*4) : 0 <sup>≤</sup> *<sup>γ</sup><sup>j</sup>* <sup>≤</sup> *<sup>γ</sup>*max

<sup>4</sup> ), *<sup>Q</sup>* = (*γ*max

*<sup>α</sup> <sup>γ</sup>*3,

*inc* : *<sup>A</sup>* ∈ A

,

*<sup>j</sup>* , *<sup>j</sup>* <sup>=</sup> 3, 4

<sup>3</sup> , <sup>Γ</sup> <sup>−</sup> *<sup>γ</sup>*max

.

<sup>3</sup> ). We distinguish 3 cases:

*<sup>Q</sup>* = *proj*Ω*in*∩*r*<sup>Γ</sup> (*P*) where *proj* is the usual projection on a convex set, see Figure 2.

*<sup>i</sup>*,0 , *i* = 1, 2.

*<sup>i</sup>* <sup>=</sup> *<sup>π</sup>s*,*<sup>d</sup>*

Let us now determine *γ*ˆ*j*, *j* = 3, 4. We have to distinguish again two cases :

We have to distinguish case 2a) and 2b) for the traffic distribution function.

G = *Aγ*ˆ *<sup>T</sup>*

statistic line determined by measurements on the network.

Recall that the final fluxes should belong to the region:

Ω*out* =

<sup>4</sup> , *<sup>γ</sup>*max

*r*�

Γ2

Γ2 max

Fig. 2. *P* belongs to Ω*in* and *P* is outside Ω*in*.

b) *P* is outside Ω*in*.

**I** Γmax *out* = Γ,

**II** Γmax

*out* > Γ.

**3.2.1 Case 2a)**

In the first case *γ*ˆ*<sup>j</sup>* = *γ*max

Let us introduce the connected set

Define *<sup>P</sup>* <sup>=</sup> *<sup>r</sup><sup>α</sup>* <sup>∩</sup> *<sup>r</sup>*Γ, *<sup>R</sup>* = (<sup>Γ</sup> <sup>−</sup> *<sup>γ</sup>*max

As for the algorithm (RA1) *π*ˆ*s*,*<sup>d</sup>*

$$\hat{\rho}\_{j} \in \begin{cases} [0, \sigma]\_{\prime} & \text{if } 0 \le \rho\_{j,0} \le \sigma \\ \{\rho\_{j,0}\} \cup \left[0, \pi(\rho\_{j,0})\right]\_{\prime} & \text{if } \sigma \le \rho\_{j,0} \le 1 \end{cases} \quad j = n+1, \ldots, n+m \tag{15}$$

*and on each incoming line Ii, i* = 1, ..., *n, the solution consists of the single wave* (*ρi*,0, *ρ*ˆ*i*)*, while on each outgoing line Ij*, *j* = *n* + 1, ..., *n* + *m, the solution consists of the single wave* (*ρ*ˆ*j*, *ρj*,0)*. Moreover <sup>π</sup>*ˆ*i*(*t*, ·,*s*, *<sup>d</sup>*) = *<sup>π</sup>s*,*<sup>d</sup> <sup>i</sup> for every t* ≥ 0, *i* ∈ {1, ..., *n*},*s* ∈ S, *d* ∈ D *and*

$$\hat{\pi}\_{\dot{\jmath}}(t, a\_{\dot{\jmath}}+, s, d) = \frac{\sum\_{i=1}^{n} \alpha\_{f}^{i, s, d, j} \pi\_{i}^{s, d}(t, b\_{i}-, s, d) f(\hat{\rho}\_{i})}{f(\hat{\rho}\_{\dot{\jmath}})}$$

*for every t* ≥ 0, *j* ∈ {*n* + 1, ..., *n* + *m*},*s* ∈ S, *d* ∈ D.

#### **3.2 Algorithm (RA2)**

To solve RPs according to (RA2) we need some additional parameters called priority and traffic distribution parameters. For simplicity of exposition, consider, junction *J* of 2 × 2 type. In this case we have only one priority parameter *q* ∈ ]0, 1[ and one traffic distribution parameter *<sup>α</sup>* <sup>∈</sup> ]0, 1[. We denote with (*ρ*1,0, *<sup>ρ</sup>*2,0, *<sup>ρ</sup>*3,0, *<sup>ρ</sup>*4,0) and (*πs*,*<sup>d</sup>* 1,0, *<sup>π</sup>s*,*<sup>d</sup>* 2,0, *<sup>π</sup>s*,*<sup>d</sup>* 3,0, *<sup>π</sup>s*,*<sup>d</sup>* 4,0) the initial data.

In order to maximize the number of packets through the junction over incoming and outgoing lines we define

$$
\Gamma = \min \left\{ \Gamma\_{in}^{\max} , \Gamma\_{out}^{\max} \right\} \,\,\,\,
$$

where Γmax *in* <sup>=</sup> *<sup>γ</sup>*max <sup>1</sup> <sup>+</sup> *<sup>γ</sup>*max <sup>2</sup> and <sup>Γ</sup>max *out* = *<sup>γ</sup>*max <sup>3</sup> <sup>+</sup> *<sup>γ</sup>*max <sup>4</sup> . Thus we want to have Γ as flux through the junction.

One easily see that to solve the RP, it is enough to determine the fluxes *γ*ˆ*<sup>i</sup>* = *f*(*ρ*ˆ*i*), *i* = 1, 2. In fact, to have simple waves with the appropriate velocities, i.e. negative on incoming lines and positive on outgoing ones, we get the constraints (14), (15). Observe that we compute *γ*ˆ*<sup>i</sup>* = *f*(*ρ*ˆ*i*), *i* = 1, 2 without taking into account the type of traffic distribution function.

We have to distinguish two cases:

$$\begin{array}{c} \mathbf{I} \quad \Gamma\_{in}^{\max} = \Gamma\_{\prime} \\ \mathbf{II} \quad \Gamma\_{in}^{\max} > \Gamma\_{\prime} \end{array}$$

In the first case we set *γ*ˆ*<sup>i</sup>* = *γ*max *<sup>i</sup>* , *i* = 1, 2.

Let us analyze the second case in which we use the priority parameter *q*. Not all packets can enter the junction, so let *C* be the amount of packets that can go through: *qC* packets come from first incoming line and (1 − *q*)*C* packets from the second. In the space (*γ*1, *γ*2), define the following lines:

$$r\_{\eta}:\gamma\_2 = \frac{1-q}{q}\gamma\_{1\prime} \quad r\_{\Gamma}:\gamma\_1 + \gamma\_2 = \Gamma\_{\prime}$$

and *P* the point of intersection of *rq* and *r*Γ. Recall that the final fluxes should belong to the region:

$$\Omega\_{\rm in} = \left\{ (\gamma\_1, \gamma\_2) : 0 \le \gamma\_i \le \gamma\_i^{\rm max}, i = 1, 2 \right\}.$$

We distinguish two cases:

a) *P* belongs to Ω*in*,

Fig. 2. *P* belongs to Ω*in* and *P* is outside Ω*in*.

b) *P* is outside Ω*in*.

10 Telecommunications Networks

*and on each incoming line Ii, i* = 1, ..., *n, the solution consists of the single wave* (*ρi*,0, *ρ*ˆ*i*)*, while on each outgoing line Ij*, *j* = *n* + 1, ..., *n* + *m, the solution consists of the single wave* (*ρ*ˆ*j*, *ρj*,0)*. Moreover*

To solve RPs according to (RA2) we need some additional parameters called priority and traffic distribution parameters. For simplicity of exposition, consider, junction *J* of 2 × 2 type. In this case we have only one priority parameter *q* ∈ ]0, 1[ and one traffic distribution

In order to maximize the number of packets through the junction over incoming and outgoing

<sup>3</sup> <sup>+</sup> *<sup>γ</sup>*max

One easily see that to solve the RP, it is enough to determine the fluxes *γ*ˆ*<sup>i</sup>* = *f*(*ρ*ˆ*i*), *i* = 1, 2. In fact, to have simple waves with the appropriate velocities, i.e. negative on incoming lines and positive on outgoing ones, we get the constraints (14), (15). Observe that we compute *γ*ˆ*<sup>i</sup>* = *f*(*ρ*ˆ*i*), *i* = 1, 2 without taking into account the type of traffic distribution function.

Let us analyze the second case in which we use the priority parameter *q*. Not all packets can enter the junction, so let *C* be the amount of packets that can go through: *qC* packets come from first incoming line and (1 − *q*)*C* packets from the second. In the space (*γ*1, *γ*2), define

and *P* the point of intersection of *rq* and *r*Γ. Recall that the final fluxes should belong to the

<sup>Ω</sup>*in* <sup>=</sup> {(*γ*1, *<sup>γ</sup>*2) : 0 <sup>≤</sup> *<sup>γ</sup><sup>i</sup>* <sup>≤</sup> *<sup>γ</sup>*max

*<sup>q</sup> <sup>γ</sup>*1, *<sup>r</sup>*<sup>Γ</sup> : *<sup>γ</sup>*<sup>1</sup> <sup>+</sup> *<sup>γ</sup>*<sup>2</sup> <sup>=</sup> <sup>Γ</sup>,

*<sup>i</sup>* , *i* = 1, 2} .

*in* , <sup>Γ</sup>max *out* } ,

<sup>Γ</sup> <sup>=</sup> min {Γmax

*out* = *<sup>γ</sup>*max

*<sup>i</sup>* , *i* = 1, 2.

*rq* : *<sup>γ</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>q</sup>*

, *if <sup>σ</sup>* <sup>≤</sup> *<sup>ρ</sup>j*,0 <sup>≤</sup> 1, *<sup>j</sup>* <sup>=</sup> *<sup>n</sup>* <sup>+</sup> 1, ..., *<sup>n</sup>* <sup>+</sup> *<sup>m</sup>*, (15)

1,0, *<sup>π</sup>s*,*<sup>d</sup>*

2,0, *<sup>π</sup>s*,*<sup>d</sup>*

<sup>4</sup> . Thus we want to have Γ as flux through

3,0, *<sup>π</sup>s*,*<sup>d</sup>*

4,0) the initial

*<sup>i</sup>* (*t*, *bi*−,*s*, *d*)*f*(*ρ*ˆ*i*)

*f*(*ρ*ˆ*j*)

[0, *<sup>σ</sup>*], *if* <sup>0</sup> <sup>≤</sup> *<sup>ρ</sup>j*,0 <sup>≤</sup> *<sup>σ</sup>*,

*<sup>i</sup> for every t* ≥ 0, *i* ∈ {1, ..., *n*},*s* ∈ S, *d* ∈ D *and*

*n* ∑ *i*=1 *α i*,*s*,*d*,*j <sup>J</sup> <sup>π</sup>s*,*<sup>d</sup>*

0, *τ*(*ρj*,0)

*ρ*ˆ*<sup>j</sup>* ∈

*<sup>π</sup>*ˆ*i*(*t*, ·,*s*, *<sup>d</sup>*) = *<sup>π</sup>s*,*<sup>d</sup>*

**3.2 Algorithm (RA2)**

data.

lines we define

*in* <sup>=</sup> *<sup>γ</sup>*max

We have to distinguish two cases:

In the first case we set *γ*ˆ*<sup>i</sup>* = *γ*max

<sup>1</sup> <sup>+</sup> *<sup>γ</sup>*max

where Γmax

the junction.

**I** Γmax *in* = Γ,

**II** Γmax

region:

*in* > Γ.

the following lines:

We distinguish two cases:

a) *P* belongs to Ω*in*,

{*ρj*,0} ∪

*for every t* ≥ 0, *j* ∈ {*n* + 1, ..., *n* + *m*},*s* ∈ S, *d* ∈ D.

*π*ˆ*j*(*t*, *aj*+,*s*, *d*) =

parameter *<sup>α</sup>* <sup>∈</sup> ]0, 1[. We denote with (*ρ*1,0, *<sup>ρ</sup>*2,0, *<sup>ρ</sup>*3,0, *<sup>ρ</sup>*4,0) and (*πs*,*<sup>d</sup>*

<sup>2</sup> and <sup>Γ</sup>max

In the first case we set (*γ*ˆ1, *γ*ˆ2) = *P*, while in the second case we set (*γ*ˆ1, *γ*ˆ2) = *Q*, with *<sup>Q</sup>* = *proj*Ω*in*∩*r*<sup>Γ</sup> (*P*) where *proj* is the usual projection on a convex set, see Figure 2.

As for the algorithm (RA1) *π*ˆ*s*,*<sup>d</sup> <sup>i</sup>* <sup>=</sup> *<sup>π</sup>s*,*<sup>d</sup> <sup>i</sup>*,0 , *i* = 1, 2.

Let us now determine *γ*ˆ*j*, *j* = 3, 4. We have to distinguish again two cases :

**I** Γmax *out* = Γ, **II** Γmax *out* > Γ.

In the first case *γ*ˆ*<sup>j</sup>* = *γ*max *<sup>j</sup>* , *j* = 3, 4. Let us determine *γ*ˆ*<sup>j</sup>* in the second case, using the traffic distribution parameter *α*. Since not all packets can go on the outgoing transmission lines, we let *C* be the amount that goes through. Then *αC* packets go on the outgoing line *I*<sup>3</sup> and (1 − *α*)*C* on the outgoing line *I*4. Consider the space (*γ*3, *γ*4) and define the following lines:

$$r\_{\alpha}: \gamma\_4 = \frac{1-\alpha}{\alpha}\gamma\_3.$$

$$r\_{\Gamma}: \gamma\_3 + \gamma\_4 = \Gamma.$$

We have to distinguish case 2a) and 2b) for the traffic distribution function.

#### **3.2.1 Case 2a)**

Let us introduce the connected set

$$\mathcal{G} = \left\{ A \hat{\gamma}\_{inc}^T : A \in \mathcal{A} \right\} \ \_{\prime}$$

and *G*<sup>1</sup> and *G*<sup>2</sup> its endpoints. Since in case 2a) we have an infinite number of matrices *A*, each of one determines a line *rα*, we choose the most "natural" line *r<sup>α</sup>* , i.e. the one nearest to the statistic line determined by measurements on the network.

Recall that the final fluxes should belong to the region:

$$\Omega\_{\rm out} = \left\{ (\gamma\_3, \gamma\_4) : 0 \le \gamma\_j \le \gamma\_j^{\rm max}, j = 3, 4 \right\} \dots$$

Define *<sup>P</sup>* <sup>=</sup> *<sup>r</sup><sup>α</sup>* <sup>∩</sup> *<sup>r</sup>*Γ, *<sup>R</sup>* = (<sup>Γ</sup> <sup>−</sup> *<sup>γ</sup>*max <sup>4</sup> , *<sup>γ</sup>*max <sup>4</sup> ), *<sup>Q</sup>* = (*γ*max <sup>3</sup> , <sup>Γ</sup> <sup>−</sup> *<sup>γ</sup>*max <sup>3</sup> ). We distinguish 3 cases:

We analyze some models used in literature to evaluate the packets loss rate with the aim to

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 433

Let us consider the transmission of two consecutive routers. The node that transmits packets is called *sender*, while the receiving one is said *receiver*. Among the nodes, there is a link or channel, with limited capacity. Assume that the sender and the receiver are synchronized each other, i.e. the receiver is able to process in real time all packets, sent by the sender. In few words, no packets are lost. The packets loss can occur only on the link, due to its finite capacity. Under the zero buffer hypotheses the loss rate is defined as the proportional excess of offered traffic over the available capacity. If *R* is the sender bit rate and *C* is the link capacity, we have a loss if *R* > *C*. The model is said *proportional-excess* or briefly *P*/*E* and suppose

> 0, *R* < *C*, *R*−*C*

In Figure 3, loss probability for *P*/*E* model (continuous curve) and FD model (dashed curve) are shown, assuming *C* = *σ* = 1/2. For values *C* < *ρ* < 2*C*, the FD model overestimates the

*0.1 0.3 0.7 0.9* <sup>Ρ</sup>

Observe that the *P*/*E* model is not realistic. In fact, the sender and the receiver are never synchronized each other and whatever transmission protocol is used by the transport layer, the receiver has a finite length buffer, where the packets wait to be processed and eventually sent to the next node. Thus queueing models are needed, to infer about network performance.

Queueing models are good at predicting loss in a network with many independent users, probably using different applications. Consider the traffic from TCP sources that send packets through a bottleneck link. The traffic is aggregated and used as an arrival process for the link. The arrival process, being the aggregation of independent sources, is approximated as a Poisson process, and the aggregated throughput is used as the rate of the Poisson process (see Wierman et al. (2003)). These considerations justify the assumption that the times between the packets arrivals are exponentially distributed. Depending on the hypothesis on the length of

*C*�Σ Ρ*max*

*<sup>C</sup>* , *<sup>R</sup>* <sup>&</sup>gt; *<sup>C</sup>*. (17)

*p* =

*0.2 0.4 0.6 0.8 1 <sup>p</sup>* <sup>Ρ</sup>

Fig. 3. Loss probabilities. Dashed line: FD model. Continuous line: P/E model.

compare its behaviour with the function depicted in Figure 1.

**4.1.1 The proportional-excess model**

deterministic arrivals. The packets bit rate is:

loss probability.

**4.1.2 Models with finite capacity**

a) G ∩ Ω*out* ∩ *r*<sup>Γ</sup> �= ∅,

b) G ∩ Ω*out* ∩ *r*<sup>Γ</sup> = ∅ and *γ*3(*G*1) < *γ*3(*R*),

c) G ∩ <sup>Ω</sup>*out* <sup>∩</sup> *<sup>r</sup>*<sup>Γ</sup> <sup>=</sup> <sup>∅</sup> and *<sup>γ</sup>*3(*G*1) <sup>&</sup>gt; *<sup>γ</sup>*max <sup>3</sup> .

If the set G has a priority over the line *r*<sup>Γ</sup> we set (*γ*ˆ3, *γ*ˆ4) in the following way. In case a) we define (*γ*ˆ3, *<sup>γ</sup>*ˆ4) = *proj*G∩Ω*out*∩*r*<sup>Γ</sup> (*P*), in case b) (*γ*ˆ3, *<sup>γ</sup>*ˆ4) = *<sup>R</sup>*, and finally in case c) (*γ*ˆ3, *<sup>γ</sup>*ˆ4) = *<sup>Q</sup>*. Otherwise, if *r*<sup>Γ</sup> has a priority over G we set (*γ*ˆ3, *γ*ˆ4) = min *γ*∈Ω*out* F(*γ*,*rα*, G) where F is a convex functional which depends on *γ*, *r<sup>α</sup>* and on the set G of the routing standards.

The vector *π*ˆ *<sup>s</sup>*,*<sup>d</sup> <sup>i</sup>* , *j* = 3, 4 are computed in the same way as for the algorithm (RA1).

#### **3.2.2 Case 2b)**

In case 2b) we have a unique matrix *A*. The fluxes on outgoing lines are computed as in the case without sources and destinations.

We distinguish two cases:


In the first case we set (*γ*ˆ3, *γ*ˆ4) = *P*, while in the second case we set (*γ*ˆ3, *γ*ˆ4) = *Q*, where *Q* = *proj*Ω*adm* (*P*). Again, we can extend to the case of *m* outgoing lines.

Finally we define *π*ˆ *<sup>s</sup>*,*<sup>d</sup> <sup>i</sup>* , *j* = 3, 4 as in the case 2a):

$$\hat{\pi}\_{\dot{f}}(t, a\_{\dot{f}}+, s, d) = \frac{\sum\_{i=1}^{n} \alpha\_{f}^{i, s, d, j} \pi\_{i}^{s, d}(t, b\_{i}-, s, d) f(\beta\_{i})}{f(\hat{\rho}\_{\dot{f}})}$$

for every *t* ≥ 0, *j* ∈ {*n* + 1, ..., *n* + *m*},*s* ∈ S, *d* ∈ D.

Once solutions to RPs are given, one can use a Wave Front Tracking algorithm to construct a sequence of approximate solutions.

#### **4. Model assumptions**

The aim of this section is to verify that the assumptions underlying the data networks fluid-dynamic model (shortly FD model) are correct. Here we focus on the fixed-point models to describe TCP, and considering various set-ups with TCP traffic in a single bottleneck topology, we investigate queueing models for estimating packet loss rate. In what follows we suppose *ρmax* = 1 and *σ* = <sup>1</sup> 2 .

#### **4.1 Loss probability function**

It is reasonable to assume that the loss probability function *p* is null for some interval, which is a right neighborhood of zero. This means that at low densities no packet is lost. Then *p* should be increasing, reaching the value 1 at the maximal density, the situation of complete stuck. With the above assumptions the loss probability function in (4) can be written as:

$$p\left(\rho\right) = \begin{cases} 0, & 0 \le \rho \le 1/2, \\ \frac{2\rho - 1}{\rho}, & 1/2 \le \rho \le 1. \end{cases} \tag{16}$$

We analyze some models used in literature to evaluate the packets loss rate with the aim to compare its behaviour with the function depicted in Figure 1.

#### **4.1.1 The proportional-excess model**

12 Telecommunications Networks

If the set G has a priority over the line *r*<sup>Γ</sup> we set (*γ*ˆ3, *γ*ˆ4) in the following way. In case a) we define (*γ*ˆ3, *<sup>γ</sup>*ˆ4) = *proj*G∩Ω*out*∩*r*<sup>Γ</sup> (*P*), in case b) (*γ*ˆ3, *<sup>γ</sup>*ˆ4) = *<sup>R</sup>*, and finally in case c) (*γ*ˆ3, *<sup>γ</sup>*ˆ4) = *<sup>Q</sup>*.

*<sup>i</sup>* , *j* = 3, 4 are computed in the same way as for the algorithm (RA1).

In case 2b) we have a unique matrix *A*. The fluxes on outgoing lines are computed as in the

In the first case we set (*γ*ˆ3, *γ*ˆ4) = *P*, while in the second case we set (*γ*ˆ3, *γ*ˆ4) = *Q*, where

Once solutions to RPs are given, one can use a Wave Front Tracking algorithm to construct a

The aim of this section is to verify that the assumptions underlying the data networks fluid-dynamic model (shortly FD model) are correct. Here we focus on the fixed-point models to describe TCP, and considering various set-ups with TCP traffic in a single bottleneck topology, we investigate queueing models for estimating packet loss rate. In what follows

It is reasonable to assume that the loss probability function *p* is null for some interval, which is a right neighborhood of zero. This means that at low densities no packet is lost. Then *p* should be increasing, reaching the value 1 at the maximal density, the situation of complete stuck. With the above assumptions the loss probability function in (4) can be written as:

2*ρ* −1

0, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> 1/2,

*<sup>ρ</sup>* , 1/2 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> 1. (16)

*n* ∑ *i*=1 *α i*,*s*,*d*,*j <sup>J</sup> <sup>π</sup>s*,*<sup>d</sup>* *γ*∈Ω*out*

*<sup>i</sup>* (*t*, *bi*−,*s*, *d*)*f*(*ρ*ˆ*i*)

*f*(*ρ*ˆ*j*)

F(*γ*,*rα*, G) where F is a convex

<sup>3</sup> .

functional which depends on *γ*, *r<sup>α</sup>* and on the set G of the routing standards.

*Q* = *proj*Ω*adm* (*P*). Again, we can extend to the case of *m* outgoing lines.

*<sup>i</sup>* , *j* = 3, 4 as in the case 2a):

*π*ˆ*j*(*t*, *aj*+,*s*, *d*) =

2 .

*p* (*ρ*) =

for every *t* ≥ 0, *j* ∈ {*n* + 1, ..., *n* + *m*},*s* ∈ S, *d* ∈ D.

sequence of approximate solutions.

we suppose *ρmax* = 1 and *σ* = <sup>1</sup>

**4.1 Loss probability function**

**4. Model assumptions**

a) G ∩ Ω*out* ∩ *r*<sup>Γ</sup> �= ∅,

The vector *π*ˆ *<sup>s</sup>*,*<sup>d</sup>*

**3.2.2 Case 2b)**

b) G ∩ Ω*out* ∩ *r*<sup>Γ</sup> = ∅ and *γ*3(*G*1) < *γ*3(*R*), c) G ∩ <sup>Ω</sup>*out* <sup>∩</sup> *<sup>r</sup>*<sup>Γ</sup> <sup>=</sup> <sup>∅</sup> and *<sup>γ</sup>*3(*G*1) <sup>&</sup>gt; *<sup>γ</sup>*max

case without sources and destinations.

We distinguish two cases:

a) *P* belongs to Ω, b) *P* is outside Ω.

Finally we define *π*ˆ *<sup>s</sup>*,*<sup>d</sup>*

Otherwise, if *r*<sup>Γ</sup> has a priority over G we set (*γ*ˆ3, *γ*ˆ4) = min

Let us consider the transmission of two consecutive routers. The node that transmits packets is called *sender*, while the receiving one is said *receiver*. Among the nodes, there is a link or channel, with limited capacity. Assume that the sender and the receiver are synchronized each other, i.e. the receiver is able to process in real time all packets, sent by the sender. In few words, no packets are lost. The packets loss can occur only on the link, due to its finite capacity. Under the zero buffer hypotheses the loss rate is defined as the proportional excess of offered traffic over the available capacity. If *R* is the sender bit rate and *C* is the link capacity, we have a loss if *R* > *C*. The model is said *proportional-excess* or briefly *P*/*E* and suppose deterministic arrivals. The packets bit rate is:

$$p = \begin{cases} 0, & \text{R} < \text{C}, \\ \frac{\text{R} - \text{C}}{\text{C}}, & \text{R} > \text{C}. \end{cases} \tag{17}$$

In Figure 3, loss probability for *P*/*E* model (continuous curve) and FD model (dashed curve) are shown, assuming *C* = *σ* = 1/2. For values *C* < *ρ* < 2*C*, the FD model overestimates the loss probability.

Fig. 3. Loss probabilities. Dashed line: FD model. Continuous line: P/E model.

Observe that the *P*/*E* model is not realistic. In fact, the sender and the receiver are never synchronized each other and whatever transmission protocol is used by the transport layer, the receiver has a finite length buffer, where the packets wait to be processed and eventually sent to the next node. Thus queueing models are needed, to infer about network performance.

#### **4.1.2 Models with finite capacity**

Queueing models are good at predicting loss in a network with many independent users, probably using different applications. Consider the traffic from TCP sources that send packets through a bottleneck link. The traffic is aggregated and used as an arrival process for the link. The arrival process, being the aggregation of independent sources, is approximated as a Poisson process, and the aggregated throughput is used as the rate of the Poisson process (see Wierman et al. (2003)). These considerations justify the assumption that the times between the packets arrivals are exponentially distributed. Depending on the hypothesis on the length of packets arriving to the queue the data transmission can be modelled with different queueing models, as *M*/*D*/1/*B* and *M*/*M*/1/*B*, characterized by deterministic and exponentially distributed lengths, respectively, and a buffer with capacity *B* − 1. From the queue length distribution, known in closed formulas or iteratively in the finite buffer case, expected time in queue and in the system, as well as packet loss rate can be derived. In what follows we denote the arrival intensity by *λ*, the service intensity by *μ* and define the load as *ρ* = *λ*/*μ*.

#### 4.1.2.1 Fixed packets dimension

In a scenario where all senders use the same data packets size, the queueing model *M*/*D*/1/*B* is the most natural choice. The probability that the buffer is full gives the loss rate:

$$p(\rho) = \frac{1 + (\rho - 1)\,\alpha\_B\left(\rho\right)}{1 + \rho\alpha\_B\left(\rho\right)},\tag{18}$$

*0.3 0.6* <sup>Σ</sup> *1.4 1.8* <sup>Ρ</sup>

*<sup>ρ</sup><sup>B</sup>* (<sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*) <sup>1</sup> <sup>−</sup> *<sup>ρ</sup>B*+<sup>1</sup> <sup>=</sup>

Dashed lines: *M*/*M*/1/*B*. Continuous line: P/E model.

events that are experienced by multiple users.

0.02 0.04 0.06 0.08 0.1 0.12 0.14 *<sup>p</sup>* <sup>Ρ</sup>

Fig. 6. Comparison of different queueing models.

lim *B*→∞ *0.02 0.04 0.06 0.08 0.10 <sup>p</sup>* <sup>Ρ</sup>

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 435

Fig. 5. Left: Loss bit rate for different values of the buffer. Right: Loss probability function.

In fact, the loss probability of the FD model represents for *σ* = 1 (up to a scale factor equal to

The loss probability for the queueing model (dashed line) and the *P*/*E* one (continuous line) is shown in Figure 5, right. The two curves almost match for small bit rate values, i.e. in the load range 0.9*σ* < *ρ* < 1.1*σ*. For greater loads values, the *P*/*E* model overestimates the loss

Theoretical and simulative studies pointed out that *M*/*D*/1/*B* and *M*/*M*/1/*B* queueing models give good prediction of the loss rate in network with many independent users performing short file transfers (shorts FTP). In literature other queueing models have been considered to describe different scenarios, as bach arrivals. For a comparison among different models see Figure 6, where the packet loss rate for *M*/*D*/1/*B*, *M*/*M*/1/*B*, *M*2/*M*/1/*B*, *M*5/*M*/1/*B* and the *P*/*E* models are reported for the case *B* = 100 and loads in the interval 0.8 < *ρ* < 1.1. Observe that *Mr*/*M*/1/*B* denotes a queue with Poisson batch arrivals of size *r* and describes the fact that TCP traffic is likely to be quite bursty due to synchronized loss

0.85 0.9 0.95 1 1.05 1.1

Significant difference are restricted to the range 0.9*σ* < *ρ* < 1.1*σ*. As the load increases above 1.1 the loss estimates become very close in the different queueing models. Any of these models

*<sup>M</sup><sup>5</sup><sup>M</sup><sup>1</sup><sup>B</sup> <sup>M</sup><sup>2</sup><sup>M</sup><sup>1</sup><sup>B</sup> <sup>M</sup><sup>D</sup><sup>1</sup><sup>B</sup> <sup>M</sup><sup>M</sup><sup>1</sup><sup>B</sup> <sup>P</sup><sup>E</sup>*

Ρ

0, 0 <sup>&</sup>lt; *<sup>ρ</sup>* <sup>≤</sup> 1,

*ρ*−1 *<sup>ρ</sup>* , *ρ* > 1.

*0.9* <sup>Σ</sup> *1.1* <sup>Ρ</sup>

*0.1 0.2 0.3 0.4 0.5 <sup>p</sup>* <sup>Ρ</sup>

2) a limit case of (19):

probability.

where

Fig. 4. Loss rates. Dashed line: *M*/*D*/1/*B* model. Continuous line: FD model.

Figure 4 shows a comparison among the loss rate (16) and (18), assuming *B* = 10. However, an *M*/*D*/1/*B* queue predicts a lower loss rate and higher throughput than is seen in the true network. This is due to fact that in real routers packet sizes are not always fixed to the maximum segment size, therefore packet sizes are more variable than a deterministic distribution.

4.1.2.2 Exponentially distributed packets size

Assume the packet size is exponentially distributed. This assumption is true if we consider the total amount of traffic as the superposition of traffic fluxes, coming from different TCP sources, each configured to use its own packet size. The *M*/*M*/1/*B* queue is a good approximation of the simulated bottleneck link shared among TCP sources under any traffic load (Wierman et al. (2003)). The loss rate for the *M*/*M*/1/*B* queueing model is:

$$p(\rho) = \frac{\rho^B \left(1 - \rho\right)}{1 - \rho^{B+1}}.\tag{19}$$

In Figure 5, left, the loss bit rate for different values of the buffer (*B* = 10, 20, 30) is reported. Notice that, increasing the *B* values, dashed lines tend to the continuous one.

14 Telecommunications Networks

packets arriving to the queue the data transmission can be modelled with different queueing models, as *M*/*D*/1/*B* and *M*/*M*/1/*B*, characterized by deterministic and exponentially distributed lengths, respectively, and a buffer with capacity *B* − 1. From the queue length distribution, known in closed formulas or iteratively in the finite buffer case, expected time in queue and in the system, as well as packet loss rate can be derived. In what follows we denote

In a scenario where all senders use the same data packets size, the queueing model *M*/*D*/1/*B*

*<sup>p</sup>*(*ρ*) = <sup>1</sup> <sup>+</sup> (*<sup>ρ</sup>* <sup>−</sup> <sup>1</sup>) *<sup>α</sup><sup>B</sup>* (*ρ*)

*<sup>k</sup>* (*<sup>B</sup>* <sup>−</sup> *<sup>k</sup>* <sup>−</sup> <sup>1</sup>)

*0.7 0.9 1.1* <sup>Ρ</sup>

Figure 4 shows a comparison among the loss rate (16) and (18), assuming *B* = 10. However, an *M*/*D*/1/*B* queue predicts a lower loss rate and higher throughput than is seen in the true network. This is due to fact that in real routers packet sizes are not always fixed to the maximum segment size, therefore packet sizes are more variable than a deterministic

Assume the packet size is exponentially distributed. This assumption is true if we consider the total amount of traffic as the superposition of traffic fluxes, coming from different TCP sources, each configured to use its own packet size. The *M*/*M*/1/*B* queue is a good approximation of the simulated bottleneck link shared among TCP sources under any traffic

*<sup>p</sup>*(*ρ*) = *<sup>ρ</sup><sup>B</sup>* (<sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*)

In Figure 5, left, the loss bit rate for different values of the buffer (*B* = 10, 20, 30) is reported.

Σ

<sup>1</sup> <sup>+</sup> *ρα<sup>B</sup>* (*ρ*) , (18)

<sup>1</sup> <sup>−</sup> *<sup>ρ</sup>B*+<sup>1</sup> . (19)

*<sup>k</sup> ρ<sup>k</sup> <sup>k</sup>*! , *<sup>B</sup>* <sup>≥</sup> 2.

the arrival intensity by *λ*, the service intensity by *μ* and define the load as *ρ* = *λ*/*μ*.

is the most natural choice. The probability that the buffer is full gives the loss rate:

*<sup>e</sup>ρ*(*B*−*k*−1) (−1)

4.1.2.1 Fixed packets dimension

*α<sup>B</sup>* (*ρ*) =

4.1.2.2 Exponentially distributed packets size

*B*−2 ∑ *k*=0

*0.05*

Fig. 4. Loss rates. Dashed line: *M*/*D*/1/*B* model. Continuous line: FD model.

load (Wierman et al. (2003)). The loss rate for the *M*/*M*/1/*B* queueing model is:

Notice that, increasing the *B* values, dashed lines tend to the continuous one.

*0.1*

*0.15*

*0.2 <sup>p</sup>* <sup>Ρ</sup>

where

distribution.

Fig. 5. Left: Loss bit rate for different values of the buffer. Right: Loss probability function. Dashed lines: *M*/*M*/1/*B*. Continuous line: P/E model.

In fact, the loss probability of the FD model represents for *σ* = 1 (up to a scale factor equal to 2) a limit case of (19):

$$\lim\_{B \to \infty} \frac{\rho^B \left(1 - \rho\right)}{1 - \rho^{B+1}} = \begin{cases} 0, & 0 < \rho \le 1, \\ \frac{\rho - 1}{\rho}, & \rho > 1. \end{cases}$$

The loss probability for the queueing model (dashed line) and the *P*/*E* one (continuous line) is shown in Figure 5, right. The two curves almost match for small bit rate values, i.e. in the load range 0.9*σ* < *ρ* < 1.1*σ*. For greater loads values, the *P*/*E* model overestimates the loss probability.

Theoretical and simulative studies pointed out that *M*/*D*/1/*B* and *M*/*M*/1/*B* queueing models give good prediction of the loss rate in network with many independent users performing short file transfers (shorts FTP). In literature other queueing models have been considered to describe different scenarios, as bach arrivals. For a comparison among different models see Figure 6, where the packet loss rate for *M*/*D*/1/*B*, *M*/*M*/1/*B*, *M*2/*M*/1/*B*, *M*5/*M*/1/*B* and the *P*/*E* models are reported for the case *B* = 100 and loads in the interval 0.8 < *ρ* < 1.1. Observe that *Mr*/*M*/1/*B* denotes a queue with Poisson batch arrivals of size *r* and describes the fact that TCP traffic is likely to be quite bursty due to synchronized loss events that are experienced by multiple users.

Fig. 6. Comparison of different queueing models.

Significant difference are restricted to the range 0.9*σ* < *ρ* < 1.1*σ*. As the load increases above 1.1 the loss estimates become very close in the different queueing models. Any of these models

see Figure 1. For the *P*/*E* model, we get

*0.1 0.2 0.3 0.4 0.5 0.6 <sup>f</sup>* <sup>Ρ</sup>

values smaller than the threshold.

consider two quantities:

define the following:

*<sup>f</sup>*(*ρ*) = *<sup>ρ</sup>v*¯, 0 *<sup>ρ</sup> <sup>σ</sup>*, (2*σ*−*ρ*)*v*¯*<sup>ρ</sup>*

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 437

<sup>Σ</sup> <sup>Ρ</sup>*max 0.4 0.8 1.2 1.6* <sup>Ρ</sup>

*0.2*

The flux in the *P*/*E* model and *M*/*M*/1/*B* model are depicted in Figure 8. Note the effects of a finite buffer on the maximal value of the flux. If *B* tends to infinity, the flux best approximates the FD model flux. For small *B* values, the maximal flux decreases and the load value in which the maximum is attained is shifted on the right due to the fact that packets are lost for load

We have a network (I, J ), with nodes of at most 2 × 2 type, and an initial data *ρ*<sup>0</sup> = (*ρi*,0)*i*=1,...,*N*. The evolution is determined by equation (9) on each line *Ii* and by Riemann Solvers *RSJ*, depending on priority and traffic distribution parameters, *q* and *α*, respectively. For the definition of *RSJ* see the case when the traffic distribution function is of type 2b).

We now consider *α* and *q* as controls. To measure the efficiency of the network, it is natural to

Clearly, to optimize 1) and 2) is the same if we refer to a single packet, but the averaged values may be very different (since there is a nonlinear relation among the two quantities). As the model consider macroscopic quantities, we can estimate the averages integrating over time and space the average velocity and the reciprocal of average velocity, respectively. We thus

> *Ii*

 *Ii* *v*(*ρi*(*t*, *x*)) *dx*,

1 *<sup>v</sup>*(*ρi*(*t*, *<sup>x</sup>*)) *dx*,

*0.4*

*0.6*

*0.8*

*1 <sup>f</sup>* <sup>Ρ</sup>

*0.2 0.4 0.6 0.8* <sup>Ρ</sup>

**5. Optimal control problems for telecommunication networks**

1) The average velocity at which packets travel through the network. 2) The average time taken by packets from source to destination.

> *J*1(*t*) = ∑ *i*

*J*2(*t*) = ∑ *i*

Now we state optimal control problems on the network.

Fig. 8. Flux. Left: P/E model. Right: *M*/*M*/1/*B* (for *B* = 5, *B* = 15, *B* = 25).

*<sup>σ</sup>* , *<sup>σ</sup> <sup>ρ</sup> <sup>ρ</sup>*max. (22)

Σ Ρ*max*

predict the loss rate equally well. However, under low loss environments, the best queueing model depends on the type of transfers by TCP sources, i.e. persistent or transient. It is shown in Olsen (2003) that *M*/*D*/1/*B* queues estimations of the loss rate can be used for transient sources. However, for sources with a slightly longer on and off periods, *M*/*M*/1/*B* queues best predict the loss rate, and for (homogeneous) persistent sources, *Mr*/*M*/1/*B* queues give better performance inferences, due to the traffic burstiness stemming from the TCP slow-start and source synchronization effect. Even if some models are more appropriate in situations of low load, others when the load is heavy, Figure 6 shows that the assumption on the loss probability function of the FD model is valid.

#### **4.2 Velocity**

The loss probability, influencing the average transmission time, has effects on the average velocity of packets:

$$v(\rho) = \bar{v} \left( 1 - p(\rho) \right) \dots$$

The behaviour of the average velocity in the FD model

$$v\left(\rho\right) = \begin{cases} \bar{v}, & 0 \le \rho \le /2, \\ \bar{v}\frac{1-\rho}{\rho}, 1/2 \le \rho \le 1, \end{cases} \tag{20}$$

is depicted in Figure 1. Notice that the velocity is constant if the system is free (no losses). Over the threshold, losses occur, and the average travelling time increasing reduces the velocity. The average packet velocity for the *P*/*E* model and the *M*/*M*/1/*B* model is plotted in Figure 7. Such two curves fit the curve of the FD model, confirming the goodness of its assumptions.

Fig. 7. Average velocity. Left: P/E model. Right: *M*/*M*/1/*B* model.

#### **4.3 Flux**

Once the velocity function is known, the flux is given by:

$$f(\rho) = v(\rho)\rho.$$

In case of the FD model

$$f\left(\rho\right) = \begin{cases} \overline{v}\rho, & 0 \le \rho \le 1/2, \\ \overline{v}(1-\rho), & 1/2 \le \rho \le 1, \end{cases} \tag{21}$$

see Figure 1. For the *P*/*E* model, we get

16 Telecommunications Networks

predict the loss rate equally well. However, under low loss environments, the best queueing model depends on the type of transfers by TCP sources, i.e. persistent or transient. It is shown in Olsen (2003) that *M*/*D*/1/*B* queues estimations of the loss rate can be used for transient sources. However, for sources with a slightly longer on and off periods, *M*/*M*/1/*B* queues best predict the loss rate, and for (homogeneous) persistent sources, *Mr*/*M*/1/*B* queues give better performance inferences, due to the traffic burstiness stemming from the TCP slow-start and source synchronization effect. Even if some models are more appropriate in situations of low load, others when the load is heavy, Figure 6 shows that the assumption on the loss

The loss probability, influencing the average transmission time, has effects on the average

*v*(*ρ*) = *v*¯ (1 − *p*(*ρ*)).

is depicted in Figure 1. Notice that the velocity is constant if the system is free (no losses). Over the threshold, losses occur, and the average travelling time increasing reduces the velocity. The average packet velocity for the *P*/*E* model and the *M*/*M*/1/*B* model is plotted in Figure 7. Such two curves fit the curve of the FD model, confirming the goodness of its assumptions.

Ρ

*f*(*ρ*) = *v*(*ρ*)*ρ*.

*<sup>v</sup>*¯*ρ*, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> 1/2,

*0.2 0.4 0.6 0.8 1 <sup>v</sup>* <sup>Ρ</sup>

Σ Ρ*max 0.4 1.2 1.6* <sup>Ρ</sup>

*v*¯ 1−*ρ*

*<sup>v</sup>*¯, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> /2,

*<sup>ρ</sup>* , 1/2 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> 1, (20)

Σ Ρ*max*

*<sup>v</sup>*¯(<sup>1</sup> <sup>−</sup> *<sup>ρ</sup>*), 1/2 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> 1, (21)

probability function of the FD model is valid.

The behaviour of the average velocity in the FD model

*v* (*ρ*) =

*0.2 0.4 0.6 0.8*

Once the velocity function is known, the flux is given by:

Fig. 7. Average velocity. Left: P/E model. Right: *M*/*M*/1/*B* model.

*f* (*ρ*) =

*0.2 0.4 0.6 0.8 1 <sup>v</sup>* <sup>Ρ</sup>

**4.3 Flux**

In case of the FD model

**4.2 Velocity**

velocity of packets:

$$f(\rho) = \begin{cases} \rho \overline{v}, & 0 \le \rho \le \sigma, \\ \frac{(2\sigma - \rho)\mathfrak{a}\rho}{\sigma}, \sigma \le \rho \le \rho\_{\text{max}}. \end{cases} \tag{22}$$

Fig. 8. Flux. Left: P/E model. Right: *M*/*M*/1/*B* (for *B* = 5, *B* = 15, *B* = 25).

The flux in the *P*/*E* model and *M*/*M*/1/*B* model are depicted in Figure 8. Note the effects of a finite buffer on the maximal value of the flux. If *B* tends to infinity, the flux best approximates the FD model flux. For small *B* values, the maximal flux decreases and the load value in which the maximum is attained is shifted on the right due to the fact that packets are lost for load values smaller than the threshold.

#### **5. Optimal control problems for telecommunication networks**

Now we state optimal control problems on the network.

We have a network (I, J ), with nodes of at most 2 × 2 type, and an initial data *ρ*<sup>0</sup> = (*ρi*,0)*i*=1,...,*N*. The evolution is determined by equation (9) on each line *Ii* and by Riemann Solvers *RSJ*, depending on priority and traffic distribution parameters, *q* and *α*, respectively. For the definition of *RSJ* see the case when the traffic distribution function is of type 2b).

We now consider *α* and *q* as controls. To measure the efficiency of the network, it is natural to consider two quantities:


Clearly, to optimize 1) and 2) is the same if we refer to a single packet, but the averaged values may be very different (since there is a nonlinear relation among the two quantities). As the model consider macroscopic quantities, we can estimate the averages integrating over time and space the average velocity and the reciprocal of average velocity, respectively. We thus define the following:

$$J\_1(t) = \sum\_{i} \int\_{I\_i} v(\rho\_i(t, \mathbf{x})) \, d\mathbf{x}\_i$$

$$J\_2(t) = \sum\_{i} \int\_{I\_i} \frac{1}{v(\rho\_i(t, \mathbf{x}))} \, d\mathbf{x}\_i$$

<sup>1</sup>−*<sup>ρ</sup><sup>ϕ</sup> <sup>ρ</sup><sup>ϕ</sup> <sup>H</sup> sϕ* 

with:

with

(8). Let

*s<sup>ϕ</sup>* =

*s<sup>ψ</sup>* =

Then, for *T* sufficiently big,

**(iv)** *if s*<sup>3</sup> = −*s*<sup>4</sup> = −1*, α* ∈

**(v)** *if s*<sup>3</sup> = −*s*<sup>4</sup> = +1*, α* ∈

, *ϕ* = 1, 2, *v*

<sup>−</sup>1, if *ρϕ*,0 <sup>≤</sup> <sup>1</sup>

+1 if *ρϕ*,0 > <sup>1</sup>

<sup>−</sup>1, if *ρψ*,0 <sup>&</sup>lt; <sup>1</sup>

<sup>+</sup>1 if *ρψ*,0 <sup>≥</sup> <sup>1</sup>

*q<sup>ϕ</sup>* =

 *ρψ* = *H*

−*s<sup>ψ</sup>* + <sup>1</sup>−*<sup>ρ</sup><sup>ψ</sup> ρψ H sψ* 

<sup>2</sup> and <sup>Γ</sup> <sup>=</sup> <sup>Γ</sup>*in*, or *ρϕ*,0 <sup>≤</sup> <sup>1</sup>

<sup>2</sup> and <sup>Γ</sup> <sup>=</sup> <sup>Γ</sup>*out*, or *ρψ*,0 <sup>≥</sup> <sup>1</sup>

*<sup>t</sup>* (*<sup>ρ</sup><sup>x</sup>*) <sup>=</sup> *<sup>ρ</sup><sup>x</sup>*

*<sup>β</sup>*<sup>−</sup> <sup>=</sup> <sup>Γ</sup> <sup>−</sup> *<sup>γ</sup>*max 3 *γ*max 3

*<sup>p</sup>*<sup>−</sup> <sup>=</sup> <sup>Γ</sup> <sup>−</sup> *<sup>γ</sup>*max 1 *γ*max 1

*J*2(*T*) *depend only on α. The optimal values for J*1(*T*) *are the following:*

**(ii)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *<sup>β</sup>*−*β*<sup>+</sup> <sup>=</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup>

 0, <sup>1</sup> 1+*β*<sup>+</sup> 

 <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1

*If* Γ = Γ*in, the optimal values for J*2(*T*) *are the following:* **(i)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +<sup>1</sup> *or sc* <sup>=</sup> <sup>−</sup>*sd* <sup>=</sup> <sup>−</sup>1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *<sup>α</sup>* <sup>=</sup> <sup>1</sup>

**(iv)** *if s*<sup>3</sup> <sup>=</sup> <sup>−</sup>*s*<sup>4</sup> <sup>=</sup> <sup>−</sup>1*, and* <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup>

**(ii)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup>

**(iii)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and* <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, <sup>α</sup>* <sup>∈</sup>

**(i)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *<sup>β</sup>*−*β*<sup>+</sup> <sup>&</sup>gt; 1, *or* <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, <sup>α</sup>* <sup>∈</sup>

**(iii)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *<sup>β</sup>*−*β*<sup>+</sup> <sup>&</sup>lt; <sup>1</sup>*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup>

<sup>1</sup> <sup>−</sup> *<sup>q</sup>*, if *<sup>ϕ</sup>* <sup>=</sup> 2, *αψ* <sup>=</sup>

*q*, if *ϕ* = 1,

<sup>2</sup> , *<sup>q</sup>ϕ*<sup>Γ</sup> <sup>&</sup>lt; *<sup>γ</sup>*max

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 439

<sup>2</sup> , *αψ*<sup>Γ</sup> <sup>&</sup>lt; *<sup>γ</sup>*max

function and *s<sup>ϕ</sup>* and *s<sup>ψ</sup>* are determined by the solution to the RP at *J*:

<sup>2</sup> , or *ρϕ*,0 <sup>≤</sup> <sup>1</sup>

<sup>2</sup> , or *ρϕ*,0 <sup>≥</sup> <sup>1</sup>

, *ψ* = 3, 4, where *H*(*x*) is the Heavyside

*<sup>ψ</sup>* and <sup>Γ</sup> <sup>=</sup> <sup>Γ</sup>*in*. *<sup>ψ</sup>* <sup>=</sup> 3, 4,

*<sup>ϕ</sup>* and Γ = Γ*out*,

 0, <sup>1</sup> 1+*β*<sup>+</sup> ;

*in the cases: <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, 1 <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> 1;

*in the cases: <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, 1 <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> 1.

2 ;

 0, <sup>1</sup> 1+*β*<sup>+</sup> ;

 <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 ;

*<sup>ϕ</sup>* and <sup>Γ</sup> <sup>=</sup> <sup>Γ</sup>*out*, *<sup>ϕ</sup>* <sup>=</sup> 1, 2,

<sup>2</sup> , *<sup>q</sup>ϕ*<sup>Γ</sup> <sup>=</sup> *<sup>γ</sup>*max

*<sup>ψ</sup>* and Γ = Γ*in*,

<sup>2</sup> , *αψ*<sup>Γ</sup> <sup>=</sup> *<sup>γ</sup>*max

*<sup>J</sup>*1(*T*) = <sup>2</sup> [*<sup>v</sup>* (*<sup>ρ</sup>*<sup>1</sup>) <sup>+</sup> *<sup>v</sup>* (*<sup>ρ</sup>*<sup>2</sup>) <sup>+</sup> *<sup>v</sup>* (*<sup>ρ</sup>*<sup>3</sup>) <sup>+</sup> *<sup>v</sup>* (*<sup>ρ</sup>*<sup>4</sup>)] ; (23)

*<sup>H</sup>* (*sx*) <sup>+</sup> *<sup>ρ</sup><sup>x</sup>* [*<sup>H</sup>* (−*sx*) <sup>−</sup> *<sup>H</sup>* (*sx*)] .

, *<sup>β</sup>*<sup>+</sup> <sup>=</sup> *<sup>γ</sup>*max 4 <sup>Γ</sup> <sup>−</sup> *<sup>γ</sup>*max 4 ,

, *<sup>p</sup>*<sup>+</sup> <sup>=</sup> *<sup>γ</sup>*max 2 <sup>Γ</sup> <sup>−</sup> *<sup>γ</sup>*max 2 .

> 0, <sup>1</sup> 1+*β*<sup>+</sup> ∪ <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 ;

We want to maximize the cost *J*1(*T*) and to minimize the cost *J*2(*T*) with respect to the parameters *α* and *q*. In Marigo (2006) and Cascone et al. (2007), you can find a similar approach for telecommunication networks and road networks, respectively, modelled with flux function

**Theorem 8.** *Consider a junction J of* 2 × 2 *type. If* Γ = Γ*in* = Γ*out and T is sufficiently big, the cost functionals J*1(*T*) *and J*2(*T*) *depend neither on α nor q. If* Γ = Γ*in, the cost functionals J*1(*T*) *and*

> 0, <sup>1</sup> 1+*β*<sup>+</sup> ;

 <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 ;

*<sup>J</sup>*2(*T*) = *<sup>t</sup>*(*<sup>ρ</sup>*<sup>1</sup>) <sup>+</sup> *<sup>t</sup>*(*<sup>ρ</sup>*<sup>2</sup>) <sup>+</sup> *<sup>t</sup>* (*<sup>ρ</sup>*<sup>3</sup>) <sup>+</sup> *<sup>t</sup>*(*<sup>ρ</sup>*<sup>4</sup>), (24)

 *α*, if *ψ* = 3, 1 − *α*, if *ψ* = 4.

and, to obtain finite values, we assume that the optimization horizon is given by [0, *T*] for some *T* > 0.

Notice that this corresponds to the following operation:


The value of such functionals depends on the order in which averages and integrations are taken.

Summarizing, we get the following optimal control problems:

**Data.** Network (I, J ); initial data *ρ*¯ = (*ρ*¯*i*)*i*=1,...,*N*; optimization horizon [0, *T*], *T* > 0.

**Dynamics.** Equation (9) on each line *I* ∈ I and Riemann Solver *RSJ* for each *J* ∈ J , depending on controls *α* and *q*.

**Control Variables.** Traffic distribution parameter *t* �→ *αJ*(*t*) and priority parameter *t* �→ *qJ*(*t*), i.e. two controls for every node *J* ∈ J .

**Control Space.** {(*αJ*, *qJ*) : *<sup>J</sup>* ∈ J , *<sup>α</sup>J*, *qJ* <sup>∈</sup> *<sup>L</sup>*∞([0, *<sup>T</sup>*], [0, 1])}.

**Cost functions.** Integrated functionals:

$$\max \int\_{0}^{T} J\_1(t) \, dt, \qquad \min \int\_{0}^{T} J\_2(t) \, dt.$$

**Definition 7.** *We call (Pi) the optimal control problem referred to the functional Ji:*

$$(\mathbf{P\_1}) \quad \max\_{(\mathbf{a}, \mathbf{q})} J\_{1\prime} \text{ subject to (9)}.$$

$$(\mathbf{P\_2}) \quad \min\_{(\mathbf{a}, \mathbf{q})} J\_{2\prime} \text{ subject to (9)}.$$

The direct solution of problems **(P***i***)** corresponds to a centralized approach. We propose the alternative approach of decentralized algorithm more precisely:

Step 1 For every node *J* and Riemann Solver *RSJ*, solve the simplified optimal control problem:

 $\max\left(or\min\right)$   $f\_l(T)\_\nu$ 

for *T* sufficiently big, on the network formed only by *J* with constant initial data, taking approximate solutions when there is lack of existence.


We consider a single node *J* with incoming lines, labelled by 1 and 2, and with outgoing lines, labelled by 3 and 4.

Since *<sup>ρ</sup>* <sup>=</sup> *<sup>γ</sup>*, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> <sup>1</sup> <sup>2</sup> , and *<sup>ρ</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>γ</sup>*, <sup>1</sup> <sup>2</sup> <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> 1, we have that *<sup>v</sup> ρϕ* = *H* −*s<sup>ϕ</sup>* +

<sup>1</sup>−*<sup>ρ</sup><sup>ϕ</sup> <sup>ρ</sup><sup>ϕ</sup> <sup>H</sup> sϕ* , *ϕ* = 1, 2, *v ρψ* = *H* −*s<sup>ψ</sup>* + <sup>1</sup>−*<sup>ρ</sup><sup>ψ</sup> ρψ H sψ* , *ψ* = 3, 4, where *H*(*x*) is the Heavyside function and *s<sup>ϕ</sup>* and *s<sup>ψ</sup>* are determined by the solution to the RP at *J*:

$$\begin{split} s\_{\boldsymbol{\Psi}} &= \begin{cases} -1, \text{ if } \rho\_{\boldsymbol{\Psi},0} \le \frac{1}{2} \text{ and } \Gamma = \Gamma\_{\text{in}}, \text{ or } \rho\_{\boldsymbol{\Psi},0} \le \frac{1}{2}, q\_{\boldsymbol{\Psi}}\Gamma = \gamma\_{\boldsymbol{\Psi}}^{\text{max}} \text{ and } \Gamma = \Gamma\_{\text{out}},\\ +1 \text{ if } \rho\_{\boldsymbol{\Psi},0} > \frac{1}{2}, \text{ or } \rho\_{\boldsymbol{\Psi},0} \le \frac{1}{2}, q\_{\boldsymbol{\Psi}}\Gamma < \gamma\_{\boldsymbol{\Psi}}^{\text{max}} \text{ and } \Gamma = \Gamma\_{\text{out}}, \end{cases} \\ s\_{\boldsymbol{\Psi}} &= \begin{cases} -1, \text{ if } \rho\_{\boldsymbol{\Psi},0} < \frac{1}{2}, \text{ or } \rho\_{\boldsymbol{\Psi},0} \ge \frac{1}{2}, a\_{\boldsymbol{\Psi}}\Gamma < \gamma\_{\boldsymbol{\Psi}}^{\text{max}} \text{ and } \Gamma = \Gamma\_{\text{in}},\\ +1 \text{ if } \rho\_{\boldsymbol{\Psi},0} \ge \frac{1}{2} \text{ and } \Gamma = \Gamma\_{\text{out}}, \text{ or } \rho\_{\boldsymbol{\Psi},0} \ge \frac{1}{2}, a\_{\boldsymbol{\Psi}}\Gamma = \gamma\_{\boldsymbol{\Psi}}^{\text{max}} \text{ and } \Gamma = \Gamma\_{\text{in}}. \end{cases} \end{split}$$

with:

18 Telecommunications Networks

and, to obtain finite values, we assume that the optimization horizon is given by [0, *T*] for

The value of such functionals depends on the order in which averages and integrations are

**Dynamics.** Equation (9) on each line *I* ∈ I and Riemann Solver *RSJ* for each *J* ∈ J ,

**Control Variables.** Traffic distribution parameter *t* �→ *αJ*(*t*) and priority parameter *t* �→ *qJ*(*t*),

*J*1(*t*) *dt*, min

The direct solution of problems **(P***i***)** corresponds to a centralized approach. We propose the

Step 1 For every node *J* and Riemann Solver *RSJ*, solve the simplified optimal control

Step 2 Apply the obtained optimal control at every time *t* in the optimization horizon and at

Notice that, for *T* sufficiently big, we can assume that the datum is constant on each line: this

We consider a single node *J* with incoming lines, labelled by 1 and 2, and with outgoing lines,

<sup>2</sup> <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> 1, we have that *<sup>v</sup>*

 *ρϕ* = *H*

−*s<sup>ϕ</sup>* +

max (*or* min) *Ji*(*T*), for *T* sufficiently big, on the network formed only by *J* with constant initial data, taking

**Definition 7.** *We call (Pi) the optimal control problem referred to the functional Ji:*

(**P1**) max (*α*,*q*)

(**P2**) min (*α*,*q*)  *T* 0

*J*1, *subject to (9).*

*J*2, *subject to (9).*

*J*2(*t*) *dt*.


**Data.** Network (I, J ); initial data *ρ*¯ = (*ρ*¯*i*)*i*=1,...,*N*; optimization horizon [0, *T*], *T* > 0.

Notice that this corresponds to the following operation:


Summarizing, we get the following optimal control problems:

**Control Space.** {(*αJ*, *qJ*) : *<sup>J</sup>* ∈ J , *<sup>α</sup>J*, *qJ* <sup>∈</sup> *<sup>L</sup>*∞([0, *<sup>T</sup>*], [0, 1])}.

max *T* 0

alternative approach of decentralized algorithm more precisely:

approximate solutions when there is lack of existence.

every node *J*, taking the value at *J* on each line as initial data.

<sup>2</sup> , and *<sup>ρ</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>γ</sup>*, <sup>1</sup>

some *T* > 0.

taken.

depending on controls *α* and *q*.

problem:

labelled by 3 and 4. Since *<sup>ρ</sup>* <sup>=</sup> *<sup>γ</sup>*, 0 <sup>≤</sup> *<sup>ρ</sup>* <sup>≤</sup> <sup>1</sup>

strongly simplifies the approach.

i.e. two controls for every node *J* ∈ J .

**Cost functions.** Integrated functionals:

$$\eta\_{\varphi} = \begin{cases} \eta\_{\prime} & \text{if } \varphi = 1, \\ 1 - \eta\_{\prime} \text{ if } \varphi = 2, \end{cases} \quad \mathfrak{a}\_{\psi} = \begin{cases} \mathfrak{a}\_{\prime} & \text{if } \psi = 3, \\ 1 - \mathfrak{a}\_{\prime} \text{ if } \psi = 4. \end{cases}$$

Then, for *T* sufficiently big,

$$J\_1(T) = 2\left[\boldsymbol{v}\left(\widehat{\rho}\_1\right) + \boldsymbol{v}\left(\widehat{\rho}\_2\right) + \boldsymbol{v}\left(\widehat{\rho}\_3\right) + \boldsymbol{v}\left(\widehat{\rho}\_4\right)\right];\tag{23}$$

$$J\_2(T) = t\left(\widehat{\rho}\_1\right) + t\left(\widehat{\rho}\_2\right) + t\left(\widehat{\rho}\_3\right) + t\left(\widehat{\rho}\_4\right) \,. \tag{24}$$

with

$$t\left(\widehat{\rho}\_{\mathbf{x}}\right) = \frac{\widehat{\rho}\_{\mathbf{x}}}{H\left(\mathbf{s}\_{\mathbf{x}}\right) + \widehat{\rho}\_{\mathbf{x}}\left[H\left(-\mathbf{s}\_{\mathbf{x}}\right) - H\left(\mathbf{s}\_{\mathbf{x}}\right)\right]}.$$

We want to maximize the cost *J*1(*T*) and to minimize the cost *J*2(*T*) with respect to the parameters *α* and *q*. In Marigo (2006) and Cascone et al. (2007), you can find a similar approach for telecommunication networks and road networks, respectively, modelled with flux function (8). Let

$$\begin{array}{c} \beta^{-} = \frac{\Gamma - \gamma\_{3}^{\max}}{\gamma\_{3}^{\max}}, \ \beta^{+} = \frac{\gamma\_{4}^{\max}}{\Gamma - \gamma\_{4}^{\max}},\\ p^{-} = \frac{\Gamma - \gamma\_{1}^{\max}}{\gamma\_{1}^{\max}}, \ p^{+} = \frac{\gamma\_{2}^{\max}}{\Gamma - \gamma\_{2}^{\max}}. \end{array}$$

**Theorem 8.** *Consider a junction J of* 2 × 2 *type. If* Γ = Γ*in* = Γ*out and T is sufficiently big, the cost functionals J*1(*T*) *and J*2(*T*) *depend neither on α nor q. If* Γ = Γ*in, the cost functionals J*1(*T*) *and J*2(*T*) *depend only on α. The optimal values for J*1(*T*) *are the following:*

**(i)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *<sup>β</sup>*−*β*<sup>+</sup> <sup>&</sup>gt; 1, *or* <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, <sup>α</sup>* <sup>∈</sup> 0, <sup>1</sup> 1+*β*<sup>+</sup> ; **(ii)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *<sup>β</sup>*−*β*<sup>+</sup> <sup>=</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup> 0, <sup>1</sup> 1+*β*<sup>+</sup> ∪ <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 ; **(iii)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *<sup>β</sup>*−*β*<sup>+</sup> <sup>&</sup>lt; <sup>1</sup>*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup> <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 ; **(iv)** *if s*<sup>3</sup> = −*s*<sup>4</sup> = −1*, α* ∈ 0, <sup>1</sup> 1+*β*<sup>+</sup> *in the cases: <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, 1 <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> 1; **(v)** *if s*<sup>3</sup> = −*s*<sup>4</sup> = +1*, α* ∈ <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 *in the cases: <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, 1 <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> 1. *If* Γ = Γ*in, the optimal values for J*2(*T*) *are the following:* **(i)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +<sup>1</sup> *or sc* <sup>=</sup> <sup>−</sup>*sd* <sup>=</sup> <sup>−</sup>1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *<sup>α</sup>* <sup>=</sup> <sup>1</sup> 2 ; **(ii)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup> 0, <sup>1</sup> 1+*β*<sup>+</sup> ; **(iii)** *if s*<sup>3</sup> <sup>=</sup> *<sup>s</sup>*<sup>4</sup> = +1*, and* <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, <sup>α</sup>* <sup>∈</sup> <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 ; **(iv)** *if s*<sup>3</sup> <sup>=</sup> <sup>−</sup>*s*<sup>4</sup> <sup>=</sup> <sup>−</sup>1*, and* <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup> 0, <sup>1</sup> 1+*β*<sup>+</sup> ;

Line *ρi*,0 *ρbi*,0 Line *ρi*,0 *ρbi*,0 Line *ρi*,0 *ρbi*,0 1 0.4 0.4 21 0.3 / 41 0.1 / 2 0.35 0.35 22 0.2 / 42 0.1 / 3 0.3 / 23 0.1 / 43 0.25 0 4 0.2 / 24 0.1 / 44 0.3 / 5 0.35 0.35 25 0.2 / 45 0.4 0.4 6 0.2 0 26 0.1 / 46 0.3 0.3 7 0.25 / 27 0.2 / 47 0.2 / 8 0.4 0.4 28 0.25 / 48 0.4 0 9 0.35 0.35 29 0.2 0 49 0.35 / 10 0.3 / 30 0.4 / 50 0.3 0 11 0.2 / 31 0.35 0.35 51 0.2 / 12 0.1 / 32 0.3 0.3 52 0.1 0 13 0.1 / 33 0.2 / 53 0.1 / 14 0.25 / 34 0.35 / 54 0.2 0 15 0.3 / 35 0.2 / 55 0.1 / 16 0.4 0.4 36 0.25 / 56 0.2 0 17 0.3 0 37 0.4 / 57 0.25 / 18 0.2 / 38 0.35 / 58 0.2 0 19 0.4 0.4 39 0.3 / 59 0.15 0 20 0.35 0.35 40 0.2 / 60 0.15 0

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 441

Table 1. Initial conditions and boundary data for the lines of the network for case A.

Table 2. Initial conditions and boundary data for the lines of the network for case C.

1. at each node parameters, that optimize the cost functionals *J*<sup>1</sup> and *J*<sup>2</sup> (*optimal case*);

of the simulation process.

simulation cases.

2. random *α* and *q* parameters (*static random case*) chosen in a random way at the beginning of the simulation process (for each simulation case, 100 static random simulations are made); 3. dynamic random parameters (*dynamic random case*) which change randomly at every step

In the following pictures, we show the values of the functionals *J*<sup>1</sup> and *J*2, computed on the whole network, as function of time. A legend for every picture indicates the different

The algorithm of optimization, which is of local type, can be applied to complex networks, without compromising the possibility of a global optimization. This situation is evident if we

Line *ρi*,0 *ρbi*,0 Line *ρi*,0 *ρbi*,0 Line *ρi*,0 *ρbi*,0 1 0.4 0.4 19 0.4 0.4 48 0.5 0.7 2 0.5 0.5 20 0.5 0.5 50 0.5 0.7 5 0.5 0.5 29 0.4 0.7 52 0.4 0.7 6 0.4 0.7 31 0.4 0.4 54 0.5 0.7 8 0.4 0.4 32 0.4 0.4 56 0.4 0.7 9 0.5 0.5 43 0.4 0.7 58 0.5 0.7 16 0.4 0.4 45 0.4 0.4 59 0.5 0.7 17 0.4 0.7 46 0.5 0.5 60 0.5 0.7

**(v)** *if s*<sup>3</sup> <sup>=</sup> <sup>−</sup>*s*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *or* <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup> <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 .

*If* Γ = Γ*out, the cost functionals J*1(*T*) *and J*2(*T*) *depend only on q. The optimal values for J*1(*T*) *and J*2(*T*) *are the same for α when* Γ = Γ*in, if we substitute α with q, β*<sup>−</sup> *with p*−*, and β*<sup>+</sup> *with p*+*.*

#### **5.1 A case study**

In what follows, we report the simulation results of a test telecommunication network, that consists of nodes of 2 × 2 type. The network, represented in Figure 9, is characterized by:


Fig. 9. Network with 24 nodes.

We distinguish three case studies, that can be called, case A, B, and C. In Table 1, we report the initial conditions *ρi*,0 and the boundary data (if necessary) *ρbi*,0 for case A.

As for case B, instead, we consider the same initial conditions of case A, but boundary data equal to 0.75.

Table 2 contains initial and boundary conditions for case C. An initial condition of 0.75 is assumed for the inner lines of the network, that are not present in Table 2.

As in Bretti et al. (2006), we consider approximations obtained by the numerical method of Godunov (Godunov (1959)), with space step Δ*x* = 0.0125 and time step determined by the CFL condition (Godlewsky et al. (1996)). The telecommunication network is simulated in a time interval [0, *T*], where *T* = 50 min. We study four simulation cases, choosing the flux function (7) or the flux function (8):

20 Telecommunications Networks

*If* Γ = Γ*out, the cost functionals J*1(*T*) *and J*2(*T*) *depend only on q. The optimal values for J*1(*T*) *and J*2(*T*) *are the same for α when* Γ = Γ*in, if we substitute α with q, β*<sup>−</sup> *with p*−*, and β*<sup>+</sup> *with p*+*.*

In what follows, we report the simulation results of a test telecommunication network, that consists of nodes of 2 × 2 type. The network, represented in Figure 9, is characterized by:

• 36 inner lines: 3, 4, 7, 10, 11, 12, 13, 14, 15, 18, 21, 22, 23, 24, 25, 26, 27, 28, 30, 33, 34, 35, 36, 37,

**15**

**16**

**19**

*17*

**20**

**22**

*43 44 57*

**23**

*58*

**24**

*59*

*60*

**21**

**17**

**18**

We distinguish three case studies, that can be called, case A, B, and C. In Table 1, we report

As for case B, instead, we consider the same initial conditions of case A, but boundary data

Table 2 contains initial and boundary conditions for case C. An initial condition of 0.75 is

As in Bretti et al. (2006), we consider approximations obtained by the numerical method of Godunov (Godunov (1959)), with space step Δ*x* = 0.0125 and time step determined by the CFL condition (Godlewsky et al. (1996)). The telecommunication network is simulated in a time interval [0, *T*], where *T* = 50 min. We study four simulation cases, choosing the flux

the initial conditions *ρi*,0 and the boundary data (if necessary) *ρbi*,0 for case A.

assumed for the inner lines of the network, that are not present in Table 2.

 <sup>1</sup> <sup>1</sup>+*β*<sup>−</sup> , 1 .

**(v)** *if s*<sup>3</sup> <sup>=</sup> <sup>−</sup>*s*<sup>4</sup> = +1*, and <sup>β</sup>*<sup>−</sup> <sup>≤</sup> <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*+, *or* <sup>1</sup> <sup>≤</sup> *<sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*+*, or <sup>β</sup>*<sup>−</sup> <sup>≤</sup> *<sup>β</sup>*<sup>+</sup> <sup>≤</sup> <sup>1</sup>*, <sup>α</sup>* <sup>∈</sup>

• 12 incoming lines: 1, 2, 5, 8, 9, 16, 19, 20, 31, 32, 45, 46; • 12 outgoing lines: 6, 17, 29, 43, 48, 50, 52, 54, 56, 58, 59, 60;

**1**

**7**

**6**

**11**

**12**

**13**

**14**

**8**

**9**

**10**

**2**

**3**

**4**

**5**

*46*

Fig. 9. Network with 24 nodes.

function (7) or the flux function (8):

equal to 0.75.

38, 39, 40, 41, 42, 44, 47, 49, 51, 53, 55, 57.

**5.1 A case study**

• 24 nodes;


Table 1. Initial conditions and boundary data for the lines of the network for case A.


Table 2. Initial conditions and boundary data for the lines of the network for case C.


In the following pictures, we show the values of the functionals *J*<sup>1</sup> and *J*2, computed on the whole network, as function of time. A legend for every picture indicates the different simulation cases.

The algorithm of optimization, which is of local type, can be applied to complex networks, without compromising the possibility of a global optimization. This situation is evident if we

As for the dynamic random simulation, its behaviour looks very similar to the optimal one for cases A and B (for case C, only *J*<sup>2</sup> presents optimal and dynamic random configurations, that are very similar). Hence, we could ask if it is possible to avoid the optimization of the network, and operate in dynamic random conditions. Indeed, this last case originates strange phenomena, that cannot be modelled, hence it is preferred to avoid such a situation for telecommunication network design. To give a confirmation of this intuition, focus the attention on line 13, that is completely inside the network and it is strongly influence by the dynamics at various nodes. In Figure 13, we see that, using optimal parameters, the density on line 13 shows a smoother profile than the one obtained through a dynamic random simulation.

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 443

*0.2 0.4 0.6 0.8 <sup>1</sup> <sup>x</sup>*

Fig. 13. Behaviour of the density on line 13 of the network of Figure 9, for *t* = 10, flux function (7), case C, in optimal and dynamic random simulations. Dashed line: optimal

A fluid-dynamic model for data networks has been described. The main advantages of this approach, with respect to existing ones, can be summarized as follows. The fluid-dynamic models are completely evolutive, thus they are able to describe the traffic situation of a network every instant of time, overcoming the difficulties encountered by many static models. An accurate description of queues formation and evolution on the network is possible. The theory permits the development of efficient numerical schemes for very large networks. The model is based on packets conservation at intermediate time scales, whose flux is determined via a loss probability function (at fast time scales) and on a semilinear equation for the evolution of the percentage of packets going from an assigned source to a given destination. The choice of the loss probability function is of paramount importance in order to achieve a feasible model. The fluid dynamic model has been compared with those obtained using various queueing paradigms, from proportional/excess to models with finite capacity, including different distributions for packet sizes. The final result is that such models give rise to velocity profiles and flux functions which are quite similar to the fluid dynamic ones. In order to solve dynamics at node,Riemann Solvers have been defined considering different traffic distribution functions (which indicate for each junction *J* the outgoing direction of traffic that started at source *s*, has *d* as final destination and reached *J* from an assigned incoming road) and rules RA1 and RA2. The algorithm RA1, already used for road traffic models, requires the definition of a traffic distribution matrix, whose coefficients describe the percentage of packets, forwarded from incoming lines to outgoing ones. Using the algorithm

*0.55 0.6 0.65 0.7 0.75* <sup>Ρ</sup>*10,x*

simulation for *J*2; solid line: dynamic random simulation.

**6. Conclusions**

Fig. 10. *J*<sup>1</sup> for flux function (8), case A, and zoom around the optimal and dynamic random case (right).

Fig. 11. *J*<sup>2</sup> for flux function (8), case B, and zoom around the optimal and dynamic random case (right).

Fig. 12. *J*<sup>1</sup> and *J*<sup>2</sup> for flux function (7), case C.

consider the behaviour of *J*<sup>1</sup> for case A and *J*<sup>2</sup> for case B. For cases A and B, the cost functionals simulated with flux function (7) are constant, which is not surprising since the initial data on the lines is less than <sup>1</sup> <sup>2</sup> . In case C, we present the behaviour of the cost functionals *J*<sup>1</sup> and *J*<sup>2</sup> for flux function (7). Boundary data are of Dirichlet type (unlike case A and B where we have considered Neumann boundary conditions) and the network is simulated with high incoming fluxes for the incoming lines and high initial conditions for inner lines. We can see, from Figure 12, that *J*<sup>1</sup> and *J*<sup>2</sup> are not constant as in cases A and B. Moreover, we have to take in mind that we have two different optimization algorithms for *J*<sup>1</sup> and *J*2. Notice that the dynamic random case follows the optimal case for *J*<sup>2</sup> and not for *J*1. Indeed, the optimal algorithm for *J*<sup>1</sup> presents an interesting aspect. When simulation begins, it is worst than the static random configuration. In the steady state, instead, the optimal configuration is the highest.

As for the dynamic random simulation, its behaviour looks very similar to the optimal one for cases A and B (for case C, only *J*<sup>2</sup> presents optimal and dynamic random configurations, that are very similar). Hence, we could ask if it is possible to avoid the optimization of the network, and operate in dynamic random conditions. Indeed, this last case originates strange phenomena, that cannot be modelled, hence it is preferred to avoid such a situation for telecommunication network design. To give a confirmation of this intuition, focus the attention on line 13, that is completely inside the network and it is strongly influence by the dynamics at various nodes. In Figure 13, we see that, using optimal parameters, the density on line 13 shows a smoother profile than the one obtained through a dynamic random simulation.

Fig. 13. Behaviour of the density on line 13 of the network of Figure 9, for *t* = 10, flux function (7), case C, in optimal and dynamic random simulations. Dashed line: optimal simulation for *J*2; solid line: dynamic random simulation.

#### **6. Conclusions**

22 Telecommunications Networks

*J2*

*38.225 38.25 38.275 38.3 38.325 38.35 38.375*

Fig. 10. *J*<sup>1</sup> for flux function (8), case A, and zoom around the optimal and dynamic random

*100.35 100.4 100.45 100.5 J2*

Fig. 11. *J*<sup>2</sup> for flux function (8), case B, and zoom around the optimal and dynamic random

consider the behaviour of *J*<sup>1</sup> for case A and *J*<sup>2</sup> for case B. For cases A and B, the cost functionals simulated with flux function (7) are constant, which is not surprising since the initial data on

for flux function (7). Boundary data are of Dirichlet type (unlike case A and B where we have considered Neumann boundary conditions) and the network is simulated with high incoming fluxes for the incoming lines and high initial conditions for inner lines. We can see, from Figure 12, that *J*<sup>1</sup> and *J*<sup>2</sup> are not constant as in cases A and B. Moreover, we have to take in mind that we have two different optimization algorithms for *J*<sup>1</sup> and *J*2. Notice that the dynamic random case follows the optimal case for *J*<sup>2</sup> and not for *J*1. Indeed, the optimal algorithm for *J*<sup>1</sup> presents an interesting aspect. When simulation begins, it is worst than the static random

configuration. In the steady state, instead, the optimal configuration is the highest.

<sup>2</sup> . In case C, we present the behaviour of the cost functionals *J*<sup>1</sup> and *J*<sup>2</sup>

*<sup>38</sup> <sup>40</sup> <sup>42</sup> <sup>44</sup> <sup>46</sup> <sup>48</sup> <sup>50</sup><sup>t</sup> min*

*12.82 12.84 12.86 12.88 12.9t min*

*<sup>5</sup> <sup>10</sup> <sup>15</sup> <sup>20</sup> <sup>25</sup> <sup>30</sup><sup>t</sup> min*

*static random dynamic random optimal*

*dynamic random optimal*

*dynamic random optimal*

*<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup><sup>t</sup> min*

*<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup><sup>t</sup> min*

*<sup>5</sup> <sup>10</sup> <sup>15</sup> <sup>20</sup> <sup>25</sup> <sup>30</sup><sup>t</sup> min*

Fig. 12. *J*<sup>1</sup> and *J*<sup>2</sup> for flux function (7), case C.

*static random dynamic random optimal*

*static random dynamic random optimal*

*static random dynamic random optimal*

case (right).

case (right).

the lines is less than <sup>1</sup>

A fluid-dynamic model for data networks has been described. The main advantages of this approach, with respect to existing ones, can be summarized as follows. The fluid-dynamic models are completely evolutive, thus they are able to describe the traffic situation of a network every instant of time, overcoming the difficulties encountered by many static models. An accurate description of queues formation and evolution on the network is possible. The theory permits the development of efficient numerical schemes for very large networks. The model is based on packets conservation at intermediate time scales, whose flux is determined via a loss probability function (at fast time scales) and on a semilinear equation for the evolution of the percentage of packets going from an assigned source to a given destination.

The choice of the loss probability function is of paramount importance in order to achieve a feasible model. The fluid dynamic model has been compared with those obtained using various queueing paradigms, from proportional/excess to models with finite capacity, including different distributions for packet sizes. The final result is that such models give rise to velocity profiles and flux functions which are quite similar to the fluid dynamic ones.

In order to solve dynamics at node,Riemann Solvers have been defined considering different traffic distribution functions (which indicate for each junction *J* the outgoing direction of traffic that started at source *s*, has *d* as final destination and reached *J* from an assigned incoming road) and rules RA1 and RA2. The algorithm RA1, already used for road traffic models, requires the definition of a traffic distribution matrix, whose coefficients describe the percentage of packets, forwarded from incoming lines to outgoing ones. Using the algorithm

Cascone, A.; D'Apice, C.; Piccoli, B. & Raritá, L. (2007). Optimization of traffic on road

Simulation and Optimal Routing of Data Flows Using a Fluid Dynamic Approach 445

Cascone, A.; Marigo, A.; Piccoli, B. & Rarità, L. (2010). Decentralized optimal routing for

Coclite, G.; Garavello, M. & and Piccoli, B. (2005). Traffic Flow on a Road Network, *SIAM Journal on Mathematical Analysis*, Vol. 36, 1862–1886, ISNN 0036-1410. Dafermos, C. (1999). *Hyperbolic Conservation Laws in Continuum Physics*, Springer-Verlag, ISBN

Daganzo, C. (1997). *Fundamentals of Transportation and Traffic Operations*, Pergamon-Elsevier,

D'Apice, C.; Manzo, R. & Piccoli, B. (2006). Packet flow on telecommunication networks, *SIAM Journal on Mathematical Analysis*, Vol. 38, No. 3, 717–740, ISNN 0036-1410. D'Apice, C.; Manzo, R. & Piccoli, B. (2008). A fluid dynamic model for telecommunication

D'Apice, C.; Manzo, R. & Piccoli, B. (2010). On the validity of fluid-dynamic models for data

Garavello, M. & Piccoli, B. (2006). *Traffic flow on networks*, AIMS Series on Applied

Godlewsky E. & Raviart, P. (1996). *Numerical Approximation of Hyperbolic Systems of Conservation Laws*, Springer Verlag, ISBN 978-0-387-94529-3, Heidelberg. Garavello, M. & Piccoli, B. (2005). Source-Destination Flow on a Road Network, *Communication in Mathematical Sciences*, Vol. 3, 261–283, ISSN 1539-6746. Godunov, S. K. (1959). A difference method for numerical calculation of discontinuous

Holden, H. & Risebro, N. H. (1995). A mathematical model of traffic flow on a network of

Marigo, A. (2006). Optimal distribution coefficients for telecommunication networks, *Networks*

Kelly, F.; Maulloo, A. K. & Tan, D. K. H. (1998). Rate control in communication networks:

Lighthill, M. J. & Whitham, G. B. (1955). On kinetic waves. II. Theory of Traffic Flows on Long

Newell, G. F. (1980). *Traffic Flow on Transportation Networks*, MIT Press, ISBN 0262140322,

Olsén, J. (2003). On Packet Loss Rates used for TCP Network Modeling, *Technical Report*,

Richards, P. I. (1956). Shock Waves on the Highway, *Oper. Res.*, Vol. 4, 42–51, ISSN 0030-364X. Serre, D. (1999). *Systems of conservation laws I and II*, Cambridge University Press, ISBN

networks with sources and destinations, *SIAM Journal on Applied Mathematics (SIAP)*,

Mathematics, vol. 1, American Institute of Mathematical Sciences, ISBN 1601330006,

solutions of the equations of hydrodynamics, *Mat. Sb.*, Vol. 47, 271–306, ISSN

unidirectional roads, *SIAM Journal on Mathematical Analysis*, Vol. 26, 999–1017, ISNN

shadow prices, proportional fairness and stability, *Journal of the Operational Research*

Crowded Roads, *Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.*, Vol. 229, 317–345, doi:

*(DCDS - B)*, Vol. 13, No. 1, 59–78, ISSN 15313492.

354064914X, New York.

United States.

0368-8666.

0036-1410.

ISBN 0080427855, Oxford.

Vol. 68, No. 4, 981–1003, ISSN 0036-1399.

networks, *Journal of Networks*, submitted, ISSN 1796-2056.

*and Heterogeneous Media*, Vol. 1, 315–336, ISNN 1556-1801.

*Society*, Vol. 49, 237–252, ISSN 0160-5682.

10.1098/rspa.1955.0089.

Cambridge (MA,USA).

521582334, 521633303, Cambridge.

Uppsala University.

networks, *Mathematical Models in Applied Sciences*, Vol. 17, 1587–1617, ISSN 0218-2025.

packets flow on data networks, *Discrete and Continuous Dynamical Systems - Series B*

RA2, not considered for urban traffic as redirections are not expected from modelling point of view (except in particular cases, as strong congestions or road closures), priority parameters, indicating priorities among flows of incoming lines, and distribution coefficients have to be assigned.

The main differences between the two algorithms are the following. The first one simply sends each packet to the outgoing line which is naturally chosen according to the final packet destination. The algorithm is blind to possible overloads of some outgoing lines and, by some abuse of notation, is similar to the behaviour of a "switch". The second algorithm, on the contrary, sends packets to outgoing lines in order to maximize the flux both on incoming and outgoing lines, thus taking into account the loads and possibly redirecting packets. Again by some abuse of notation, this is similar to a "router" behaviour. Hence, RA1 forwards packets on outgoing lines without considering the congestion phenomena, unlike RA2. Observe that a routing algorithm of RA1 type working through a routing table, according to which flows are sent with prefixed probabilities to the outgoing links, is of "distance vector" type. Reverse, an algorithm of RA2 type can redirect packets on the basis of link congestions, so it works on the link states (hence on their congestions) and so it is of "link-state" type.

The performance analysis of the networks was made through the use of different cost functionals, measuring average velocity and average travelling time, using the model consisting of the conservation law. The optimization is over parameters, which assign priority among incoming lines and traffic distribution among outgoing lines. A complete solution is provided in a simple case, and then used as local optimal choice for a complex test network. Three different choices of parameters have been considered: locally optimal, static random, and dynamic random (changing in time). The local optimal outperforms the others. Then, the behaviour of packets densities on the lines, that permits to rule out the dynamic random case has been analyzed.

All the optimization results have been obtained using a decentralized approach, i.e. an approach which sets local optimal parameters for each junction of the network. The cooperative aspect of such decentralized approach is the following. When a router optimizes the (local) functionals, it takes into considerations entering and exiting lines. Such lines reach other nodes, which benefit from the optimal choice. This in fact reflects in good global behavior as showed by simulations, described below. In future we aim to extend the optimization results to more general junctions and to explore global optimization techniques.

#### **7. References**


24 Telecommunications Networks

RA2, not considered for urban traffic as redirections are not expected from modelling point of view (except in particular cases, as strong congestions or road closures), priority parameters, indicating priorities among flows of incoming lines, and distribution coefficients have to be

The main differences between the two algorithms are the following. The first one simply sends each packet to the outgoing line which is naturally chosen according to the final packet destination. The algorithm is blind to possible overloads of some outgoing lines and, by some abuse of notation, is similar to the behaviour of a "switch". The second algorithm, on the contrary, sends packets to outgoing lines in order to maximize the flux both on incoming and outgoing lines, thus taking into account the loads and possibly redirecting packets. Again by some abuse of notation, this is similar to a "router" behaviour. Hence, RA1 forwards packets on outgoing lines without considering the congestion phenomena, unlike RA2. Observe that a routing algorithm of RA1 type working through a routing table, according to which flows are sent with prefixed probabilities to the outgoing links, is of "distance vector" type. Reverse, an algorithm of RA2 type can redirect packets on the basis of link congestions, so it works on

The performance analysis of the networks was made through the use of different cost functionals, measuring average velocity and average travelling time, using the model consisting of the conservation law. The optimization is over parameters, which assign priority among incoming lines and traffic distribution among outgoing lines. A complete solution is provided in a simple case, and then used as local optimal choice for a complex test network. Three different choices of parameters have been considered: locally optimal, static random, and dynamic random (changing in time). The local optimal outperforms the others. Then, the behaviour of packets densities on the lines, that permits to rule out the dynamic random case

All the optimization results have been obtained using a decentralized approach, i.e. an approach which sets local optimal parameters for each junction of the network. The cooperative aspect of such decentralized approach is the following. When a router optimizes the (local) functionals, it takes into considerations entering and exiting lines. Such lines reach other nodes, which benefit from the optimal choice. This in fact reflects in good global behavior as showed by simulations, described below. In future we aim to extend the optimization results to more general junctions and to explore global optimization techniques.

Alderson, D.; Chang, H.; Roughan, M.; Uhlig, S. & Willinger, W. (2007). The many facets

Baccelli, F.; Chaintreau, A.; De Vleeschauwer, D. & McDonald, D. (2006). HTTP turbulence,

Baccelli, F; Hong, D. & Liu, Z. (2001). Fixed points methods for the simulation of the sharing

Report RR-4154, INRIA, Le Chesnay Cedex, France), ISBN 0249-6399. Bretti, G.; Natalini, R. & Piccoli, B. (2006). Numerical approximations of a traffic flow model on networks, *Networks and Heterogeneous Media*, Vol. 1, 57–84, ISSN 1556-1801. Bressan, A. (2000). *Hyperbolic Systems of Conservation Laws - The One-dimensional Cauchy*

*Networks and Heterogeneous Media*, Vol. 1, 1–40, ISSN 1556-1801.

*Problem*, Oxford University Press, ISBN 0198507003, Oxford.

of internet topology and traffic, *Networks and Heterogenous Media*, Vol. 1, Issue 4,

of a local loop by large number of interacting TCP connections, *Proceedings of the ITC Specialist Conference on Local Loop*, 1–27, Barcelona, Spain, (also available in Technical

the link states (hence on their congestions) and so it is of "link-state" type.

assigned.

has been analyzed.

**7. References**

569–600, ISSN 1556-1801.


26 Telecommunications Networks

446 Telecommunications Networks – Current Status and Future Trends

Tanenbaum, A. S. (2003). *Computer Networks*, Prentice Hall, ISBN 0130661023, Upper Saddle

Wierman, A.; Osogami, T. & Olsén, J. (2003). A Unified Framework for Modeling TCP-Vegas,

Willinger, W. & Paxson, V. (1998). Where Mathematics meets the Internet, *Notices of the AMS*,

TCP-SACK, and TCP-Reno, *Proceedings of the IEEE/ACM International Symposium on modeling, Analysis and Simulation of Computer and Telecommunication Systems (MASCOTS)*, 269–278, ISBN 0-7695-2039-1, Orlando, Florida, October 2003, Los

River.

Alamitos, California, Washington.

Vol. 45, 961–970, ISSN 0002-9920.
