MODELLING TRANSIENT STATES IN QUEUEING MODELS OF COMPUTER NETWORKS: - - PowerPoint PPT Presentation

MODELLING TRANSIENT STATES IN QUEUEING MODELS OF COMPUTER NETWORKS: A FEW PRACTICAL ISSUES (engineering approach)

Tadeusz Czachórski

Institute of Theoretical and Applied Informatics

• f Polish Academy of Sciences,

IITiS PAN, Gliwice, Poland EuroNGI, EuroFGi, Engineering of Future Internet (of W. Burakowski)

ACROSS, Amsterdam, 13th – 14th February 2014

Why transient states analysis is needed ?

1 2 3 4 5 6 7 8 5000 10000 15000 20000 25000 30000 Bellcore

Real traffic intensity (Bellcore)

2

History

Origins: Agner Krarup Erlang, models to describe the Copenhagen telephone exchange.

• A. K. Erlang, The Theory of Probabilities and Telephone

Conversations, Nyt Tidsskrift Matematik, no. B 20, pp. 33-39, 1909,

• A. K. Erlang, Solutions of Some Problems in the Theory of

Probabilities of Significance in Automatic Telephone Exchnges, Electroteknikeren, no. 13, pp. 5-13, 1917.

• E. C. Molina, Application of the Theory of Probability to Telephone

Tranking Problems, Bell Syst. Tech. J., no. 6, pp. 461-494, 1927. Kolmogorow, Khinchin, Crommelin, Palm contributed

3

"customers" input flow

• utput flow

queue server

Known: − arrival pattern, e.g. interarrival time distribution − service time, e.g. service time distribution − queueing discipline − queue size limitations To determine: − queue distribution (or its moments) − waiting time distribution (or its moments) − losses (if limited queue)

Queueing model, generic problem

4

Kendall’s (David George Kendall) notation, the standard system used to describe and classify a queueing node:

A/B/c/N/H/D

• A denotes the type of interarrival times distribution,
• B – service time distribution, the same notattion as for A
• c – the number of servers
• N – the number of places in the system
• H – the calling population
• D – the queue’s discipline

5

Interpretation of customers and servers

• clients are people and service points are real-life locations such as

post office, counters, supermarket checkouts, corridors at a railway station to cross,

• clients might be ships docking at ports (here, "service time" is the

unloading time), parcels sent to warehouses (service time is the storage time), or vehicles arriving at junctions (the time taken to cross), water resevoirs (clients are cubic meters of incoming flows),

• in computer systems performed processes are clients and system

resources (central processing unit, external memories, CPU, transmission lines, terminals) are severs; service time is the time a process needs a resource,

• in telecommunication networks customers are packets (datagrams,

etc.) – block of data being sent to to a certain destination, the service time is the time needed to emit them at a node. The delays and losses introduced there influence the quality of service

6

input packets

• utput packets

queues (buffer overflow) lost packets

Router

7

Decade later: computer networks

Sender Receiver

• ther flows
• ther flows

router1 routerM

Transmission time = propagation time (fixed) + queueing times in routers (random, depending on current load) Quality of service = f ( transmission time, jitter(variability), loss probability )

8

2 4 6 8 10 12 14 10 20 30 40 50 60 70 80 90 100 Liczba pakietow/Jedn. czasu Czas [s], Jednostka czasu = 0.1 s

Self-similarity: Real traffic intensity, time unit = 0.1 s

9

10 20 30 40 50 60 100 200 300 400 500 600 700 800 900 1000 Liczba pakietow/Jedn. czasu Czas [s], Jednostka czasu = 1 s

Self-similarity: Real traffic intensity, time unit = 1 s

10

100 200 300 400 500 600 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Liczba pakietow/Jedn. czasu Czas [s], Jednostka czasu = 10 s

Self-similarity: Real traffic intensity, time unit = 10 s

11

1000 2000 3000 4000 5000 6000 7000 10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 Liczba pakietow/Jedn. czasu Czas [s], Jednostka czasu = 100 s

Self-similarity: Real traffic intensity, time unit = 100 s

12

Definition of self-similarity of the continuous time stochastic process:

A stochastic process X(t) is self-similar with Hurst parameter H, if for a positive factor g the process X(gt)/gH has the same distribution as the

• riginal process X(t).

In particular: E[X(t)] = E[X(gt)] gH V ar[X(t)] = V ar[X(gt)] g2H Rx(t, s) = Rx(gt, gs) g2H , where Rx(t, s) = E[X(t)X(s)]

13

Self-similarity of the discrete time stochastic process:

A sequence X1, X2, . . . Xk, . . . is self-similar, if for aggregated sequence X(m)

k

= 1 m(Xkm−m+1 + ... + Xkm), where m > 1 the variance is equal to V ar[X(m)] = V ar[X] mβ , gdzie 0 < β < 1, H = 1 − β/2 RX(m)(k) = RX(k).

14

Features of long-range dependent process (self-similar)

• for large m the variance of aggregated process X(m)

k

is equal to: V ar[X]m−β, 0 < β < 1,

• ∞

k=−∞ Cov(k) is divergent,

• spectrum S(ω) = ∞

−∞ R(k)e−jωk at ω = 0 is singular. 15

Features of short-range dependent process

• for large m the variance of aggregated process X(m)

k

is equal to:

V ar[X] m

• ∞

k=−∞ Cov(k) is convergent

• spectrum S(ω) = ∞

−∞ R(k)e−jωk at ω = 0 is finite. 16

Complexity of Internet, need of large topologies, see Caida www.caida.org/projects

17

• Transient queue analysis is needed to model time-dependent flows

and the dynamics of changes of router queues in computer networks.

