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

in case of B µ 1 (1 − a 1 ) µ 1 (1 − a 1 ) µ 1 (1 − a 1 ) ✛ ✘ ✛ ✘ ✛ ✘ ✛ ✘ ✡ ✓ ✡ ✓ ✴ ✢ ✡ ✢ ✡ λ λ λ ✲ ✲ ✲ ✲ . . . . . . 0 , 0 1 , 1 2 , 1 3 , 1 ✚ ✙ ✚ ✙ ✚ ✙ ✚ ✙ ■ ❅ ■ ❅ ■ ❅ ❅ ❅ ❅ ❅ ❅ ❅ µ 1 a 1 µ 1 a 1 µ 1 a 1 ❅ ❅ ❅ µ 2 µ 2 µ 2 ✛ ✘ ✛ ✘ ❅ ✛ ✘ ❅ ❅ ❄ ❄ ❄ ❅ ❅ ❅ λ λ ✲ ✲ ✲ 2 , 2 3 , 2 1 , 2 . . . ✚ ✙ ✚ ✙ ✚ ✙ and in case of A λ 1 (1 − a 1 ) λ 1 (1 − a 1 ) ❙ ❙ ❙ ❙ ✛ ✘ ✛ ✘ ✛ ✘ ✛ ✘ ✛ ✘ ✇ ❙ ✇ ❙ λ 1 a 1 λ 2 λ 1 a 1 λ 2 ✲ ✲ ✲ ✲ ✲ . . . . . . 1 , 0 2 , 0 1 , 1 2 , 1 1 , 2 ✚ ✙ ✚ ✙ ✚ ✙ ✚ ✙ ✚ ✙ ❏ ❪ ❙ ♦ ❪ ❏ ❙ ❏ ❏ ❙ µ µ µ 26

Cox distributions are convenient, because for ε > 0 there is such a Cox distribution widh f C ( x ) , that for any practically used pdf f X ( x ) � ∞ � � � f X ( x ) − f C ( x ) � dx < ε � � 0 and � ∞ � � x [ f X ( x ) − f C ( x )] dx � < ε , � � � 0 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 C 25 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 ) = π ( t ) Q , (1) dt have the analytical solution: π ( t ) = π (0) e Q t , (2) where π ( t ) is the probability vector and π (0) is the initial condition. 29

However, it is not easy to compute the expression ∞ ( Q t ) k e Q t = � . (3) k ! k =0 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: K m ( P , v 1 ) ≡ span { v 1 , P v 1 , ..., P m − 1 v 1 } , v 1 - initial values. Approximation of the eigenvalue of P can be obtained from the eigenvalue of upper Hessenberg matrix: H m = V T where V m = [ v 1 , v 2 , ..., v m ] . m P V m , • the state probability vector is approximated by βV m e H m τ e 1 , where β = � v � 2 , v = V m e 1 , in vector e i 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]. of 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 < W 0 then W := W + 1 else W := W + 1 /W • after timeout: W := W/ 2 35

RED 1 K λ µ p_max source max_th min_th min_th max_th K 36

