Perfect Sampling of Queuing Networks with Complex Routing - PowerPoint PPT Presentation

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Markovian Modelling Theorem (Markov process) If { ξ n } is a sequence of iid random variables , the process { X n } is a homogeneous discrete time Markov chain. Random Iterated system of functions The trajectory X n is the successive application of random functions taken in the set { Φ( ., ξ ) , ξ ∈ E} according a probability measure on E [Diaconis and Friedman 98] Perfect Sampling of Queuing Networks with Complex Routing 19 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Coupling Inequality Typical trajectory States Stationary version τ 0 time Coupling time After τ the two processes are not distinguishable, then stationary Scheme used to prove Markov convergence (coupling inequality) | P ( X n ∈ A ) − π A | � P ( τ � n ) Perfect Sampling of Queuing Networks with Complex Routing 20 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Forward Sampling : avoid initial state dependence Forward coupling Example 1 − p p State Steady-state ? 1 Time f 1 f 3 f 7 f 6 f 4 f 3 Always couple in the blue state Does not guarantee the steady state ! Perfect Sampling of Queuing Networks with Complex Routing 21 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect Sampling : Backward Idea Set dynamic In what state could I be at time n = 0 ? X 0 ∈ X = Z 0 ∈ Φ( X , e − 1 ) = Z 1 ∈ Φ(Φ( X , e − 2 ) , e − 1 ) = Z 2 . . . . . . ∈ Φ(Φ( · · · Φ( X , e − n ) , · · · ) , e − 2 ) , e − 1 ) = Z n Perfect Sampling of Queuing Networks with Complex Routing 22 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect sampling : Backward Idea X s s e c o P r y X r n a o i a t t S X Time X 0 All the trajectories collapse − i Z 0 = X − j Z i Z j Z − τ ∗ = { X 0 } − τ ∗ Perfect Sampling of Queuing Networks with Complex Routing 23 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect Sampling : Convergence Theorem Theorem Provided some condition on the events the sequence of sets {Z n } n ∈ N is decreasing to a single state, stationary distributed. τ ∗ = inf { n ∈ N ; Card ( Z n ) = 1 } . backward coupling time The set of possible states at time 0 is decreasing with regards to n Perfect Sampling of Queuing Networks with Complex Routing 24 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect Sampling : Coupling Condition Theorem Suppose that the set of events is finite. Then the two conditions are equivalent: τ ∗ < + ∞ almost surely; There exist a finite sequence of events with positive probability S = { e 1 , · · · , e M } such that | Φ( X , e 1 → M ) | = 1 . The sequence S is called a synchronizing pattern (synchronizing word, renovating event,...) Perfect Sampling of Queuing Networks with Complex Routing 25 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect sampling : Coupling Condition (proof) Proof ⇒ If τ ∗ < + ∞ almost surely there is a trajectory that couples in a finite time. This finite trajectory is a synchronizing pattern. ⇐ Suppose there is a synchronizing pattern with length M . Because the sequence of events is iid, it occurs almost surely on every trajectory. Applying Borel-Cantelli lemma gives the result. The forward and backward coupling time have the same distribution τ ∗ has an exponentially dominated distribution tail P ( τ ∗ > M . n ) � (1 − P ( e 1 → M )) n . Practically efficient Perfect Sampling of Queuing Networks with Complex Routing 26 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect sampling : convergence theorem (proof 1) Proof based on the shift property First, because τ < + ∞ and the ergodicity of the chain there exists N 0 s.t. | P (Φ( X , e 1 → n ) = { x } ) − π x | � ǫ. But the sequence of events is iid (stationary) then P (Φ( X , e 1 → n ) = { x } ) = P (Φ( X , e − n +1 → 0 ) = { x } ) τ ∗ < + ∞ then there exists N 1 such that P ( τ ∗ � N 1 ) � ǫ ; then P (Φ( X , e − n +1 → 0 ) = { x } ) = P (Φ( X , e − n +1 → 0 ) = { x } , τ ∗ < N 1 ) + P (Φ( X , e − n +1 → 0 ) = { x } , τ ∗ � N 1 ) , = P (Φ( X , e − τ ∗ → 0 ) = { x } , τ ∗ < N 1 ) + ǫ ′ , = P ( X 0 = x , τ ∗ < N 1 ) + ǫ ′ = P ( X 0 = x ) + ǫ ′′ . Perfect Sampling of Queuing Networks with Complex Routing 27 / 115 By the dominated convergence theorem the result follows

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect Sampling : Convergence Theorem (proof 2) Proof based on the coupling property [Haggstrom] Consider N such that P ( τ ∗ � N ) � ǫ then consider a process { X n } with the same events e − N +1 → 0 but with X − N +1 generated according π . The process { X n } is stationary. On the event ( τ ∗ < N ) we have X 0 = X 0 and P ( X 0 � = X 0 ) � P ( τ ∗ � N ) � ǫ (coupling inequality) . Finally P ( X 0 = x ) − π x = P ( X 0 = x ) − P ( X 0 = x ) � P ( X 0 � = X 0 ) � ǫ ; π x − P ( X 0 = x ) = P ( X 0 = x ) − P ( X 0 = x ) � P ( X 0 � = X 0 ) � ǫ ; and the result follows. Perfect Sampling of Queuing Networks with Complex Routing 28 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect Sampling : Algorithm Backward algorithm Trajectories Representation : transition fonction States τ∗ Xn +1 = Φ( Xn , en +1) . 1100 1010 for all x ∈ X do 1001 y ( x ) ← x 1000 end for 0110 repeat 0101 e ← Random event (); for all x ∈ X do 0100 z ( x ) ← y (Φ( x , e )); 0011 end for 0010 y ← z 0001 until All y ( x ) are equal 0000 return y ( x ) Time Convergence : If the algorithm stops, the returned −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 0 value is steady state distributed U8 U7 U6 U5 U4 U3 U2 U1 Coupling time: τ < + ∞ , properties of Φ Mean time complexity c Φ mean computation cost of Φ( x , e ) C � Card ( X ) . E τ. c Φ . Perfect Sampling of Queuing Networks with Complex Routing 29 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Perfect Reward Sampling Backward reward Trajectories Representation : transition fonction States Xn +1 = Φ( Xn , en +1) . 1100 1010 Arbitrary reward function 1001 Cost for all x ∈ X do 1000 1 y ( x ) ← x 0110 end for 0 0101 repeat 0100 e ← Random event (); for all x ∈ X do 0011 y ( x ) ← y (Φ( x , e )); 0010 end for 0001 until All Reward ( y ( x )) are equal 0000 return Reward ( y ( x )) Time Convergence : If the algorithm stops, the returned −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 0 value is steady state reward distributed U8 U7 U6 U5 U4 U3 U2 U1 τ r � τ < + ∞ Coupling time: Mean time complexity C Reward � Card ( X ) . E τ. c Φ Depends on the reward function . Perfect Sampling of Queuing Networks with Complex Routing 30 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Inverse of PDF P ( X � x ) Cumulative distribution function 1 p 3 p 2 p 1 0 1 2 3 · · · K − 1 K x Generation Inverse function algorithm s=0; k=0; Divide [0 , 1[ in intervals with length p k Find the interval in which Random falls u=random() Returns the index of the interval while u > s do Computation cost : O ( E X ) steps k=k+1 Memory cost : O (1) s=s+ p k end while return k Perfect Sampling of Queuing Networks with Complex Routing 32 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Searching optimization Optimization methods pre-compute the pdf in a table rank objects by decreasing probability use a dichotomy algorithm use a tree searching algorithm (optimality = Huffmann coding tree) Comments - Depends on the usage of the generator (repeated use or not) - pre-computation usually O ( K ) could be huge - Perfect Sampling of Queuing Networks with Complex Routing 33 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Generation : Visual representation [0 , 1] partitionning 2 12 0 1 3 12 1 2 11 12 4 7 3 1 12 12 12 12 7 12 3 4 3 12 1 6 12 12 f 1 f 2 f 3 f 4 f 5 f 6 f 7 f 8 Random iterated system of functions Function f 1 f 2 f 3 f 4 f 5 f 6 f 7 f 8 1 1 1 1 1 4 2 1 Probability 12 12 12 12 12 12 12 12 Stochastic matrix P = ⇒ simulation algorithm = RIFS Perfect Sampling of Queuing Networks with Complex Routing 34 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis The coupling problem τ estimation 1 1 1 2 2 2 1 2 0 0 1 1 0 1 0 1 Couples with probability 1 Couples with probability 1 Never couples 2 τ = 1 E τ = 2 τ = ∞ Perfect Sampling of Queuing Networks with Complex Routing 35 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis General problem Objective Given a stochastic matrix P = (( p i , j )) build a system of function ( f θ , θ ∈ Θ) and a probability distribution ( p θ , θ ∈ Θ) such that : the RIFS implements the transition matrix P , 1 ensures coupling in finite time 2 achieve the “best” mean coupling time : tradeoff between 3 - choice of the transition function according to (( p θ )), - computation of the transition Remarks Usual method | Θ | = number of non-negative elements of P = O ( n 2 ) choice in O (log n ) Perfect Sampling of Queuing Networks with Complex Routing 36 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Non sparse matrices Rearranging the system 2 12 0 1 0 1 3 12 1 2 11 12 4 7 3 1 12 12 12 12 7 12 3 4 3 12 1 6 12 12 f 1 f 2 f 3 f 4 f 5 f 6 f 7 f 8 ”Synchronizing” transformation Perfect Sampling of Queuing Networks with Complex Routing 37 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Non sparse matrices Convergence rate When at least one column is non-negative ⇒ one step coupling. The RIFS ensures coupling and the coupling time τ is upper bounded by a geometric distribution with rate � min p i , j i j number of transition functions : could be more than the number of non-negative elements at most n 2 Perfect Sampling of Queuing Networks with Complex Routing 38 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Aliasing technique Initialization Alias and threshold tables K objects while L � = ∅ do Extract k ∈ L list L= ∅ ,U= ∅ ; Extract i ∈ U for k=1; k � K; k++ do S[k]=P[k] P[k]= p k if P[k] � 1 A[k]=i K then P[i] = P[i] - ( 1 K -P[k]) U=U+ { k } ; if P[i] � 1 K then else U=U+ { i } ; L=L+ { k } ; else end if L=L+ { i } ; end for end if end while Combine uniform and alias value when rejection Perfect Sampling of Queuing Networks with Complex Routing 39 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Aliasing technique : generation 1/8 Perfect Sampling of Queuing Networks with Complex Routing 40 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Aliasing technique : generation Generation k=alea(K) 1 if Random . K � S[k] then return k else return A[k] end if Complexity Computation time : - O ( K ) for pre-computation - O (1) for generation Memory : - threshold O ( K ) (real numbers as probability) - alias O ( K ) (integers indexes in a tables) Perfect Sampling of Queuing Networks with Complex Routing 41 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Sparse matrices Rearranging the system 2 12 0 1 0 1 3 12 1 2 11 12 4 7 3 1 12 12 12 12 7 12 3 4 3 12 1 6 12 12 f 1 f 2 f 3 f 4 f 5 f 6 f 7 f 8 Aliasing transformation Complexity M = maximum out degree of states p θ uniform on { 1 , · · · , M } , threshold comparison O (1) to compute one transition Combination with “Synchronizing” techniques Perfect Sampling of Queuing Networks with Complex Routing 42 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Uniform-binary decomposition Uniform superposition 2 12 0 1 3 12 A 1 A 2 A 3 A 4 1 2 11 12 4 7 3 1 12 12 12 12 7 12 3 4 3 12 1 6 12 12 Aliasing transformation 2 3 1 1 2 1 2 1 2 1 2 1 1 2 1 3 1 1 1 1 1 1 3 3 3 2 2 1 1 3 2 3 3 3 3 4 4 4 4 1 1 1 1 1 A 1 3 A 3 A 2 A 4 Decomposition M 1 X P = Pi , Pi : stochastic matrix with at most 2 non negative elements per row M i =1 Perfect Sampling of Queuing Networks with Complex Routing 43 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Coupling property Exchange of columns or thresholds give an equivalent representative ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� Equivalent ��� ��� ��� ��� representation ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� Spanning tree Irreducibility = ⇒ there is a spanning tree going to a single state where coupling occurs. P ( τ ∗ < + ∞ ) = 1 . τ is geometrically bounded, so τ ∗ and τ ∗ C . Perfect Sampling of Queuing Networks with Complex Routing 44 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Ψ software Perfect Sampling of Queuing Networks with Complex Routing 45 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Example Random transition coefficients: Number of states 10 100 500 1000 3000 Mean coupling time 3.1 4.5 5.3 5.7 6.1 Mean execution time µ s 3 17 170 360 1100 Pentium III 700MHz and 256Mb memory. Sample size 10000. Remarks: - very small coupling time - Coefficients : same order of magnitude, aliasing enforces coupling Comparison with birth and death process : Number of states 10 100 500 1000 3000 Mean coupling time 41 557 2850 5680 17000 Mean execution time µ s 28 1800 88177 366000 3.5s Remarks: - large coupling time - sparse matrix, large graph diameter Perfect Sampling of Queuing Networks with Complex Routing 46 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Overflow model Model Parameters K servers, λ µ1 priority on overflows µ2 input rate λ , µ3 different service rate Overflow state ( x 1 , · · · , x K ), x i ∈ { 0 , 1 } , µ K size ∼ 130000 rejection low diameter non product-form structure, Coupling time distribution 0.2 Statistics Sample size : 10000 Parameter Value 0.15 minimum 113 Density 0.1 maximum 1794 median 465 0.05 mean 498 Std 180 0 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 exponential tail, low mean value Coupling time (iterations) Perfect Sampling of Queuing Networks with Complex Routing 47 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Overflow model (2) Marginal distribution Occupied servers 1 0.3 Sample size 10000 0.8 Sample size = 10000 0.25 Utilization of server 0.6 0.2 Probability 0.15 0.4 0.1 0.2 0.05 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 K 0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Index of the server Number of occupied servers Marginal probability estimation Marginal distribution of the occupied servers P ( X i = 1) Perfect Sampling of Queuing Networks with Complex Routing 48 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Overflow model (3) Reward coupling 2500 Maximum Sample size 10000 occupation Quartile 3 state 2000 Median mean value Coupling time (iterations) Quartile 1 load Minimum 1500 1000 500 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Marginal law (server number) Reward coupling time - gain 20% for the first marginals - utilization : best reduction Perfect Sampling of Queuing Networks with Complex Routing 49 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Monotonicity and perfect sampling : idea ( X , ≺ ) partially ordered set (lattice) Typically componentwise ordering on products of intervals min = (0 , · · · , 0) and Max = ( C 1 , · · · , C n ) . An event e is monotone if Φ( ., e ) is monotone on X If all events are monotone then X 0 ∈ Z n ⊂ [Φ( min , e − n → 0 ) , Φ( Max , e − n → 0 )] ⇒ 2 trajectories Perfect Sampling of Queuing Networks with Complex Routing 51 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis The Doubling Scheme Complexity Need to store the backward sequence of events Consider 2 trajectories issued from { min , Max } at time − n and test if coupling One step backward ⇒ 2 . (1 + 2 + · · · + τ ∗ ) = τ ∗ ( τ ∗ + 1) = O ( τ ∗ 2 ) calls to the transition function. Consider 2 trajectories issued from { min , Max } at time − 2 k and test if coupling Doubling step backward ⇒ 2 . (1 + 2 + · · · + 2 k ) = 2 k +2 − 2 calls to the transition function, with k such that 2 k − 1 < τ ∗ � 2 k , Number of calls : O ( τ ∗ ) Perfect Sampling of Queuing Networks with Complex Routing 52 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Monotonicity and Perfect Sampling Trajectories Monotone PS States Doubling scheme : : Maximum M n=1;R[1]=Random event; repeat n=2.n; y ( min ) ← min Generated y ( Max ) ← Max State for i=n downto n/2+1 do R[i]=Random event; end for for i=n downto 1 do 2 1 y ( min ) ← Φ( y ( min ) , R [ i ]) 0 : minimum y ( Max ) ← Φ( y ( Max ) , R [ i ]) end for −32 −16 −8 −4 −2 −1 0 until y ( min ) = y ( Max ) return y ( min ) Mean time complexity 4 C m � 2 . (2 . E τ ) . c Φ . Reduction factor : Card ( X ) . Perfect Sampling of Queuing Networks with Complex Routing 53 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Event Modelling Multidimensional state space : X = X 1 × · · · × X K with X i = { 0 , · · · , C i } . Event e : ❀ transition function Φ( ., e ); (skip rule) ❀ Poisson process λ e States Time Events e1 e2 e3 e4 1781-1840 Perfect Sampling of Queuing Networks with Complex Routing 54 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Event modelling Uniformization λ e and P ( event e ) = λ e � Λ = Λ ; e Trajectory : { e n } n ∈ Z i.i.d. sequence . ⇒ Homogeneous Discrete Time Markov Chain [Bremaud 99] X n +1 = Φ( X n , e n +1 ). Generation among a small finite space E : O (1) Perfect Sampling of Queuing Networks with Complex Routing 55 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Index Routing in Queuing Networks Index functions for event e I e For queue i i : { 0 , · · · , C i } − → O (totally ordered set) . Property : ∀ x i , x j I e i ( x i ) � = I e j ( x j ) . ex: inverse of a priority,... Routing algorithm: if x origin > 0 then { a client is available in the origin queue } x origin = x origin − 1; { the client is removed from the origin queue } j = argmin i I e i ( x i ); { computation of the destination } if j � = -1 then x j = x j +1; { arrival of the client in queue j } { in the other case, the client goes out of the network } end if end if Perfect Sampling of Queuing Networks with Complex Routing 56 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Monotonicity of Index Routing Policies Proposition If all index functions I e i are monotone then event e is monotone. Proof : Let x ≺ y two states and let be an index routing event. Let i be the origin queue for the event. j x = argmin j I e j ( x j ) and j y = argmin j I e j ( y j ) Case 1 x i = y i = 0 nothing happens and Φ( x , e ) = x ≺ y = Φ( y , e ) Case 2 x i = 0 , y i > 0 then Φ( x , e ) = x ≺ y − e i + e j y = Φ( y , e ) Case 3 x i > 0 , y i > 0 then I e j x ( x j x ) < I e j y ( x j y ) � I e j y ( y j y ) < I e j x ( y j x ); then x j x < y j x and Φ( x , e ) = x − e i + e j x � y − e i � y − e i + e j y = Φ( y , e ) Perfect Sampling of Queuing Networks with Complex Routing 57 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Monotonicity of Routing Examples [Glasserman and Yao] All of these events could be expressed as index based routing policies : - external arrival with overflow and rejection - routing with overflow and rejection or blocking - routing to the shortest available queue - routing to the shortest mean available response time - general index policies [Palmer-Mitrani] - rerouting inside queues ... Perfect Sampling of Queuing Networks with Complex Routing 58 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Monotonicity of Routing : Examples State dependent routing Stateless routing Join the shortest queue Overflow routing � x j � prio ( j ) if x j < C j ; if x j < C j ; I e I e j ( x j ) = j ( x j ) = + ∞ elsewhere; + ∞ elsewhere I e − 1 = max C j . I e − 1 = max C j . j j Join the shortest response time Routing with blocking � x j +1 � prio ( j ) if x j < C j ; if x j < C j ; I e I e j ( x j ) = j ( x j ) = µ j + ∞ elsewhere + ∞ elsewhere; I e i = max C j . I e − 1 = max C i . j i Perfect Sampling of Queuing Networks with Complex Routing 59 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Coupling Experiment Feed-forward queuing model C 1 λ 3 C λ C 0 3 λ 1 0 λ λ 2 5 λ 4 C 2 Estimation of E τ : τ 400 350 300 250 200 150 100 50 0 0 1 2 3 4 λ Perfect Sampling of Queuing Networks with Complex Routing 60 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Main Result Bound on coupling time K C i + C 2 Λ � i E τ � , Λ i 2 i =1 - Λ : global event rate in the network, - Λ i the rate of events affecting Q i - C i is the capacity of Queue i . Sketch of the proof - Explicit computation for the M / M / 1 / C - Computable bounds for the M / M / 1 / C - Bound with isolated queues Perfect Sampling of Queuing Networks with Complex Routing 61 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Explicit Computation for the M / M / 1 / C E τ b = E min( h 0 → C , h C → 0 ) λ Absorbing time in a finite Markov chain; p = λ + µ = 1 − q 0,C Level 2 q p p 0,C−1 1,C Level 3 q p q p p 0,C−2 1,C−1 2,C Level 4 q p q q p p p 0,C−3 1,C−2 2,C−1 3,C Level 5 p q q q q q p p p p p 0,1 1,2 C−2,C−1 C−1,C Level C+1 q p q q q q 0,0 C,C Level C+2 Explicit recurrence equations Case λ = µ E τ b = C + C 2 . 2 Perfect Sampling of Queuing Networks with Complex Routing 62 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Computable bounds for M / M / 1 / C If the stationary distribution is concentrated on 0 ( λ < µ ), E τ b � E h 0 → C is an accurate bound . Theorem The mean coupling time E τ b of a M/M/1/C queue with arrival rate λ and service rate µ is bounded using p = λ/ ( λ + µ ) = 1 − q. E τ b � C 2 + C Critical bound: ∀ p ∈ [0 , 1] , . 2 C ) q (1 − ( q p ) E τ b � if p > 1 C Heavy traffic Bound: 2 , p − q − . ( p − q ) 2 C ) p (1 − ( p q ) E τ b � if p < 1 C Light traffic bound: 2 , q − p − . ( q − p ) 2 Perfect Sampling of Queuing Networks with Complex Routing 63 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Computable Bounds for M / M / 1 / C Example with C = 10 120 C + C 2 100 Light traffic heavy traffic bound bound 80 C + C 2 2 60 40 20 E τ b 0 0 0.2 0.4 0.6 0.8 1 p Perfect Sampling of Queuing Networks with Complex Routing 64 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Example for tandem queues 6 = C 0 6 X 0 5 = X 0 5 0 = 5 0 4 4 3 = C 1 X 1 3 = C 1 0 = 3 2 2 X 1 1 = 2 1 1 0 0 Time Time − τ b 0 − τ b 1 ( s 0 = 2) 0 0 Coupling of Queue 0 Coupling of queue 1 conditionned by state of queue 0 6 = C 0 X 0 5 0 = 5 4 X 1 0 = 3 3 = C 1 X 1 1 = 2 2 1 0 Time − τ b 0 − τ b τ b 1 ( s 0 = 5) 1 ( s 0 = 5) 0 Then τ b � st ∞ τ b 1 + τ b 0 , normalized Perfect Sampling of Queuing Networks with Complex Routing 65 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Bound with Isolated Queues Theorem In an acyclic stable network of K M / M / 1 / C i queues with Bernoulli routing and loss if overflow, the coupling time from the past satisfies in expectation, � C i   � p i K − 1 p i (1 − ) Λ C i q i E [ τ b ] � − �   ( q i − p i ) 2 ℓ i + µ i q i − p i   i =0 K − 1 Λ � ( C i + C 2 i ) . � ℓ i + µ i i =0 Perfect Sampling of Queuing Networks with Complex Routing 66 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Conjecture for General Networks 800 700 B 2 (conjecture) 600 B 3 (conjecture) 500 400 B 1 (proven) 300 200 B 1 ∧ B 2 ∧ B 3 100 E τ b 0 λ 5 0 0.5 1 1.5 2 2.5 3 3.5 4 Extension to cyclic networks, Generalization to several types of events Application : Grid and call centers Perfect Sampling of Queuing Networks with Complex Routing 67 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Software architecture Aim of the software General architecture finite capacity queuing network simulator rare events estimation (rejection, blocking,...) statistical guarantees (independence of samples) ⇒ Simulation kernel open source (C, GPL licence) extensible library of events multiplatforms (Linux (Debian), mac OSX,...) Perfect Sampling of Queuing Networks with Complex Routing 68 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Queueing Network Description Constrained communications Description file # Number of queues µ 1 λ 3 # Queues capacities 1 1 50 # queues minimal initial state µ 2 Capacity C 0 0 0 # queues maximal initial state ν 1 1 50 # Number of events Overflow 4 # Index file - N for No index file Events types File: N # table of events type action 1 Server departure # id type rate nbq origin d1 d2 d3 d4 2 External arrival to the first empty room in the list DQ 0 2 0.8 5 -1 : 0 1 2 -1 3 Multi-server departure to DQ 4 Join the shortest queue in DQ 1 1 0.6 2 0 : -1 5 Index routing according an index table 6 Routing to the first empty room in the list DQ and overflow 2 1 0.4 2 1 : -1 7 Routing to the first empty room in the list DQ and 3 7 2.0 5 2 : 0 1 2 -1 blocking in the origin queue Perfect Sampling of Queuing Networks with Complex Routing 69 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Simulation control and output Control parameters Output # Sample number # P.S.I.2 version 4.4.4 10000 # Data Network model # Number of Antithetic variable # Number of queues 1 ... # Size of maximal trajectory # Parameters 3000000 # Sample number # Random generator seed # 10000 5 # Number of Antithetic variates # Coupling file name ... File: No file # ============= 0 [ [ 0 1 10 ] ] 1 [ [ 1 1 13 ] ] Statistical analysis 2 [ [ 1 1 2 ] ] 0.11 0.1 3 [ [ 1 1 33 ] ] Probability distribution in the buffer 0.09 0.08 ... 0.07 0.06 9999 [ [ 1 1 2 ] ] 0.05 0.04 # Size 10000 Sampling time : 0.03 0.02 0.01 3809.202000 micro-seconds 0 0 10 20 30 40 50 # Seed Value 5 Perfect Sampling of Queuing Networks with Complex Routing 70 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Example Delta interconnection network, C = 10 ρ = 0 . 9 0 8 16 24 1 9 17 25 2 10 18 26 3 11 19 27 4 12 20 28 5 13 21 29 6 14 22 30 7 15 23 31 9999 [ [ 0 2 5 7 2 8 7 4 0 7 10 3 3 2 1 5 0 0 6 3 3 6 0 3 9 1 2 4 3 1 3 6 ] ] # Size 10000 Sampling time : 4302.413600 micro-seconds Perfect Sampling of Queuing Networks with Complex Routing 71 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Monotonous reward sample First server analysis 92 [ [ 0 ] ] 93 [ [ 1 ] ] 94 [ [ 1 ] ] 95 [ [ 1 ] ] 96 [ [ 1 ] ] 97 [ [ 1 ] ] 98 [ [ 1 ] ] 99 [ [ 1 ] ] # Size 100 Sampling time : 36.230000 micro-seconds Time reduction 99 [ [ 1 1 1 ] ] # Size 100 Sampling time : 308.100000 micro-seconds Perfect Sampling of Queuing Networks with Complex Routing 72 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Coupling time study (doubling scheme) Coupling time for each queue Perfect Simulation (with doubling period) started 0 6 2 1 6 [ 1 0 9 ] 1 7 5 1 7 [ 1 0 8 ] 2 8 2 6 8 [ 1 1 7 ] 3 8 1 2 8 [ 1 0 8 ] 4 8 1 1 8 [ 1 1 9 ] 5 7 7 4 7 [ 0 1 0 ] 6 6 1 2 6 [ 1 0 10 ] 7 8 1 1 8 [ 1 1 1 ] Carefull : number of steps to couple ( τ = 2nb steps) Last queue have the largest coupling time. Perfect Sampling of Queuing Networks with Complex Routing 73 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Coupling Time Study (one step scheme) Coupling time for each queue Perfect Simulation started 0 26 1 1 26 [ 1 0 10 ] 1 58 1 1 58 [ 0 0 10 ] 2 100 2 1 100 [ 1 0 8 ] 3 91 1 51 91 [ 1 1 2 ] 4 114 1 1 114 [ 1 1 6 ] 5 210 1 1 210 [ 1 1 9 ] Distribution of the rewards coupling time Perfect Sampling of Queuing Networks with Complex Routing 74 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Download : http://gforge.inria.fr/projects/psi Perfect Sampling of Queuing Networks with Complex Routing 75 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Synthesis The psi2 unix command USAGE : psi2 unix [-ipo] argument [-hdtv] -i : input file in ext directory -p : parameter file in ext directory -o : output file in ext directory By default, output file has outputtest.txt name in ext directory -h : help. -d : With details on output file -t : Without doubling period, show coupling time of each queue -v : version Enjoy ! Perfect Sampling of Queuing Networks with Complex Routing 76 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Priority Servers Events Erlang model Event type Rate Origin Destination list Arrivals Servers Output − 1 Q 1 ; Q 2 ; Q 3 ; − 1 Arrival λ − 1 Departure µ 1 Q 1 Overflow − 1 on next free server Departure µ 2 Q 2 − 1 Departure µ 3 Q 3 Rejection if all servers are buzy Results X = { 0 , 1 } 3 Validation χ 2 test E = { e 0 , e 1 , e 2 , e 3 } K = 30 µ i decreasing Saturation probability 0 . 0579 ± 4 . 710 − 4 Simulation time 0 . 4 ms Card ( X ) = 2 K τ = 577 Perfect Sampling of Queuing Networks with Complex Routing 78 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Priority servers Saturation probability Erlang model Proba Prob 0.001 0.6 Arrivals Servers 0.5 1e−04 Output 0.4 0.3 1e−05 0.2 Overflow 0.1 on next free server 0 1e−06 λ 0 10 20 30 40 50 60 70 λ 8 8.5 9 9.5 10 10.5 11 11.5 Rejection if all servers Coupling time are buzy s µ X = { 0 , 1 } 40 400 350 µ 1 = 1 , Coupling time 300 µ 2 = 0 . 8 , 250 Reward coupling time µ 3 = 0 . 5 200 Sample size 5 . 10 6 150 Card ( X ) = 2 K 100 50 0 0 10 20 30 40 50 60 70 λ Perfect Sampling of Queuing Networks with Complex Routing 79 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Line of Servers Tandem queues Events λ C1 µ1 Event type Rate Origin Dest. list reject Arrival λ − 1 Q 1 ; − 1 C2 µ2 Routing/block µ 1 Q 1 Q 2 ; Q 1 C3 µ3 Routing/block µ 2 Q 2 Q 3 ; Q 2 · · · · · · · · · · · · C4 µ4 − 1 Departure µ 5 Q 5 µ5 C5 X = { 0 , · · · , 100 } 5 Results E = { e 0 , · · · , e 5 } C = 100 λ = 0 . 9; µ = 1 p = 1 2 Blocking probability b 1 = 0 . 34, b 2 = 0 . 02 b 3 = 0 . 02, Card ( X ) = C K b 4 , = 0 . 02. Simulation time < 1 ms Perfect Sampling of Queuing Networks with Complex Routing 80 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Multistage network Delta network Events 0 8 16 24 Event type Rate Origin Dest. list 1 9 17 25 2 10 18 26 Arrival λ − 1 Q i ; − 1 3 11 19 27 1 Routing/rejection 2 µ Q i Q j ; − 1 4 12 20 28 5 13 21 29 · · · · · · · · · · · · 6 14 22 30 Departure µ Q k − 1 7 15 23 31 X = { 0 , · · · , 100 } 32 Results E = { e 0 , · · · , e 64 } C = 100 λ = 0 . 9; µ = 1 Card ( X ) = C K Loss rate Simulation time 135 ms τ ≃ 400000 Perfect Sampling of Queuing Networks with Complex Routing 81 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Multistage network Queue length and saturation proba at level 3 E(N_34) Delta network 0.0014 7 0.0012 6 0 8 16 24 0.001 5 1 9 17 25 0.0008 4 2 10 18 26 0.0006 3 3 11 19 27 0.0004 2 4 12 20 28 1 0.0002 5 13 21 29 0 λ 0 λ 0 0.2 0.4 0.6 0.8 0.58 0.6 0.62 0.64 0.66 0.68 0.7 0.72 0.74 0.76 0.78 6 14 22 30 7 15 23 31 Coupling time X = { 0 , · · · , 100 } 32 s µ 18000 E = { e 0 , · · · , e 64 } Global coupling 16000 Reward at least 1 queue saturated (3rd level) 14000 12000 Card ( X ) = C K Reward queue 31 saturated 10000 8000 Sample size 100000 6000 4000 2000 0 0 0.2 0.4 0.6 0.8 λ Perfect Sampling of Queuing Networks with Complex Routing 82 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Resource Broker Grid model Q1 µ1 Q2 µ2 Q10 Q3 µ3 Input rates Q4 µ4 Allocation strategy State dependent allocation Resource Q5 µ5 Broker Q11 Index based routing : destination λ Q6 µ6 minimize a criteria Q7 µ7 Overflow Q12 Q8 µ8 Q9 µ9 Problem Optimization of throughput, response time,... Comparison of policies, analysis of heuristics ... Perfect Sampling of Queuing Networks with Complex Routing 83 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Routing Customers in Parallel Queues The problem: Find a routing policy maximizing the expected (discounted) throughput of the system. Several variations on this problem depend on the information available to the controller: current size of all queues (and size of the arriving batch). The applications: improve batch schedulers for cluster and grid infrastructures. Assert the value of information in such cases. Perfect Sampling of Queuing Networks with Complex Routing 84 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Index policies for routing Optimal routing policy problem is still open for n different M / M / 1 Heuristic : index policy inspired from the Multi-Armed Bandit ⇒ free parameter and compute an equilibrium point. [Mitrani 2005] for routing and repair problems. µ b λ µ Capacity C µ W S servers W is the rejection cost (free parameter). Theorem There is an optimal policy of threshold type: there exists θ such that :Reject if x � θ and accept otherwise. - θ does not depend on C as long as C > θ (including if C is infinite). - θ is a non-decreasing function of W . Perfect Sampling of Queuing Networks with Complex Routing 85 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Index policies for Routing(II) Computation of θ ( W ) linear system of corresponding to Bellman’s equation, after uniformization. 30 ’TW.txt’ 25 BestThreshold(W) 20 15 10 5 0 0 2 4 6 8 10 12 14 16 18 20 W Index function I ( x ) = inf { W | θ ( W ) = x } . Indifference case : when queue size is x , rejecting or accepting the next batch are both optimal choices if the rejection cost is I ( x ). Perfect Sampling of Queuing Networks with Complex Routing 86 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Some numerical experiments(I) 1-9_1-9_100.txt 1-9_9-9.95_100.txt 5 1.5 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.45 4.5 JSQ-mu JSQ-mu JSQ-mu2 JSQ-mu2 1.4 4 Seq Seq *optim* 1.35 *optim* 3.5 1.3 3 1.25 1.2 2.5 1.15 2 1.1 1.5 1.05 1 1 0.5 0.95 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load 9-1_1-9_100.txt 9-1_9.0-9.95_100.txt 1.8 1.4 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.7 JSQ-mu 1.35 JSQ-mu JSQ-mu2 JSQ-mu2 1.6 Seq 1.3 Seq *optim* *optim* 1.5 1.25 1.4 1.2 1.3 1.15 1.2 1.1 1.1 1.05 1 1 0.9 0.95 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load Several cases with two queues with respective parameters ( µ = 9 , 1), C = 100 Perfect Sampling of Queuing Networks with Complex Routing 87 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Some numerical experiments(III) 8-2_1-9_100.txt 8-2_9.0-9.95_100.txt 1.2 1.14 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.18 JSQ-mu JSQ-mu 1.12 JSQ-mu2 JSQ-mu2 1.16 Seq Seq 1.1 1.14 *optim* *optim* 1.12 1.08 1.1 1.06 1.08 1.04 1.06 1.04 1.02 1.02 1 1 0.98 0.98 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load 8-2_1-9_50.txt 8-2_9.0-9.95_50.txt 1.2 1.14 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.18 JSQ-mu JSQ-mu 1.12 JSQ-mu2 JSQ-mu2 1.16 Seq Seq 1.1 1.14 *optim* *optim* 1.12 1.08 1.1 1.06 1.08 1.06 1.04 1.04 1.02 1.02 1 1 0.98 0.98 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load Several cases with two queues with respective parameters ( µ = 8 , 2), C = 100 . Perfect Sampling of Queuing Networks with Complex Routing 88 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Some numerical experiments(IV) 8.1-1.9_1-9_100.txt 8.1-1.9_9.0-9.95_100.txt 1.25 1.16 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ JSQ-mu 1.14 JSQ-mu 1.2 JSQ-mu2 JSQ-mu2 Seq 1.12 Seq *optim* *optim* 1.15 1.1 1.08 1.1 1.06 1.05 1.04 1.02 1 1 0.95 0.98 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load 7-3_1-9_100.txt 7-3_9.0-9.95_100.txt 1.05 1.045 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.04 JSQ-mu JSQ-mu 1.04 JSQ-mu2 JSQ-mu2 1.035 Seq Seq *optim* *optim* 1.03 1.03 1.025 1.02 1.02 1.015 1.01 1.01 1.005 1 1 0.99 0.995 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load Some other cases Perfect Sampling of Queuing Networks with Complex Routing 89 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Some numerical experiments(V) 8-1-1_1-9_50.txt 8-1-1_9.0-9.95_50.txt 1.9 1.45 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.8 JSQ-mu 1.4 JSQ-mu JSQ-mu2 JSQ-mu2 1.7 1.35 Seq Seq *optim* *optim* 1.6 1.3 1.5 1.25 1.4 1.2 1.3 1.15 1.2 1.1 1.1 1.05 1 1 0.9 0.95 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load 7-2-1_1-9_50.txt 7-2-1_9.0-9.95_50.txt 1.3 1.25 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ JSQ-mu JSQ-mu 1.25 JSQ-mu2 1.2 JSQ-mu2 Seq Seq 1.2 *optim* *optim* 1.15 1.15 1.1 1.1 1.05 1.05 1 1 0.95 0.95 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load Three queues. Respectively, µ = 8 , 1 , 1 and µ = 7 , 2 , 1. Perfect Sampling of Queuing Networks with Complex Routing 90 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Numerical experiments(VI) 6-3-1_1-9_50.txt 6-3-1_9.0-9.95_50.txt 1.2 1.14 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ JSQ-mu JSQ-mu 1.12 JSQ-mu2 JSQ-mu2 1.15 Seq Seq 1.1 *optim* *optim* 1.08 1.1 1.06 1.05 1.04 1.02 1 1 0.95 0.98 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load 6-2-2_1-9_50.txt 6-2-2_9.0-9.95_50.txt 1.1 1.07 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.06 JSQ-mu JSQ-mu 1.08 JSQ-mu2 JSQ-mu2 1.05 Seq Seq 1.06 *optim* *optim* 1.04 1.03 1.04 1.02 1.02 1.01 1 1 0.99 0.98 0.98 0.96 0.97 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load Now, µ = 6 , 3 , 1 and µ = 6 , 2 , 2. Perfect Sampling of Queuing Networks with Complex Routing 91 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Numerical experiments(VII) 5-4-1_1-9_50.txt 5-4-1_9.0-9.95_50.txt 1.16 1.12 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.14 JSQ-mu JSQ-mu 1.1 JSQ-mu2 JSQ-mu2 1.12 Seq Seq *optim* *optim* 1.08 1.1 1.06 1.08 1.06 1.04 1.04 1.02 1.02 1 1 0.98 0.98 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load 5-3-2_1-9_50.txt 5-3-2_9.0-9.95_50.txt 1.025 1.02 Average waiting time (ratio with optim) JSQ Average waiting time (ratio with optim) JSQ 1.02 1.015 JSQ-mu JSQ-mu JSQ-mu2 JSQ-mu2 1.015 1.01 Seq Seq 1.01 *optim* *optim* 1.005 1.005 1 1 0.995 0.995 0.99 0.99 0.985 0.985 0.98 0.98 0.975 0.975 0.97 0.97 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load Now, µ = 5 , 4 , 1 and µ = 5 , 3 , 2. Perfect Sampling of Queuing Networks with Complex Routing 92 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Robustness of Index policies The index policy was computed for λ = 5 or 9 and used over the whole range λ = 1 to 10. 9-1_1-9_50.txt.1 9-1_9.0-9.95_50.txt.1 1.8 1.45 mu : 9, 1 JSQ mu : 9, 1 JSQ Average waiting time (ratio with optim) Average waiting time (ratio with optim) 1.7 c : 1, 1 JSQ-mu 1.4 c : 1, 1 JSQ-mu Nmax : 50 JSQ-mu2 Nmax : 50 JSQ-mu2 1.35 1.6 Seq Seq lambda index : 5 lambda index : 5 *optim* *optim* 1.3 1.5 1.25 1.4 1.2 1.3 1.15 1.2 1.1 1.1 1.05 1 1 0.9 0.95 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load 9-1_1-9_50.txt.2 9-1_9.0-9.95_50.txt.2 1.8 1.45 mu : 9, 1 JSQ mu : 9, 1 JSQ Average waiting time (ratio with optim) Average waiting time (ratio with optim) 1.4 1.7 c : 1, 1 JSQ-mu c : 1, 1 JSQ-mu JSQ-mu2 JSQ-mu2 Nmax : 50 Nmax : 50 1.35 1.6 lambda index : 9 Seq lambda index : 9 Seq *optim* *optim* 1.3 1.5 1.25 1.4 1.2 1.3 1.15 1.2 1.1 1.1 1.05 1 1 0.9 0.95 1 2 3 4 5 6 7 8 9 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10 Load Load Perfect Sampling of Queuing Networks with Complex Routing 93 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Adaptation to Structured Models Model Parameters set of uniformized events E = { e 1 , .., e p } global states are tuples of local states ˜ s = ( s 1 , . . . , s K ) transition function: Φ(˜ s , e i ) = ˜ r * each ˜ s ∈ X has a set of enabled events and its firing conditions and consequences Constraints Well-formed SAN models needed * exploring the subset X R (Reachable state space) State space explosion still a problem Perfect Sampling of Queuing Networks with Complex Routing 95 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Analysis of Complex Discrete Systems GRAPHICAL MODEL r = Φ(˜ ˜ s , e p ) , e p ∈ ξ s ∈ X R ˜ (STATES + TRANSITIONS) Φ(˜ s , e 1 ) Φ(˜ s , e 2 ) Φ(˜ s , e 3 ) Φ(˜ s , e 4 ) A (1) A (2) { 0;0 } { 1;0 } { 0;0 } { 0;0 } { 0;0 } { 0;1 } { 1;1 } { 0;1 } { 0;2 } { 0;1 } 0 (1) 0 (2) { 0;2 } { 1;2 } { 0;2 } { 0;2 } { 0;0 } e 4 e 5 { 1;0 } { 1;0 } { 0,2 } { 1;0 } { 1;0 } e 2 e 1 e 2 e 5 { 1;1 } { 1;1 } { 1;1 } { 1;2 } { 1;1 } 1 (1) 2 (2) 1 (2) { 1;2 } { 1;2 } { 1;2 } { 1;2 } { 1;0 } e p ∈ ξ Rates Uniformized Rates e 3 e 1 λ 1 λ 1 / ( λ 1 + λ 2 + λ 3 + λ 4 + λ 5 ) e 2 λ 2 λ 2 / ( λ 1 + λ 2 + λ 3 + λ 4 + λ 5 ) e 3 λ 3 λ 3 / ( λ 1 + λ 2 + λ 3 + λ 4 + λ 5 ) e 4 λ 4 λ 4 / ( λ 1 + λ 2 + λ 3 + λ 4 + λ 5 ) e 5 λ 5 λ 5 / ( λ 1 + λ 2 + λ 3 + λ 4 + λ 5 ) Perfect Sampling of Queuing Networks with Complex Routing 96 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis New solutions for huge SAN models Monotonicity and Perfect Simulation Idea Monotonicity property for SAN related to the analysis of structural conditions * component-wise state space formation Families of SAN models SAN models with a natural partial order (canonical) * e.g. derived from Queueing systems models [Vincent 2005] SAN models with a given component-wise partial order (non-lattice) Perfect Sampling of Queuing Networks with Complex Routing 97 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Partially Ordered State Spaces Canonical component-wise ordering 00 Queueing Network Model Equivalent MC λ µ ν 10 Event e 1 K 1 K 2 Event e 2 20 01 Equivalent SAN Model Event e 12 A (1) A (2) 11 Type Event Rate 0 0 λ loc e 1 21 02 syn e 12 µ e 12 e 1 e 2 e 12 loc e 2 ν . . . . . . 12 e 12 e 1 e 2 e 12 K 1 K 2 22 03 e 1 e 12 13 23 Perfect Sampling of Queuing Networks with Complex Routing 98 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Partially Ordered State Spaces Non-lattice component-wise ordering Find a partial order of X demands a high c.c. Possible to find extremal global states in the underlying chain * |X M | states: more than two extremal states Complexity: related to τ , but also |X M | Extremal states Component-wise formation has ordered state indexes * consider an initial state composing X M * add to X M the states without transitions to states with greater indexes Perfect Sampling of Queuing Networks with Complex Routing 99 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Non-lattice component-wise ordering Resource sharing model with reservation (Dining Philosophers) Type Event Rate loc lt i µ P (1) P ( K ) P ( i ) P ( K − 1) syn tr i λ rl 2 tr K − 1 tr i − 1 tr K − 2 syn rl i λ lr K rl 0 rl i +1 tl K T (0) T ( K ) T ( i ) T ( K − 1) loc rt K µ syn tl K λ syn lr K λ tr 1 lt 1 rt K tl K tr i lt i tr K − 1 lt K − 1 lr K tr i − 1 tr K − 2 R (0) L (0) R ( K ) L ( K ) R ( i ) L ( i ) R ( K − 1) L ( K − 1) rl 0 rl 1 lr K rl i rl K − 1 i = 2 . . . ( K − 2) Perfect Sampling of Queuing Networks with Complex Routing 100 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Non-lattice component-wise ordering e.g. three philosophers with resources reservation, graphical model of the underlying transition chain, extremal states identification 000 001 010 100 002 011 101 020 200 012 102 201 Perfect Sampling of Queuing Networks with Complex Routing 101 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis SAN Perfect Simulation Resource sharing model with reservation- K Philosophers X R X M K X PEPS* (iteration) Perfect PEPS* (sample) 8 6,561 985 43 0.003185 sec. 0.032354 sec. 10 59,049 5,741 111 0.038100 sec. 0.111365 sec. 12 531,441 33,461 289 0.551290 sec. 0.689674 sec. 14 4,782,969 195,025 755 5.712210 sec. 2.686925 sec. 16 43,046,721 1,136,689 1,975 68.704325 sec. 15.793501 sec. 18 387,420,489 6,625,109 5,169 —- 83.287321 sec. Numerical results 3.2 GHz Intel Xeon processor under Linux, 1 GByte RAM times: for one iteration on PEPS and for one sample generation on Perfect PEPS Remarks: X contraction in |X M | X limitation 6 × 10 7 onPEPS overcame Perfect Sampling of Queuing Networks with Complex Routing 102 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Non Monotonic Systems Classic non monotonous events Batch arrivals Batch services Join procedures Negative customers · · · Almost monotonous events Monotonicity is not verified at the frontier Batch arrivals : batch rejection (queue full) Batch services : system almost empty Join procedures : system almost empty Negative customers : system almost empty Perfect Sampling of Queuing Networks with Complex Routing 103 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Envelopes Aim : Build a monotonous upper and lower process Hypothesis : X is lattice, T = max X and B = min X U ( M , m , e ) = sup Φ( s , e ); m � s � M L ( M , m , e ) = m � s � M Φ( s , e ) . inf Remark : if e is monotonous U ( M , m , e ) = Φ( M , e ) and L ( M , m , e ) = Φ( m , e ) Bounding process : Y 0 = T and Z 0 = B ; Y n = U ( Y n − 1 , Z n − 1 , e 1 → n ) and Z n = L ( Y n − 1 , Z n − 1 , e 1 → n ) . Remark : Y n and Z n are not Markov chains but the couple is Markov. Y n � X n � Z n . Perfect Sampling of Queuing Networks with Complex Routing 104 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Envelopes (2) Theorem (Convergence) Either the process ( Y n , Z n ) hits the diagonal in finite time with probability 1 or nether hits the diagonal. When the process hits the diagonal the value is stationary distributed. Problems ( P 1 ) The assumption that S n hits the diagonal may not be verified. ( P 2 ) Even if convergence theorem, the coupling time could become prohibitively large. ( P 3 ) The time needed to compute U ( M , m , e ) and L ( M , m , e ) might depend on the number of states between m and M Perfect Sampling of Queuing Networks with Complex Routing 105 / 115

System Simulation Perfect Sampling Discrete Time Markov Chain Event Driven Simulation Case Studies Advanced Topics Synthesis Variance Reduction Coupled trajectories U n · · · U 5 U 4 U 3 U 2 U 1 U 0 StatesReward R 1 / A 0 StatesReward R 2 / A 0 StatesReward R A / A 0 Coupled samples hope : negatively correlated Perfect Sampling of Queuing Networks with Complex Routing 106 / 115