• It is needed in stability analysis of control at Internet connections.
• It helps to predict the probability of packet loss and the queueing

delays which are the major factors of the quality of service.

• We need efficient modelling tools which are able to cope with large

network topologies typical for modern Internet.

18

Classical queueing theory: queues M/M/1, M/GI/1, GI/M/1, etc.

The simplest queue M/M/1 : exponential distributions A and B, fA(x) = λe−λx, fB(x) = µe−µx, one server Kolmogorov equations dp(0, t) dt = − p(0, t)λ + p(1, t)µ , dp(n, t) dt = − p(n, t)[λ + µ] + p(n + 1, t)µ + + p(n − 1, t)λ , n ≥ 1.

19

• r, in steady state, when limt→∞ p(n, t) = p(n), dp(n, t)/dt = 0

= − p(0)λ + pµ , = − p(n)[λ + µ] + p(n + 1)µ + p(n − 1)λ , n ≥ 0.

✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ❆ ❆ ❆ ❯ ❆ ❆ ❆ ❯ ❆ ❆ ❆ ❑ ❆ ❆ ❆ ❑ ❆ ❆ ❆ ❑ ❆ ❆ ❆ ❯ ❆ ❆ ❆ ❯ ❆ ❆ ❆ ❑ ♣ ♣ ♣ ✲ ✛ 1 2 K − 1 K λ µ µ µ µ µ λ λ λ λ

· · ·

20

Solution

p(n, t) = e(λ+µ)t ̺

n−i 2 In−i(at) + ̺ n−i−1 2

In+i+1(at)+ +(1 − ̺)̺n

∞

• j=n+i+2

̺

−j 2 Ij(at)

  , where a = 2µ√̺, and Ik(x) is modified Bessel function of first kind,

• rder k:

Ik(x) =

∞

• m=0

( x

2)k+2m

(k + m) ! m ! . initial condition p(n, 0) = δin, δin is Kronecker delta: t = 0 i customers is present in the system. practicaly the unique applicable result concerning transient states

21

in steady state p(n) = (1 − ̺)̺n , where ̺ = λ/µ < 1 in similar way M/M/c, M/M/c/N, M/M/c/N/H systems may be analyzed and closed form solution for p(n) is known. It is valid for open or closed networks of such stations – a product form solution, Jackson’s, Gordon-Newel’s models. We have some results for M/G/1 station (Pollaczek-Khinchin formula) and G/M/1 station based on embedded Markov chains but we have practically nothing applicable for G/G/1 stations.

22

Classical results do not fit well the problem of modelling IP routers, where

• the incoming streams are not Poisson and
• the size of packets is not exponentially distributed.
• the result p(n, t; n0) refers to transient states but it is assumed that

the model parameters, λ in particular, are constant.

23

We need models

• describing constantly changing non-Poisson flows and general

distributions of service times.

• allowing to include self similar flows.
• that are scalable to meet very large topologies characteristic to the

Internet. Below we summarize our experience with three approaches:

• Markovian queues solved numerically,
• fluid-flow approximation,
• diffusion approximation.

24

Markovian queues solved numerically

If the distributions A and B are not exponantial but they are a linear combination of exponential phases, we still may construct Kolmogorov equations, but they should be solved numerically.

• eg. if we use Cox distribution

✚✙ ✛✘ ✚✙ ✛✘ ✲ ✲ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❲ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❲ ✲ ✚✙ ✛✘ ✚✙ ✛✘ ✲ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❲ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❲ ✲ ✲ ✲ ✲ a1 a2 µ1e−µ1x µ2e−µ2x µk−1e−µk−1x µke−µkx 1 − a1 1 − a2 1 − ak−1 ak−1 ♣ ♣ ♣ 1 2 k−1 k

to represent A or B, the states of Markov chain are defined by n and phase of the corresponding distribution, and for C2 we will have state transition diagrams

25

in case of B

✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✲ ✲ ✲ 0, 0 1, 1 2, 1 3, 1 1, 2 2, 2 3, 2 ❅ ❅ ❅ ❅ ❅ ❅ ■ ❅ ❅ ❅ ❅ ❅ ❅ ■ ❅ ❅ ❅ ❅ ❅ ❅ ■ ❄ ❄ ❄ ✲ ✲ ✲ ✲ λ λ λ λ λ µ2 µ2 µ2 µ1a1 µ1a1 µ1a1 . . . . . . . . . ✓ ✓ ✴ µ1(1 − a1) µ1(1 − a1) µ1(1 − a1) ✛ ✚ ✘ ✙ ✡ ✡ ✢ ✡ ✡ ✢

and in case of A

✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✛ ✚ ✘ ✙ ✲ ✲ ✲ 1, 0 2, 0 1, 1 1, 2 ✲ λ1a1 λ2 λ1a1 . . . . . . ✛ ✚ ✘ ✙ 2, 1 λ2 ✛ ✚ ✘ ✙ ✲ ❏ ❏ ❪ ❏ ❏ ❪ ❙ ❙ ❙ ♦ µ µ µ ❙ ❙ ❙ ✇ λ1(1 − a1) ❙ ❙ ❙ ✇ λ1(1 − a1) 26

Cox distributions are convenient, because for ε > 0 there is such a Cox distribution widh fC(x), that for any practically used pdf fX(x) ∞

• fX(x) − fC(x)
• dx < ε

and

• ∞

x[fX(x) − fC(x)]dx

• < ε ,

but now recipe how to choose parameters, only numerically; some tools exist, see e.g. Reinecke, P., Krauß, T.,Wolter, K.: HyperStar: Phase-Type Fitting Made Easy. In: 9th International Conference on the Quantitative Evaluation of Systems (QEST) 2012. (September 2012) 201 ˝ U202 Tool Presentation.