RED  0 if K ≤ min th       K − min th  max th − min th ∗ max p  p = if min th < K ≤ max th        1 if max th < K  moving average K  (1 − w ) ∗ K + w ∗ n      if queue is not empty   K = µ λ ∗ K (1 − w )       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 R ( t ) − W ( t ) 1 × W ( t − R ( t ) ˙ W ( t ) = R ( t − R ( t ))[ p ( t − R ( t )] 2 W ( t ) = d W ( t ) ˙ • the changes of the window size, dt • 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  − C + N ( t ) R ( t ) W ( t ) , q > 0    q ( t ) ˙ = 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 : C eff = C − u 0  − C + K ( t ) R ( t ) W ( t ) + u ( t ) , q > 0  q ( t ) ˙ = (4) max [0 , − ( C − u 0 ) + K ( t ) R ( t ) W ( t ) + u ( t )] , q = 0 .  40

Figure 1: Block diagram of differential equations 41

Linearization around working point ( W l 0 , p 0 , q 0 , u 0 ) gives transfer functions − P win ( s ) e − sR 0 p ( s ) W l ( s ) = K l ( s ) = W l ( s ) R 0 q ( s ) = P que ( s ) [ l ( s ) + u ( s )] where R 0 C 2 eff 1 2 K 2 P win ( s ) = , P que ( s ) = , C eff = C − u 0 . 2 K s + C eff 1 s + R 2 0 C eff C R 0 42

RED queue The process of taking the moving average (denoted below by x ) may be modelled as in [Misra00] dx dt = log e (1 − w ) / ∆ − log e (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 output has the form k C aqm ( s ) = L red k + s where max p k = − ln(1 − w ) / ∆ = − C ln(1 − w ) L red = . and max th − min th 43

Nyquist criterion G 0 ( jω ) – transfer function of the whole control loop. 44

Numerical example C = 12500 packets/s, K = 60 , RED parameters are p max = 0 . 01 , max th = 200 , min th = 100 , w = 0 . 0001 . Delay R is varied between 10 ms and 100ms, C eff is equal 0 . 25 C , 0 . 5 C , and 0 . 75 C . 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 C eff , 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. 1 − P D α R ( t ) + rtt − W ( t ) P D α ˙ W ( t ) = + × R ( t ) 2 × W ( t − R ( t ) − τ ) R ( t − R ( t ) − τ ) [ p ( t − R ( t ) − τ ) + P D ] . where P D 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

congested�queue 1 1 R C � q W W q � 1 _ 1 � � N R � rtt R _ � 1 2 Queuing algorithm Queuing algorithm � � Time 1 Delay (drop-tail,�RED) (drop-tail,�RED) R R+ô p wireless loss � P d model Block diagram of differential equations in case of TCP-DCR 52

Linearization of f ( x 1 , x 2 , . . . , x n ) around workung point ( x 1 , 0 , x 2 , 0 , . . . , x n, 0 ) for small changes of arguments n � ∂f ( x 1 , x 2 , . . . , x n ) � � δf = δx i � ∂x i � ( x 1 = x 1 , 0 ,x 2 = x 2 , 0 ,...,x n = x n, 0 ) i =1 at working point f 1 = f 2 = 0 , hence R 0 + rtt − W 2 1 − P D α P D α 0 + [ p 0 + P D ] = 0 R 0 2 R 0 and − C + NW 0 = 0 R 0 53

linearization � � ˙ W ( t ) = f 1 W ( t ) , W ( t − R ( t ) − τ ) , q ( t ) , q ( t − R ( t ) − τ ) , p ( t − R ( t ) − τ ) = f 1 ( W, W R + τ , q, q R + τ , p R + τ ) = 1 − P D α q/C + T p + rtt − W P D α W R + τ q/C + T p + q R + τ /C + T p + rtt [ p R + τ + P D ] 2 � � q ( t ) ˙ = f 2 q ( t ) , W ( t ) = f 2 ( q, W ) = ∂f 1 ∂f 1 ∂p R + τ δp R + τ + ∂f 1 ∂f 1 ∂f 1 δ ˙ W ∂W δW + ∂W R + τ δW R + τ + ∂q δq + ∂q R + τ δq R + τ ∂W δW + ∂f 2 ∂f 2 δ ˙ q = ∂q δq 54

where ∂f 1 − W 0 p 0 + P D = ∂W 2 R 0 + rtt ∂f 1 ∂f 1 = ∂W R + τ ∂W − W 2 ∂f 1 1 0 = ∂p R + τ 2 R 0 + rtt ∂f 1 − (1 − P D α ) /C P D α/C = − R 2 ( R 0 + rtt ) 2 ∂q 0 W 2 ∂f 1 [ p 0 + P D ] /C 0 = ( R 0 + rtt ) 2 ∂q R + τ 2 ∂f 2 N = ∂W R 0 ∂f 2 NW 0 = CR 2 ∂q 0 55

4 0.5 0.7 3 0.3 0.1 2 1 0.9 0 -4 -3 -2 -1 0 1 2 3 4 5 -1 -2 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

3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 −0.5 0.0 0.5 1.0 1.5 2.0 −0.5 −1.0 G ( jω ) , the influence of the changes of rtt = 0 . 1 RTT, . . . , 0 . 9 RTT 57

3.5 3.0 2.5 ALT 2.0 1.5 1.0 0.5 0.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -0.5 -1.0 CLT -1.5 G ( jω ) The influence of multiplicative factor a to decrease the congestion window in case of a transmission ( a = 1 / 2 in case of 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 a ik 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 : B ij = A ij · p j ( x j ) . 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 W i ( t ) at a connection i is: dW i ( t ) R i ( q ( t )) − W i ( t ) 1 W i ( t − τ ) = · R i ( q ( t − τ )) · 2 dt   � ·  1 − (1 − B ij ) (6)  j ∈ V 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 of customers present in the system at time t – by a continuous diffusion process X ( t ) ; ∂ 2 f ( x, t ; x 0 ) ∂f ( x, t ; x 0 ) = α − β ∂f ( x, t ; x 0 ) ∂x 2 ∂t 2 ∂x defines the conditional pdf f ( x, t ; x 0 ) dx = P [ x ≤ X ( t ) < x + dx | X (0) = x 0 ] . and f ( n, t ; n 0 ) ≈ p ( n, t ; n 0 ) α , β 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 • β = λ − µ , A λ 3 + σ 2 B µ 3 = C 2 • α = σ 2 A λ + C 2 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 ˆ N n is a series of random variables derived from N ( t ) : N n = N ( nt ) − ( λ − µ ) nt ˆ B µ 3 ) √ n , A λ 3 + σ 2 ( σ 2 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] � � x − x 0 − βt � � x + x 0 + βt �� f ( x, t ; x 0 ) = ∂ 2 βx α Φ − e Φ ∂x α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 ∂ 2 f ( x, t ; x 0 ) ∂f ( x, t ; x 0 ) α − β ∂f ( x, t ; x 0 ) = + ∂x 2 ∂t 2 ∂x + λ 0 p 0 ( t ) δ ( x − 1) + λ N p N ( t ) δ ( x − N + 1) , dp 0 ( t ) x → 0 [ α ∂f ( x, t ; x 0 ) = lim − βf ( x, t ; x 0 )] − λ 0 p 0 ( t ) , dt 2 ∂x dp N ( t ) x → N [ − α ∂f ( x, t ; x 0 ) = lim + βf ( x, t ; x 0 )] − λ N p N ( t ) , dt 2 ∂x where δ ( x ) is Dirac delta function. 67

Steady state  λp 0 − β (1 − e zx ) 0 < x ≤ 1 , for      λp 0  − β ( e − z − 1) e zx f ( x ) = 1 ≤ x ≤ N − 1 , for  µp N  − β ( e z ( x − N ) − 1)  N − 1 ≤ x < N , for    where z = 2 β α and p 0 , p N are determined through normalization ̺ t →∞ p 0 ( t ) = { 1 + ̺e z ( N − 1) + 1 − ̺ [1 − e z ( N − 1) ] } − 1 , p 0 = lim t →∞ p N ( t ) = ̺p 0 e z ( N − 1) . p N = lim 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 − x 0 )    ∞ n − βt ) 2 � βx ′ − ( x − x 0 − x ′ �  1 �  � n  √ exp  φ ( x, t ; x 0 ) = α 2 αt 2 Π αt n = −∞  n − βt ) 2 � βx ′′ �� − ( x − x 0 − x ′′   n − exp    α 2 αt  where x ′ n = 2 nN , x ′′ n = − 2 x 0 − x ′ n 69

f ( x, t ; ψ ) is represented as a superposition of φ ( x, t ; ψ ) � t f ( x, t ; ψ ) = φ ( x, t ; ψ ) + g 1 ( τ ) φ ( x, t − τ ; 1) dτ 0 � t + g N − 1 ( τ ) φ ( x, t − τ ; N − 1) dτ 0 • g 1 ( τ ) – probability density that a new process φ ( x, τ ; 1) is started at x = 1 • g N − 1 ( τ ) – probability density that a new process φ ( x, τ ; N − 1) is started at x = N − 1 • g 1 ( τ ) , g N − 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 ; ψ ) = ¯ ¯ g 1 ( s ) ¯ g N − 1 ( s ) ¯ φ ( x, s ; ψ ) + ¯ φ ( x, s ; 1) + ¯ φ ( 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 N � ln 2 � f ( t ) = ln 2 V i ¯ � f t i , 2 i =1 where min( i,N/ 2) k N/ 2+1 (2 k )! � V i = ( − 1) N/ 2+ i ( N/ 2 − k )! k !( k − 1)!( i − k )!(2 k − i )! . k = ⌊ i +1 2 ⌋ N is an even integer end depends on a computer precision. 71

0.8 "ERLANG_plot/ERL_400m10.0" "ERL_400m10.0_N10_A512" "ERL_400m10.0_N30_A512" 0.7 "ERL_400m10.0_N50_A512" "ERL_400m10.0_N70_A512" "ERL_400m10.0_N100_A512" "ERL_400m10.0_N120_A512" 0.6 "ERL_400m10.0_N150_A512" 0.5 0.4 0.3 0.2 0.1 0 -0.1 0 10 20 30 40 50 Quality of Stehfest’s algorithm as a function of N and the length of word, retransformation of E 400 72

0.8 ’t = 0.1’ ’t = 0.5’ ’t = 1’ 0.7 ’t = 2’ 0.6 0.5 0.4 0.3 0.2 0.1 0 0 2 4 6 8 10 G/G/1/10 station, solution f ( x, t ; 5) , ( λ = 1 , µ = 2 ) 73

0.6 ’t = 2’ ’t = 5’ ’t = 10’ ’t = 15’ ’t = 20’ 0.5 0.4 0.3 0.2 0.1 0 0 2 4 6 8 10 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 of α 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 Input rate 1.4 Lambda 1.2 1 0.8 0.6 0.4 0.2 0 0 10 20 30 40 50 60 70 80 90 100 Input traffic intensity λ ( t ) . 76

G/G/1/N Model 8 Mean Queue Length E[N] Diffusion Mean Queue Length E[N] Simulation 7 6 5 4 3 2 1 0 0 10 20 30 40 50 60 70 80 90 100 Time in time-units Mean number of customers (linear scale) as a function of time; diffusion approximation and simulation results. 77

G/G/1/N Model 100 Mean Queue Length E[N] Diffusion Logarithmic scale Mean Queue Length E[N] Simulation Logarithmic scale 10 1 0.1 0.01 0.001 0.0001 0 10 20 30 40 50 60 70 80 90 100 Time in time-units Mean number of customers (logarithmic scale) as a function of time; diffusion approximation and simulation results. 78

G/G/1/N Model 1 p(0) Diffusion p(0) Simulation 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 10 20 30 40 50 60 70 80 90 100 Probability p (0 , t ) of the empty queue, diffusion approximation and simulation results. 79

G/G/1/N Model Mean Queue Length PN Diffusion Mean Queue Length PN Simulation 0.0003 0.00025 0.0002 0.00015 0.0001 5e-005 0 0 10 20 30 40 50 60 70 80 90 100 Time in time-units Probability p ( N, t ) of the saturated queue, diffusion approximation and simulation results. 80

2 1 3 5 4 1 7 8 12 9 13 10 6 11 2 14 16 15 3 17 4 18 5 6 19 Topology of the network considered in numerical examples 81

0 − 20 20 − 25 25 − 65 65 − 70 70 − 75 75 − 100 time units 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

Switch 3 Mean Queue Node 9 Dif Mean Queue Node 9 Sim Mean Queue Node 10 Dif Mean Queue Node 10 Sim 2.5 2 1.5 1 0.5 0 0 10 20 30 40 50 60 70 80 90 100 Mean queue lengths at nodes 9 and 10, diffusion approximation and simulation results 83

Switch 5 Mean Queue Node 17 Dif Mean Queue Node 17 Sim Mean Queue Node 18 Dif 4.5 Mean Queue Node 18 Sim 4 3.5 3 2.5 2 1.5 1 0.5 0 0 10 20 30 40 50 60 70 80 90 100 Mean queue lengths at nodes 9 and 10, diffusion approximation and simulation results 84

Switch 4 LaIn Node 8 Dif LaIn Node 8 Sim LaIn Node 9 Dif LaIn Node 9 Sim 3.5 LaIn Node 19 Dif LaIn Node 19 Sim 3 2.5 2 1.5 1 0.5 0 10 20 30 40 50 60 70 80 90 100 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

60 "case1-diffusion" "case1-simulation" "case2-diffusion" "case2-simulation" 50 40 30 20 10 0 0 500 1000 1500 2000 2500 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

100 "DIFF-init-state-25" "SIM-init-state-25" 90 "DIFF-init-state-23" "SIM-init-state-23" "DIFF-init-state-21" "SIM-init-state-21" 80 70 60 50 40 30 20 10 0 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Mean queue length, the source has 31 phases, initial phase j 0 = 21 , 23 , 25 ; [ Maglaris, et al., packet video communications ] constant service time; diffusion and simulation results. 88

12 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 10 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 8 6 4 2 0 0 100 200 300 400 500 600 700 800 900 1000 Mean queue length as a function of time, C 2 A = 1 , C 2 B is changing, ̺ = 0 . 75 . 89

2 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 1.8 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 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 100 200 300 400 500 600 700 800 900 1000 Absolute error, mean queue length, as a function of time, C 2 A = 1 , C 2 B is changing 90

200 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 180 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 160 140 120 100 80 60 40 20 0 -20 0 100 200 300 400 500 600 700 800 900 1000 Relative error, mean queue length, as a function of time, C 2 A = 1 , C 2 B is changing 91

A comment to errors Mean queue length, ̺ = 0 . 75 , t = 75 time-units, see how it depends on C 2 A , C 2 B 92

Errors as a function of time and C 2 A , C 2 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 C 2 B . • Absolute error grows negatively with the growth of coefficient of variation C 2 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 ) C 2 A ( t ) + nµC 2 B , • β ( x, t ) = λ ( t ) − nµ n − 1 < x < n . for or α 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 other side at a small (e.g. 2 − 10 ) distance ε 96

barrier barrier γ L g n − ε ( t ) n ( t ) γ R n ( t ) g n + ε ( t ) x x n − ε n n + ε n Flow balance for the barrier at x = n γ R γ L n ( t ) = g n − ε ( t ) , n ( t ) = g n + ε ( t ) , n = 1 , . . . , N − 1 with two exceptions concerning flows coming from barriers at x = 0 : g 1+ ε ( t ) = γ L 1 ( t ) + g 1 ( t ) , and x = N : g N − 1 − ε ( t ) = γ R N − 1 ( t ) + g N − 1 ( t ) . 97

The density functions f i ( x, t ; ψ i ) , i = 1 , . . . N , for the intervals x ∈ ] i − 1 , i [ are: � t f 1 ( x, t ; ψ 1 ) = φ 1 ( x, t ; ψ 1 ) + g 1 − ε ( τ ) φ 1 ( x, t − τ ; 1 − ε ) dτ , 0 � t f n ( x, t ; ψ n ) = φ n ( x, t ; ψ n ) + g n − 1+ ε ( τ ) φ n ( x, t − τ ; n − 1 + ε ) dτ + 0 � t + g n − ε ( τ ) φ n ( x, t − τ ; n − ε ) dτ, 0 n = 2 , . . . N − 1 , � t f N ( x, t ; ψ N ) = φ N ( x, t ; ψ N ) + g N − 1+ ε ( τ ) φ N ( x, t − τ ; N − 1 + ε ) dτ 0 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 ): � t γ 0 ( t ) = p 0 (0) δ ( t ) + γ ψ 1 , 0 ( t ) + g 1 − ε ( τ ) γ 1 − ε, 0 ( t − τ ) dτ , 0 � t γ L 1 ( t ) = γ ψ 1 , 1 ( t ) + g 1 − ε ( τ ) γ 1 − ε, 1 ( t − τ ) dτ , 0 � t γ L n ( t ) = γ ψ n ,n ( t ) + g n − 1+ ε ( τ ) γ n − 1+ ε,n ( t − τ ) dτ + 0 � t + g n − ε ( τ ) γ n − ε,n ( t − τ ) dτ , n = 2 , . . . N − 1 0 � t γ R n ( t ) = γ ψ n +1 ,n ( t ) + g n + ε ( τ ) γ n + ε,n ( t − τ ) dτ + 0 � t + g n +1 − ε ( τ ) γ n +1 − ε,n ( t − τ ) dτ , n = 1 , . . . N − 2 0 � t γ R N − 1 ( t ) = γ ψ N ,N − 1 ( t ) + g N − 1+ ε ( τ ) γ N − 1+ ε,N − 1 ( t − τ ) dτ 0 � t γ N ( t ) = p N (0) δ ( t ) + γ ψ N ,N ( t ) + g N − 1+ ε ( τ ) γ N − 1+ ε,N ( t − τ ) dτ. 0 99

M/M/20/20 Model 30 Intensity of output stream, Diffusion Intensity of output stream, Simulation Intensity of input stream 25 20 15 10 5 0 0 10 20 30 40 50 60 70 80 90 100 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