27

Fitting histogram with C25

28

Problem: the number of states is growing very rapidly with the complexity of an object being modelled: each state of the Markov chain corresponds to one state of the system. The existing solvers as e.g. QNAP, XMARCA, PEPS, PEPSY, PRISM consider only steady state Markov chains. Theoretically, for any continuous time Markov chain the Chapman-Kolmogorov equations with transition matrix Q the solution is known dπ(t) dt = π(t)Q, (1) have the analytical solution: π(t) = π(0)eQt, (2) where π(t) is the probability vector and π(0) is the initial condition.

29

However, it is not easy to compute the expression eQt =

∞

• k=0

(Qt)k k! . (3) when Q is a large matrix, [C. Moler, Ch. Van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later!, SIAM REVIEW Vol. 45, No. 1, pp. 30-49]. The task is numerically unstable, especially for large Q. Additionally, to consier λ(t), we should make the parameters of the model piece-wise constant in small intervals and apply the solution at each of these intervals.

30

We are developping our own package OLYMP. It is a library for generating transition matrices of continuous time Markov chains (CTMC), and also for solving them. Definition Java language to define network nodes and the interactions between them. Generation of matrices is parallelized, and they are compressed on-the-fly using a dedicated compression based on finite-state automata. At the moment we are able with OLYMP to generate and solve Markov chains of the 150 million of states.

31

Implemented methods

• The basic method is Arnoldi – orthogonal projection process onto the

Krylov subspace: Km(P , v1) ≡ span{v1, P v1, ..., P m−1v1}, v1 - initial values. Approximation of the eigenvalue of P can be obtained from the eigenvalue of upper Hessenberg matrix: Hm = V T

mP V m,

where V m = [v1, v2, ..., vm].

• the state probability vector is approximated by βVmeHmτe1,

where β = v 2, v = Vme1, in vector ei the i-th element is equal 1 and the others are null.

32

We increase the size of tractable Markov chains by several orders through the use of a GPU-CPU (graphical processing unit) and a better design of computational algorithms for parallel computing and optimization of memory usage (uniformization method is implemented). GPU capabilities go far beyond the computer graphics. It is well known that a potential computational power of GPUs is much greater than that of contemporary CPUs (in a sense of the performance measured by number of floating point operations per second). Thus, it is possible to shorten the time of computations. Due to the enormous amount of the data to be processed, which does not fit in memory, methods must be developed to store vectors and matrices with intelligent management of memory.

33

Fluid-flow approximation

Fluid-flow approximation is a well known approach, e.g. [V. Misra, W.-B. Gong, D. Towsley: Fluid-based Analysis of a Network of AQM Routers Supporting TCP Flows with an Application to RED, ACM SIGCOMM 2000].

• f modelling transient behaviour where only mean values of traffic

intensity and service times are considered.

34

TCP Reno congestion window W size changes

• after receiving an ACK:

if W < W0 then W := W + 1 else W := W + 1/W

• after timeout: W := W/2

35

RED

min_th max_th K

µ λ

source min_th max_th K p_max 1

36

RED

p =                  if K ≤ minth

K−minth maxth−minth ∗ maxp

if minth < K ≤ maxth 1 if maxth < K moving average K K =                (1 − w) ∗ K + w ∗ n if queue is not empty (1 − w)

µ λ ∗ K

if queue is empty where w is a fixed (small) parameter and n is the instantaneous queue size (low pass filter). This is basic schema. Numerous variants

37

Dynamics of congestion window

˙ W(t) = 1 R(t) − W(t) 2 × W(t − R(t) R(t − R(t))[p(t − R(t)]

• ˙

W(t) = d W (t)

dt

the changes of the window size,

• R(t) round trip time RTT,
• W(t) window size at time t,
• p(t) probability of packet loss at time t,

38

bottleneck router queue

˙ q(t) =        −C + N(t) R(t) W(t), q > 0 max{0, −C + N(t) R(t) W(t)}, q = 0.

39

Unresponsive flows (UDP)

A way to incorporate in the model the unresponsive flows (UDP) noise is to see them as changes of C: Ceff = C − u0 ˙ q(t) =    −C + K(t)

R(t) W(t) + u(t),

q > 0 max[0, −(C − u0) + K(t)

R(t) W(t) + u(t)],

q = 0. (4)

40

Figure 1: Block diagram of differential equations

41

Linearization around working point (Wl0, p0, q0, u0) gives transfer functions Wl(s) = −Pwin(s)e−sR0p(s) l(s) = K R0 Wl(s) q(s) = Pque(s) [l(s) + u(s)] where Pwin(s) =

R0C2

eff

2K2

s +

2K R2

0Ceff

, Pque(s) = 1 s + Ceff

C 1 R0

, Ceff = C−u0.

42

RED queue

The process of taking the moving average (denoted below by x) may be modelled as in [Misra00] dx dt = loge(1 − w)/∆ − loge(1 − w)/∆q(t) (5) where ∆ is the interarrival time of packets and is taken as ∆ = 1/C. Hence, the transfer function of RED mechanism having changes of current queue δq at the entrance and changes of packet loss δp as the

• utput has the form

Caqm(s) = Lred k k + s where k = − ln(1−w)/∆ = −C ln(1−w) and Lred = maxp maxth − minth .

43

Nyquist criterion

G0(jω) – transfer function of the whole control loop.

44

Numerical example

C = 12500 packets/s, K = 60, RED parameters are pmax = 0.01, maxth = 200, minth = 100, w = 0.0001. Delay R is varied between 10 ms and 100ms, Ceff is equal 0.25C, 0.5C, and 0.75C.

45

Channel throughput C, changes in time

46

Congestion window, reaction on the changes of C

47

Queue length at router, reaction on the changes of C

48

Reno connections, G(jω) for various values of delay R, stable and unstable cases

49

Reno connections, G(jω) for various values of effective band Ceff, stable and unstable cases

50

Wired and wireless network, TCP-DCR

• S. Bhandarkar N. Sadry A.L.N. Reddy N. Vaidya: TCP-DCR: A novel

protocol for tolerating wireless channel errors, Technical Report TAMU-ECE-2003-01, February 2003. ˙ W(t) = 1 − PDα R(t) + PDα R(t) + rtt − W(t) 2 × ×W(t − R(t) − τ) R(t − R(t) − τ) [p(t − R(t) − τ) + PD]. where PD is congestion-independent loss probability in wireless part of the network, rtt is time after which the wireless protocol is able to recover from an error with probability α τ is an additional time the TCP protocol waits for the acknowledgement, giving thus a wireless protocol a chance to deal with the errors due to the imperfect transmission media.

51

• W
• W

q

• q

p

_

C

R 1 R 1 R 1

Time Delay R+ô

N 2 1 _ congestedqueue

1

wireless loss model Pd Queuing algorithm

(drop-tail,RED)

Queuing algorithm

(drop-tail,RED)

• rtt

R 1

Block diagram of differential equations in case of TCP-DCR

52

Linearization of f(x1, x2, . . . , xn) around workung point (x1,0, x2,0, . . . , xn,0) for small changes of arguments δf =

n

• i=1

∂f(x1, x2, . . . , xn) ∂xi

• (x1=x1,0,x2=x2,0,...,xn=xn,0)

δxi at working point f1 = f2 = 0, hence 1 − PDα R0 + PDα R0 + rtt − W 2 2R0 [p0 + PD] = 0 and −C + NW0 R0 = 0

53

linearization

˙ W(t) = f1

• W(t), W(t − R(t) − τ), q(t), q(t − R(t) − τ), p(t − R(t) − τ)
• =

f1(W, WR+τ, q, qR+τ, pR+τ) = 1 − PDα q/C + Tp + PDα q/C + Tp + rtt − W 2 WR+τ qR+τ/C + Tp + rtt[pR+τ + PD] ˙ q(t) = f2

• q(t), W(t)
• = f2(q, W)

δ ˙ W = ∂f1 ∂W δW + ∂f1 ∂WR+τ δWR+τ + ∂f1 ∂pR+τ δpR+τ + ∂f1 ∂q δq + ∂f1 ∂qR+τ δqR+τ δ ˙ q = ∂f2 ∂W δW + ∂f2 ∂q δq

54

where ∂f1 ∂W = −W0 2 p0 + PD R0 + rtt ∂f1 ∂WR+τ = ∂f1 ∂W ∂f1 ∂pR+τ = −W 2 2 1 R0 + rtt ∂f1 ∂q = −(1 − PDα)/C R2 − PDα/C (R0 + rtt)2 ∂f1 ∂qR+τ = W 2 2 [p0 + PD]/C (R0 + rtt)2 ∂f2 ∂W = N R0 ∂f2 ∂q = NW0 CR2

55

• 4
• 3
• 2
• 1

1 2 3 4 5

• 2
• 1

1 2 3 4

0.9 0.7 0.5 0.3 0.1

G(jω), the influence of the probability α = 0.1, 0.3, 0.5, 0.7, 0.9 of successful retransmission in the wireless part

56

−0.5 0.0 0.5 1.0 1.5 2.0 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

G(jω), the influence of the changes of rtt = 0.1RTT, . . . , 0.9RTT

57

• 0.5

0.0 0.5 1.0 1.5 2.0 2.5

• 1.5
• 1.0
• 0.5

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

CLT ALT

G(jω) The influence of multiplicative factor a to decrease the congestion window in case of a transmission (a = 1/2 in case

• f aggressive limited transmit NCR (ALT), and a = 2/3 in

case of careful limited transmit NCR (CLT)

58

The case of many connections and many routers

Network V be composed of routers. The instantaneous and average queue lengths in the network are noted by the vectors q and x. The network structure is represented by binary matrix A. Its rows correspond to TCP flows and the columns represent network nodes. If a flow i traverses a node k, the value of the element aik is determined as “one”, otherwise the element is set to “zero”. The matrix A and vector p(x) are used to define a new matrix B: Bij = Aij · pj(xj). The B matrix is used to calculate the total packet loss probability on the path of the entire flow.

59

The dynamics of the congestion window Wi(t) at a connection i is: dW i(t) dt = 1 Ri(q(t)) − Wi(t) 2 · Wi(t − τ) Ri(q(t − τ)) · ·  1 −

• j∈V

(1 − Bij)   (6) The rest is similar.

60

There are already solvers, e.g. [ Y. Sakumoto, R. Asai, H. Ohsaki, and M. Imase, "Design and implementation of flow-level simulator for performance evaluation of large scale networks", 15th Annual Meeting of the IEEE International Symposium on Modeling, Analysis, and Simulation of Computer and Telecommunication Systems (MASCOTS) 2007, Oct. 2007.] we made ours: Monika Nycz. The computation time depends heavily on the topology and the number of monitored nodes (all data stored).

61

Figure 2: Diagram of the tested topology. Test results: for a network of 150 000 nodes and 150 000 flows computation time in the range [2; 22] minutes, depending on the number of monitored nodes, the case of 22 minutes concerned monitoring of all 150 thousand nodes. Collecting the data from nine network key points took approximately 126 seconds. For a networks of 15 000 nodes and 15 000 flows, with monitoring nine key nodes : 11 seconds. Generally, it is also possible to monitor a network with over 1 million nodes, but the duration of the calculations is much longer. Therefore the current goal is to improve the imlementation efficiency.

62

Diffusion Approximation

Diffusion Approximation replaces the stochastic process N(t) – number

• f customers present in the system at time t – by a continuous diffusion

process X(t); ∂f(x, t; x0) ∂t = α 2 ∂2f(x, t; x0) ∂x2 − β ∂f(x, t; x0) ∂x defines the conditional pdf f(x, t; x0)dx = P[x ≤ X(t) < x + dx | X(0) = x0]. and f(n, t; n0) ≈ p(n, t; n0) α, β should reflect features of input stream, service times distribution and queueing discipline.

63

The both processes X(t) and N(t) have normally distributed changes; in case of a single server FIFO queue, the choice

• β = λ − µ,
• α = σ2

Aλ3 + σ2 Bµ3 = C2 Aλ + C2 Bµ

ensures the same ratio of time-growth of mean and variance of these distributions.

64

More formally,

• D. Iglehart, W. Whitt, Multiple Channel Queues in Heavy Traffic, Part

I-III, Advances in Applied Probability, vol. 2, pp. 150-177, 355-369, 1970: If ˆ Nn is a series of random variables derived from N(t): ˆ Nn = N(nt) − (λ − µ)nt (σ2

Aλ3 + σ2 Bµ3)√n ,

then the series is weakly convergent (in the sense of distribution) to standard Wiener process (if β > 0).

65

Boundary conditions

• reflecting barrier at x = 0,

[H. Kobayashi, Application of the diffusion approximation to queueing networks, Part 1: Equilibrium queue distributions, J.ACM, vol. 21, no. 2,

• pp. 316-328, Part 2: Nonequilibrium distributions and applications to

queueing modeling, J.ACM, vol. 21, no. 3, pp. 459-469, 1974] f(x, t; x0) = ∂ ∂x

• Φ

x − x0 − βt αt

• − e

2βx α Φ

x + x0 + βt αt

• heavy-load approximation

66

• barriers with instantaneous (elementary) jumps

[E. Gelenbe, On Approximate Computer Systems Models. Journal of ACM, 1975, vol. 22 1975, no. 2.] G/G/1/N station Barriers at x = 0, x = N. The model equations become ∂f(x, t; x0) ∂t = α 2 ∂2f(x, t; x0) ∂x2 − β ∂f(x, t; x0) ∂x + +λ0p0(t)δ(x − 1) + λNpN(t)δ(x − N + 1) , dp0(t) dt = lim

x→0 [α

2 ∂f(x, t; x0) ∂x − βf(x, t; x0)] − λ0p0(t) , dpN(t) dt = lim

x→N [−α

2 ∂f(x, t; x0) ∂x + βf(x, t; x0)] − λNpN(t) , where δ(x) is Dirac delta function.

67

Steady state f(x) =              λp0 −β (1 − ezx) for 0 < x ≤ 1 , λp0 −β (e−z − 1)ezx for 1 ≤ x ≤ N − 1 , µpN −β (ez(x−N) − 1) for N − 1 ≤ x < N , where z = 2β

α and p0, pN are determined through normalization

p0 = lim

t→∞ p0(t) = {1 + ̺ez(N−1) +

̺ 1 − ̺[1 − ez(N−1)]}−1 , pN = lim

t→∞ pN(t) = ̺p0ez(N−1) . 68

G/G/1/N Transient solution

• T. Czachórski, A method to solve diffusion equation with instantaneous

return processes acting as boundary conditions, Bulletin of Polish Academy of Sciences, Technical Sciences, vol. 41, no. 4, 1993. Diffusion on x ∈ (0, N) with 2 absorbing barriers – its pdf: φ(x, t; x0) =                δ(x − x0) 1 √ 2Π αt

∞

• n=−∞
• exp

βx′

n

α − (x − x0 − x′

n − βt)2

2αt

• − exp

βx′′

n

α − (x − x0 − x′′

n − βt)2

2αt

• where x′

n = 2nN, x′′ n = −2x0 − x′ n 69

f(x, t; ψ) is represented as a superposition of φ(x, t; ψ) f(x, t; ψ) = φ(x, t; ψ) + t g1(τ)φ(x, t − τ; 1)dτ + t gN−1(τ)φ(x, t − τ; N − 1)dτ

• g1(τ) – probability density that a new process φ(x, τ; 1) is started at

x = 1

• gN−1(τ) – probability density that a new process φ(x, τ; N − 1) is

started at x = N − 1

• g1(τ), gN−1(τ) obtained from the balance equations for the

probability mass entering barriers on the basis of first passage times 1 → 0, 1 → N and (N − 1) → 0, (N − 1) → N and sojourn times at barriers.

70

The system of equation is solved in terms of Laplace transforms, ¯ f(x, s; ψ) = ¯ φ(x, s; ψ) + ¯ g1(s) ¯ φ(x, s; 1) + ¯ gN−1(s) ¯ φ(x, s; N − 1) . then ¯ f(x, s; ψ) inversed numerically. Stehfest’s algorithm: a function f(t) is obtained from its transform ¯ f(s) for any fixed argument t as f(t) = ln 2 2

N

• i=1

Vi ¯ f ln 2 t i

• ,

where Vi = (−1)N/2+i

min(i,N/2)

• k=⌊ i+1

2 ⌋

kN/2+1(2k)! (N/2 − k)!k!(k − 1)!(i − k)!(2k − i)!. N is an even integer end depends on a computer precision.

71

• 0.1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 10 20 30 40 50 "ERLANG_plot/ERL_400m10.0" "ERL_400m10.0_N10_A512" "ERL_400m10.0_N30_A512" "ERL_400m10.0_N50_A512" "ERL_400m10.0_N70_A512" "ERL_400m10.0_N100_A512" "ERL_400m10.0_N120_A512" "ERL_400m10.0_N150_A512"

Quality of Stehfest’s algorithm as a function of N and the length of word, retransformation of E400

72

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 2 4 6 8 10 ’t = 0.1’ ’t = 0.5’ ’t = 1’ ’t = 2’

G/G/1/10 station, solution f(x, t; 5), (λ = 1, µ = 2)

73

0.1 0.2 0.3 0.4 0.5 0.6 2 4 6 8 10 ’t = 2’ ’t = 5’ ’t = 10’ ’t = 15’ ’t = 20’

G/G/1/10 station, solution f(x, t; 5), (λ = 1, µ = 2), transient state solution but for constant parameters ...

74

This solution is for transient states but for fixed values

• f α and β.

Otherwise we are obliged to divide the time interval into small subintervals where the parameters are kept constant. The solution at the end of one interval defines the initial condition for the next one. Hence, we may control the values of the parameters at each instant (e.g. they may be given by a self-similar input process)

75

Solution with time-dependent parameters, numerical example: G/G/1/30 queue

0.2 0.4 0.6 0.8 1 1.2 1.4 10 20 30 40 50 60 70 80 90 100 Input rate Lambda

Input traffic intensity λ(t).

76

1 2 3 4 5 6 7 8 10 20 30 40 50 60 70 80 90 100 Time in time-units G/G/1/N Model Mean Queue Length E[N] Diffusion Mean Queue Length E[N] Simulation

Mean number of customers (linear scale) as a function of time; diffusion approximation and simulation results.

77

0.0001 0.001 0.01 0.1 1 10 100 10 20 30 40 50 60 70 80 90 100 Time in time-units G/G/1/N Model Mean Queue Length E[N] Diffusion Logarithmic scale Mean Queue Length E[N] Simulation Logarithmic scale

Mean number of customers (logarithmic scale) as a function

• f time; diffusion approximation and simulation results.

78

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10 20 30 40 50 60 70 80 90 100 G/G/1/N Model p(0) Diffusion p(0) Simulation

Probability p(0, t) of the empty queue, diffusion approximation and simulation results.

79

5e-005 0.0001 0.00015 0.0002 0.00025 0.0003 10 20 30 40 50 60 70 80 90 100 Time in time-units G/G/1/N Model Mean Queue Length PN Diffusion Mean Queue Length PN Simulation

Probability p(N, t) of the saturated queue, diffusion approximation and simulation results.

80

2 3 9 13 12 8 7 14 16 10 5 17 18 19 15 6 11 1 4

1 2 3 4 5 6

Topology of the network considered in numerical examples

81

time units 0 − 20 20 − 25 25 − 65 65 − 70 70 − 75 75 − 100 node 8 0.8 1.7 1.7 1.7 1.7 1.3 node 9 1.8 1.8 1.5 1.0 1.0 1.0 node 19 1.0 1.0 2.5 2.5 2.5

• 1. 5

Generated flows λi, i = 8, 9, 19, as a function of time

82

0.5 1 1.5 2 2.5 3 10 20 30 40 50 60 70 80 90 100 Switch Mean Queue Node 9 Dif Mean Queue Node 9 Sim Mean Queue Node 10 Dif Mean Queue Node 10 Sim

Mean queue lengths at nodes 9 and 10, diffusion approximation and simulation results

83

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 10 20 30 40 50 60 70 80 90 100 Switch Mean Queue Node 17 Dif Mean Queue Node 17 Sim Mean Queue Node 18 Dif Mean Queue Node 18 Sim

Mean queue lengths at nodes 9 and 10, diffusion approximation and simulation results

84

0.5 1 1.5 2 2.5 3 3.5 4 10 20 30 40 50 60 70 80 90 100 Switch LaIn Node 8 Dif LaIn Node 8 Sim LaIn Node 9 Dif LaIn Node 9 Sim LaIn Node 19 Dif LaIn Node 19 Sim

Flow intensities at nodes 8, 9, and 19, diffusion approximation and simulation results

85

Self-similar, time varying input traffic

• As each small time interval we may introduce new parameters of the

diffusion process, we may give them values resulting from any autocorrelated or self-similar input proces.

• In numerical examples we took MMPP. If a modulator has M states,

and for he state i the input parameter is λi.

86

10 20 30 40 50 60 500 1000 1500 2000 2500 "case1-diffusion" "case1-simulation" "case2-diffusion" "case2-simulation"

Mean queue length, the source has two exponentially distributed phases of the same mean 1/α = 100; λ1 = 1.6, λ2 = 0.4 (average λ = 1.0, case 1) or λ2 = 0 (average λ = 0.8, case 2); constant service time; diffusion and simulation results.

87

10 20 30 40 50 60 70 80 90 100 200 400 600 800 1000 1200 1400 1600 1800 2000 "DIFF-init-state-25" "SIM-init-state-25" "DIFF-init-state-23" "SIM-init-state-23" "DIFF-init-state-21" "SIM-init-state-21"

Mean queue length, the source has 31 phases, initial phase j0 = 21, 23, 25; [ Maglaris, et al., packet video communications] constant service time; diffusion and simulation results.

88

2 4 6 8 10 12 100 200 300 400 500 600 700 800 900 1000 Mean queue length Sim EN Ro 0.75 Ca 1 Cb 1 Mean queue length Dif EN Ro 0.75 Ca 1 Cb 1 Mean queue length Sim EN Ro 0.75 Ca 1 Cb 3 Mean queue length Dif EN Ro 0.75 Ca 1 Cb 3 Mean queue length Sim EN Ro 0.75 Ca 1 Cb 5 Mean queue length Dif EN Ro 0.75 Ca 1 Cb 5 Mean queue length Sim EN Ro 0.75 Ca 1 Cb 8 Mean queue length Dif EN Ro 0.75 Ca 1 Cb 8

Mean queue length as a function of time, C2

A = 1, C2 B is changing, ̺ = 0.75. 89

• 0.2

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 100 200 300 400 500 600 700 800 900 1000 Mean queue length Absolute Error ErrA Ro 0.75 Ca 1 Cb 1 Mean queue length Absolute Error ErrA Ro 0.75 Ca 1 Cb 3 Mean queue length Absolute Error ErrA Ro 0.75 Ca 1 Cb 5 Mean queue length Absolute Error ErrA Ro 0.75 Ca 1 Cb 8

Absolute error, mean queue length, as a function of time, C2

A = 1, C2 B is changing 90

• 20

20 40 60 80 100 120 140 160 180 200 100 200 300 400 500 600 700 800 900 1000 Mean queue length Relative Error ErrR Ro 0.75 Ca 1 Cb 1 Mean queue length Relative Error ErrR Ro 0.75 Ca 1 Cb 3 Mean queue length Relative Error ErrR Ro 0.75 Ca 1 Cb 5 Mean queue length Relative Error ErrR Ro 0.75 Ca 1 Cb 8

Relative error, mean queue length, as a function of time, C2

A = 1, C2 B is changing 91

A comment to errors

Mean queue length, ̺ = 0.75, t = 75 time-units, see how it depends on C2

A, C2 B 92

Errors as a function of time and C2

A, C2 B

Absolute errors, mean queue length, ̺ = 0.75, t = 75 time-units

93

Relative errors [%], mean queue length, ̺ = 0.75, t = 75 time-units

94

Errors - conclusions

• Diffusion approximation error decreases with the growth of channel

utilization ̺.

• Absolute error (diffusion minus simulation) grows positively with the

growth of coefficient of variation C2

B.

• Absolute error grows negatively with the growth of coefficient of

variation C2

A.

• Diffusion error is relatively small in case where both coefficients

have similar values.

• Relative error does not exceed 20% which means method is well

suited for network modelling.

95

Here: state-dependent diffusion parameters, G/G/N/N station

If n − 1 < x < n, it is supposed that n channels are busy, hence we choose

• α(x, t) = λ(t)C2

A(t) + nµC2 B,

• β(x, t) = λ(t) − nµ

for n − 1 < x < n.

• r

α1, β1, α2, β2 . . . αN, βN For each subinterval a separate solution, fictive (absorbing) barriers between subintervals: probability flows are absorbed and reappear on the

• ther side at a small (e.g. 2−10) distance ε

96

γR

n (t)

n

x x

n + ε barrier barrier gn−ε(t) n n − ε γL

n (t)

gn+ε(t)

Flow balance for the barrier at x = n

γR

n (t) = gn−ε(t),

γL

n (t) = gn+ε(t),

n = 1, . . . , N − 1 with two exceptions concerning flows coming from barriers at x = 0: g1+ε(t) = γL

1 (t) + g1(t),

and x = N: gN−1−ε(t) = γR

N−1(t) + gN−1(t). 97

The density functions fi(x, t; ψi), i = 1, . . . N, for the intervals x ∈]i − 1, i[ are: f1(x, t; ψ1) = φ1(x, t; ψ1) + t g1−ε(τ)φ1(x, t − τ; 1 − ε)dτ , fn(x, t; ψn) = φn(x, t; ψn) + t gn−1+ε(τ)φn(x, t − τ; n − 1 + ε)dτ + + t gn−ε(τ)φn(x, t − τ; n − ε)dτ, n = 2, . . . N − 1, fN(x, t; ψN) = φN(x, t; ψN) + t gN−1+ε(τ)φN(x, t − τ; N − 1 + ε)dτ

98

The densities γR

n (t), γL n (t) are obtained in the similar way as in

G/G/1/N, γi,j(t) are the densities of first passage times between points (i, j): γ0(t) = p0(0)δ(t) + γψ1,0(t) + t g1−ε(τ)γ1−ε,0(t − τ)dτ , γL

1 (t)

= γψ1,1(t) + t g1−ε(τ)γ1−ε,1(t − τ)dτ , γL

n (t)

= γψn,n(t) + t gn−1+ε(τ)γn−1+ε,n(t − τ)dτ + + t gn−ε(τ)γn−ε,n(t − τ)dτ , n = 2, . . . N − 1 γR

n (t)

= γψn+1,n(t) + t gn+ε(τ)γn+ε,n(t − τ)dτ + + t gn+1−ε(τ)γn+1−ε,n(t − τ)dτ , n = 1, . . . N − 2 γR

N−1(t)

= γψN,N−1(t) + t gN−1+ε(τ)γN−1+ε,N−1(t − τ)dτ γN(t) = pN(0)δ(t) + γψN,N(t) + t gN−1+ε(τ)γN−1+ε,N(t − τ)dτ.

99

5 10 15 20 25 30 10 20 30 40 50 60 70 80 90 100 M/M/20/20 Model Intensity of output stream, Diffusion Intensity of output stream, Simulation Intensity of input stream

Input stream and output streams in case of exponential service - diffusion and simulation (100 000 repetitions), λout(t) = N

n=1 p(n, t)nµ

100

2 4 6 8 10 12 14 16 18 20 10 20 30 40 50 60 70 80 90 100 M/M/20/20 Model Mean Queue Length E[N], Diffusion Mean Queue Length E[N], Simulation

Mean number of customers as a function of time, diffusion approximation and simulation results

101

1 2 3 4 5 6 7 8 10 20 30 40 50 60 70 80 90 100 M/M/20/20 Model Intensity of lost customers, Diffusion Intensity of lost customers, Simulation

M/M/20/20, the flow of rejected customers, diffusion approximation and simulation results

102

0.05 0.1 0.15 0.2 0.25 0.3 10 20 30 40 50 60 70 80 90 100 M/M/20/20 Model Probability of saturated queue p(N,t), Diffusion Probability of saturated queue p(N,t), Simulation

M/M/20/20, probability p(20, t) that all channels are

• ccupied, diffusion approximation and simulation results

103

0.2 0.4 0.6 0.8 1 10 20 30 40 50 60 70 80 90 100 M/M/20/20 Model p(0,t), Diff. p(0,t), Sim.

M/M/20/20, probability p(0, t) that all channels are empty

104

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 5 10 15 20 M/M/20/20 Model PDF for 10 Diffusion PDF for 10 Simulation

M/M/20/20, distribution of the number of customers, t = 10 diffusion approximation and simulation results

105

0.02 0.04 0.06 0.08 0.1 0.12 5 10 15 20 M/M/20/20 Model PDF for 30 Diffusion PDF for 30 Simulation

M/M/20/20, distribution of the number of customers, t = 30 diffusion approximation and simulation results

106

0.02 0.04 0.06 0.08 0.1 0.12 0.14 5 10 15 20 M/M/20/20 Model PDF for 50 Diffusion PDF for 50 Simulation

M/M/20/20, distribution of the number of customers, t = 50 diffusion approximation and simulation results

107

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 5 10 15 20 M/M/20/20 Model PDF for 70 Diffusion PDF for 70 Simulation

M/M/20/20, distribution of the number of customers, t = 70 diffusion approximation and simulation results

108

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 10 20 30 40 50 60 70 80 90 100 M/M/2/2 Model E[N], Diffusion E[N], Simulation

M/M/2/2, mean number of customers as a function of time, diffusion approximation and simulation results

109

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10 20 30 40 50 60 70 80 90 100 M/M/2/2 Model Probability of empty queue P(0) Diffusion Probability of empty queue P(0) Simulation

M/M/2/2, p(0, t), diffusion approximation and simulation results

110

0.1 0.2 0.3 0.4 0.5 0.6 0.7 10 20 30 40 50 60 70 80 90 100 M/M/2/2 Model Probability of saturated queue P(N) Diffusion Probability of saturated queue P(N) Simulation

M/M/2/2, probability of saturated queue p(2, t), diffusion approximation and simulation results

111

a(x) d1(x) d2(x) a1(x) pBc a2(x) 1 − pBe pBe · · · 1 1 · · · Be a2(x) lower prority rejected cells 1 − pBc Bc a1(x) i i · · · · · · higher priority cells cells cells stream

Application to sliding window mechanism: during time T, Bc packets are allowed, then Bc best effort packets may enter,

• thers are rejected, correspnds to G/D/Bc/Bc model

112

5 10 15 20 25 30 10 20 30 40 50 60 70 80 90 100 M/D/20/20 Model Intensity of output stream, Diffusion Intensity of output stream, Simulation Intensity of input stream

Input stream and output streams in case of constant service, sliding window model – diffusion and simulation, λout(t) = λin(t − T)[(1 − p(N, t − T)]

113

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 5 10 15 20 M/D/20/20 Model PDF for 10 Diffusion PDF for 10 Simulation PDF for 30 Diffusion PDF for 30 Simulation

M/D/20/20, distribution of the number of customers, t = 10, 30 diffusion approximation and simulation results

114

0.2 0.4 0.6 0.8 1 5 10 15 20 M/D/20/20 Model PDF for 50 Diffusion PDF for 50 Simulation PDF for 70 Diffusion PDF for 70 Simulation

M/D/20/20, distribution of the number of customers, t = 50, 70 diffusion approximation and simulation results

115

λ

Electronic-optical access (edge) node

116

. . .

N −M +2 N −1

M 4

3 2 1

. . .

x N −M +1

1st interval barriers all packets are admitted to buffer all, except M-block, packets are admitted to buffer

• nly one-block packets are admitted to buffer

2nd 3rd Mth

N −M N −2

Diffusion subintervals

117

arrival of M-block packet, optical packet with N − 1 blocks is sent N-M+2 N-2 N-1

M 4

3 2 1

. . . . . .

x N-M N-M+1 arrival of 2-block packet, optical packet with N − 1 blocks is sent arrival of 4-block packet to empty buffer arrival of M-block packet to empty buffer arrival of 1-block packet, full optical packet is sent

Multiple jumps, but diffusion coefficients remain the same in all subintervals

118

Scalability

1 2 3 . . . . . . . . . 10 11 . . . 20 21 . . . 30 . . . . . . . . . 91 . . . 100 101 . . . 110 . . . 191 . . . 200 201 . . . 210 . . . 291 . . . 300 . . . . . . . . . . . . . . . 901 . . . 910 . . . 991 . . . 1000

1 9 1 9 1 9 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10 1 10

119

120

Conclusions

• Markovian models are flexible but have frequently enormous state

space and therefore are time and space consuming. Fluid flow approximation is the simplest, it may be easily applied to very large

• configurations. Diffusion approximation is somewhere between.
• each of them may be useful.
• there is also simulation: we have developed an extension of

OMNET++ of transient state models. (random generators were modified to make possible the changes of their parameters as a function of time, a new software was added to collect the statistics of multiple runs and to aggregate them)

• time consuming: a simulation run should be repeated sufficient

number of times (e.g. 500 thousands in our examples) and the results for a fixed time should be averaged.

• typically in some of our examples, 5 min of computation on a

standard PC station for a diffusion model is compared to 24 hours of simulations on the same machine.

121