Queueing Theory and Stochastic Processes
An academic course applying queueing systems, Poisson and Markov processes, performance metrics, capacity planning, and discrete-event simulation to computing systems.
Queueing theory analyzes waiting, service, and completion in systems with variable arrival and service processes. Stochastic processes provide their time-dependent probabilistic foundation. Models are applied to task queues, connection pools, server capacity, latency distributions, and simulation of computing systems.
Unit 1: Queueing Systems and Performance Measures
1.1. Arrival, Service, and Waiting Times
A queueing system comprises an arrival process, waiting discipline, server count, and service-time distribution. Arrival rate, service rate, average jobs in system, and response time are distinct measures. Response time includes both waiting and actual service.
1.2. Kendall Notation and System Structure
Kendall notation summarizes arrival distribution/service distribution/number of servers. Queue capacity and customer class can be appended where necessary. A single mean service rate need not characterize service-time variability across workloads.
W = W[q] + E[S]; L = λ[eff]·W1.3. Latency Distribution and Measurement Discipline
Measure arrival, service-start, and completion timestamps separately. Mean response time is not the same as p95/p99 latency; censored and timed-out requests should not be silently discarded.
1.4. Numerical and Algorithmic Example
If throughput is three jobs per second and mean work-in-system is six jobs, Little’s law gives a mean sojourn time of two seconds in steady state. Use admitted throughput, not offered traffic rejected before entry.
1.5. Arrival, Service and Departure Processes
In a queueing system, arrivals, service requirements and server capacity are distinct random processes. Arrival rate need not equal completed throughput when jobs are blocked or lost. Waiting time and in-service time are different measurements. Under appropriate long-run conditions, Little's law relates mean population and mean sojourn time for the same system boundary; measurements from different boundaries must not be mixed.
L = λₑff·W; W = Wq + E[S]1.6. Kendall Notation and Modeling Assumptions
Kendall notation A/S/c encodes interarrival and service distributions and the number of servers. M/M/1 assumes exponential interarrivals, exponential service and one server. Queue discipline, buffer capacity and calling population require explicit assumptions when not specified. An M/M/1 formula cannot be transferred to multiserver, priority or finite-buffer systems based on superficial similarity.
1.7. Measurement Boundaries and Latency Percentiles
A service may contain separate admission, worker, database-connection and downstream waiting queues. End-to-end latency can be decomposed into stages, but the percentile of a sum is generally not the sum of the component percentiles. Low average delay can coexist with high p99 latency. Sampling windows, clock synchronization and trace correlation determine the validity of observed measurements.
1.8. Computational Example and Calculation
If the system admits 2 jobs per second and mean sojourn time is 3 seconds, its long-run mean population is 6 jobs. Population includes both waiting and in-service jobs.
L = 2·3 = 6 jobs1.9. Independent Verification and Implementation Constraints
Measurements should consistently use admitted throughput and completed-job sojourn times, excluding rejected requests.
Unit 2: Random Variables and Service-Time Distributions
2.1. Mean Service Time and Variance
Service time is modeled as a nonnegative random variable. Its expectation describes mean work while variance affects delay fluctuations. The exponential distribution is memoryless, a property not shared by every real CPU or I/O workload.
2.2. Exponential, Erlang, and Heavy-Tailed Models
Deterministic, exponential, Erlang, and heavy-tailed service distributions can lead to different queue behavior. A large second moment increases waiting even if mean service time is unchanged. Mean service rate alone is inadequate for model selection.
E[S²] = Var(S) + (E[S])²2.3. Distribution Fit and Sampling
Report empirical CDFs, median, mean, variance, and high percentiles. Short observation windows may inadequately estimate rare long-service events.
2.4. Numerical and Algorithmic Example
Two service distributions may both average 10 ms while one is deterministic and another contains rare requests longer than 100 ms. Different second moments imply different waiting behavior.
2.5. Service-Time Variability and Second Moments
Two workloads with identical mean service time can exhibit very different queueing delays. Variance and the second moment are decisive in M/G/1 waiting-time analysis. Squared coefficient of variation summarizes relative service variability. Heavy-tailed workloads may contain rare but extremely long jobs that dominate averages and tail latency. Retaining only sample mean is insufficient.
E[S²] = Var(S) + (E[S])²; cₛ² = Var(S)/(E[S])²2.6. Exponential, Deterministic and Heavy-Tailed Service
Exponential service is memoryless, whereas deterministic service has a fixed duration. For equal means, deterministic variance is zero while exponential variance equals the squared mean. Parsing or media-processing workloads often vary with input size and need not follow an exponential law. Distributional fit should use quantiles, moments and class-conditional distributions rather than a convenient label alone.
2.7. Job Classes and Mixture Distributions
A mixture of short and long jobs need not behave like a single homogeneous service class. Class probabilities, means and variances determine mixture moments, including between-class variability. Priority or separate-queue design decisions can be analyzed from this mixture. Without job-class labels, observed production distributions are difficult to interpret causally.
2.8. Computational Example and Calculation
With mean service 0.1 seconds and second moment 0.03 seconds squared, variance equals 0.02 and squared coefficient of variation equals 2, higher than the exponential service value of 1.
Var(S)=0.03−0.01=0.02; cₛ²=22.9. Independent Verification and Implementation Constraints
Mean service time alone cannot characterize the waiting-time effect of service variance.
Unit 3: Poisson Processes and Arrival Models
3.1. Independent Increments and Event Counts
A homogeneous Poisson process has independent increments, and arrival counts over a time interval follow a Poisson law with parameter proportional to its length. Consecutive interarrival times are independent exponentials. Bursts and nonstationary rates can invalidate a homogeneous model.
3.2. Time-Varying Intensity and Batch Arrivals
Nonhomogeneous arrivals use intensity λ(t). Batch releases from upstream services may violate independent Poisson assumptions. Superposition and thinning results require their corresponding independence conditions.
P{N(t)=k} = exp(−λt)·(λt)ᵏ/k!3.3. Statistical Checks of Poisson Assumptions
Arrival counts are collected in equal intervals; compare mean–variance relations, autocorrelation, and time-varying intensity. If Poisson assumptions fail, simulate using a more general arrival process.
3.4. Numerical and Algorithmic Example
Under independent Poisson arrivals averaging six per minute, the probability of zero arrivals during one minute is e⁻⁶. Burst arrivals or synchronized job releases require a different arrival model.
3.5. Properties of Poisson Counting Processes
A homogeneous Poisson counting process has independent increments whose distribution depends only on interval length. Over a small interval, the one-arrival probability is approximately λΔt, with multiple arrivals of smaller order. Exponential interarrival times follow from these assumptions. Batch arrivals, periodic jobs and synchronized retries can violate independence in real systems.
P{N(t)=k} = e^(−λt)·(λt)ᵏ/k!3.6. Superposition and Thinning
Superposition of independent Poisson processes is Poisson with the sum of their rates. Independently retaining arrivals with a fixed probability produces another Poisson process through thinning. Dependent sources or time-varying selection can invalidate the simple result. Correlation between workload origins must be measured when modeling server input.
3.7. Time-Varying Arrival Intensities
Arrival intensity may vary substantially over time. A nonhomogeneous Poisson model uses λ(t), with expected event count given by its integral over the interval. Using a single global average can understate capacity needs at peak periods. Estimation from historical arrivals should distinguish time-of-day seasonality, workday effects and exceptional bursts.
E[N(b) − N(a)] = ∫ₐᵇ λ(t) dt3.8. Computational Example and Calculation
In a homogeneous Poisson process of rate 2 per second, the probability of no arrivals over 3 seconds equals e⁻⁶. Batched arrivals would require a different arrival model.
P{N(3)=0}=e⁻⁶ ≈ 0.002483.9. Independent Verification and Implementation Constraints
Poisson assumptions require testing independence rather than inferring it from the average rate alone.
Unit 4: Discrete-Time Markov Chains
4.1. Transition Matrices and the Markov Property
In a Markov chain, the next-state distribution depends on the present state but not the earlier history once the current state is known. Each transition-matrix row is a probability distribution. Existence of a stationary distribution differs from convergence toward it from every starting state.
4.2. State Design and Transition Probabilities
States may encode job count, service status, or load levels. An insufficient state summary may not be Markovian; residual service time may need to be incorporated.
πP = π; ∑[i] πᵢ = 14.3. Stationarity and Convergence Checks
Validate stochastic matrix rows, accessibility, periodicity, and recurrence conditions. Compare simulated long-run state frequencies to theoretical stationary probabilities where valid.
4.4. Numerical and Algorithmic Example
For P = [[0.8, 0.2], [0.4, 0.6]], stationary probabilities are (⅔,⅓). They characterize measured state shares only if the system’s actual transition mechanism matches the assumed matrix.
4.5. Transition Matrices and the Markov Property
In a discrete-time Markov chain, the next-state distribution depends only on the present state given that state. Rows of a transition matrix are nonnegative and sum to one. Multistep transitions follow from matrix powers. Omitting job classes or residual service times from the state may destroy the Markov property for the real queueing system.
P⁽ⁿ⁾ = Pⁿ; ∑ⱼPᵢⱼ = 14.6. Stationary Distributions and Irreducibility
A stationary distribution satisfies πP=π. A finite irreducible chain has a unique stationary distribution; convergence from arbitrary initial states also depends on periodicity conditions. Reducible models can exhibit multiple long-run regimes associated with distinct closed classes. Computed stationary vectors must be nonnegative and normalized.
4.7. Return Times and Communicating Classes
Mutual reachability partitions states into communicating classes. Recurrent, transient and absorbing classes determine long-run behavior. For a positive recurrent state, mean return time is related to the reciprocal stationary probability. A queueing model with an absorbing failure or shutdown state cannot be analyzed correctly by ignoring that class while claiming a global stationary regime.
4.8. Computational Example and Calculation
For transition matrix [[0.8,0.2],[0.4,0.6]], stationary probabilities are (⅔,⅓). Multiplying this row vector by P must reproduce it.
π=(⅔,⅓); πP=π4.9. Independent Verification and Implementation Constraints
Verification additionally checks nonnegativity and unit normalization.
Unit 5: Continuous-Time Markov Chains and Birth–Death Processes
5.1. Generator Matrices and Transition Rates
Continuous-time Markov chains use an infinitesimal generator. Off-diagonal rates are nonnegative and each row sums to zero. Birth–death processes have transitions only between neighboring states.
5.2. Birth–Death Models
For queue length, births correspond to arrivals and deaths to completions. A stationary distribution exists only under suitable stability conditions. Finite capacity and customer abandonment change the transition rates.
Qᵢᵢ = −∑[j≠i] Qᵢⱼ; πQ = 05.3. Stability and Boundary States
Simulation distinguishes exponential holding times from actual event timestamps. Verify generator row sums, including boundary states.
5.4. Numerical and Algorithmic Example
In M/M/1, transitions n→n+1 occur at rate λ and n→n−1 at rate μ for n>0. No service completion transition exists at state zero. Mishandling this boundary distorts the generator.
5.5. Generator Matrices and Transition Rates
In a continuous-time Markov chain, off-diagonal entries of generator Q are transition rates; diagonal entries make row sums zero. Short-interval transition probabilities follow from a first-order approximation based on Q. Rates are not probabilities and have inverse-time units. Mixing seconds and minutes can invalidate the entire predicted latency scale.
P(Δt) ≈ I + Q·Δt; Qᵢᵢ = −∑ⱼ≠ᵢQᵢⱼ5.6. Balance in Birth–Death Processes
In a birth–death chain where state is system population, transitions n→n+1 represent arrivals and n→n−1 service completions. Constant rates produce the special M/M/1 case. Local balance can determine stationary probabilities subject to the correct boundary conditions. State-dependent capacity requires state-dependent death rates; one constant μ may not describe every population level.
5.7. Numerical Evaluation of Matrix Exponentials
For a finite homogeneous chain, the transition matrix after time t is e^(Qt). General matrix exponential routines suit small systems, while methods such as uniformization can help with large sparse generators. The resulting matrix should have approximately unit row sums and nonnegative entries. These properties provide valuable numerical implementation invariants.
P(t) = exp(Qt)5.8. Computational Example and Calculation
For birth rate λ=2 and death rate μ=3 in the standard single-server birth–death model, utilization is ⅔ and empty-system probability is ⅓.
ρ=2/3; π₀=1−ρ=1/35.9. Independent Verification and Implementation Constraints
The result assumes a constant death rate at positive populations; state-dependent service rates require different balance equations.
Unit 6: M/M/1 and Multi-Server Queueing Models
6.1. Single-Server Queue Stability
M/M/1 assumes Poisson arrivals, exponential service, one server, and sufficient buffer capacity. Stability requires an arrival rate below the service rate. Mean response time grows sharply as utilization approaches one.
6.2. Multi-Server Service Models
M/M/c has c parallel servers with aggregate service capacity c times the per-server rate. Waiting probability can be calculated using Erlang C under its assumptions on server homogeneity and service discipline.
ρ = λ/μ < 1; W(M/M/1) = 1/(μ−λ)6.3. Comparing Analytical Results with Measurements
Compare analytical predictions to measurement only after testing arrival/service assumptions, adequate steady-state sampling, and server behavior. Matching mean response times does not establish correct tail predictions.
6.4. Numerical and Algorithmic Example
For one server with μ=10/s and λ=8/s, utilization is .8 and mean sojourn is .5 seconds. At λ=9/s it becomes one second; queue delay responds nonlinearly to utilization.
6.4.1. Illustrative Pseudocode
rho = arrival_rate / service_rate
assert rho < 1 # stable M/M/1 hypothesis
mean_system_time = 1 / (service_rate - arrival_rate)
mean_system_jobs = arrival_rate * mean_system_time
compare_to_steady_state_measurements()6.5. M/M/1 Steady-State Distribution
An M/M/1 queue with λ<μ has a geometric stationary system-population distribution with utilization ρ=λ/μ. At λ≥μ, the standard infinite-buffer model has no finite-mean stationary regime. Mean population is ρ/(1−ρ) and mean sojourn time is 1/(μ−λ). Waiting time in queue and total sojourn time must be distinguished.
P{N=n} = (1−ρ)ρⁿ; L = ρ/(1−ρ)6.6. M/M/c and Erlang C Waiting Probability
For M/M/c, λ is arrival rate, μ is service rate per server, and c is the number of parallel servers. Stability requires ρ=λ/(cμ)<1. Erlang C gives the probability that an arrival must wait under the standard infinite-buffer FCFS exponential-service assumptions. If workers share contended physical resources, treating them as independent servers can be misleading.
a = λ/μ; ρ = a/c; ρ < 16.7. Delay Near Capacity Saturation
As utilization approaches one, M/M/1 mean delay grows hyperbolically rather than linearly. Adding workers does not necessarily scale service capacity proportionally: memory bandwidth, connection pools and lock contention may become bottlenecks. Capacity tests should measure whether per-worker service time changes under load and evaluate whether the queueing assumptions remain applicable.
6.8. Computational Example and Calculation
For an M/M/1 server with λ=8 and μ=10 jobs per second, utilization is 0.8, mean sojourn time 0.5 seconds and mean waiting time in queue 0.4 seconds.
ρ=0.8; W=1/(10−8)=0.5 s; Wq=0.4 s6.9. Independent Verification and Implementation Constraints
This calculation does not establish that a real CPU and database pool behave like one independent exponential server.
Unit 7: General Service Distributions and M/G/1 Systems
7.1. Second Moment and Mean Waiting
M/G/1 retains Poisson arrivals but permits general service times. For FCFS and suitable finite-moment and stability assumptions, the Pollaczek–Khinchine formula yields mean waiting. Its dependence on the second moment explains strong variability effects.
7.2. Service Variability and FCFS
Deterministic and highly variable service can produce different waiting times even at equal mean service. Correlated service or priority scheduling invalidates direct use of the basic FCFS expression.
W[q] = λ·E[S²]/(2(1−ρ)); ρ = λ·E[S] < 17.3. Distribution Estimation and Validity Conditions
Estimating service-time variability requires sufficient observations, particularly in the tails. Queue waiting and processing time are recorded separately.
7.4. Numerical and Algorithmic Example
Let mean service be .1 s, second moment .02 s², and arrival rate eight per second. With utilization .8, the M/G/1 FCFS formula gives W[q]=8·0.02/(2·0.2)=0.4 s. Ignoring the second moment would give an incorrect prediction.
7.5. Pollaczek–Khinchine Mean Waiting Time
M/G/1 assumes Poisson arrivals, one server and independent generally distributed service times. Under FCFS and stability, mean waiting time depends on the second moment of service time. Equal mean service times can therefore produce different delays. If the second moment is infinite, the standard formula cannot provide finite mean waiting time; a mean-only M/M/1 substitution is inadequate.
Wq = λ·E[S²] / (2(1−ρ)); ρ = λ·E[S] < 17.6. Impact of Service Variability
Deterministic service has second moment E[S]², while exponential service has second moment 2E[S]². At equal λ and mean service, mean queue waiting in M/D/1 is half that of M/M/1. This highlights why file size and job type affect service-time distributions. Heavy tails and insufficient sampling can make the estimated second moment unreliable.
7.7. Independence and Service-Order Assumptions
Pollaczek–Khinchine allows general service distributions but retains assumptions including independence and FCFS. Its mean-wait formula does not directly apply to priority scheduling or processor sharing. Workloads with clusters of large jobs may exhibit service-time correlation. Autocorrelation, quantiles and class-specific measurements should guide model selection.
7.8. Computational Example and Calculation
For M/G/1 with λ=8, E[S]=0.1 seconds and E[S²]=0.02 seconds squared, utilization is 0.8 and Wq=0.4 seconds; total W=0.5 seconds.
Wq = 8·0.02/(2·0.2)=0.4 s; W=0.5 s7.9. Independent Verification and Implementation Constraints
Uncertainty in the second moment estimated from samples propagates into Wq predictions.
Unit 8: Little’s Law and Renewal Processes
8.1. Work-in-System and Throughput
Little’s law relates mean work-in-system, effective throughput, and mean sojourn time under suitable stability assumptions. Offered arrival rate can differ from effective admission or departure rate when jobs are rejected or lost. The law does not require Poisson arrivals.
8.2. Renewal–Reward Theorem
Renewal processes calculate long-run rates from independent, identically distributed cycles. The renewal-reward theorem divides expected reward per cycle by expected cycle duration. Partially observed cycles at measurement boundaries need explicit treatment.
L = λ[eff]·W; lim[t→∞] R(t)/t = E[R]/E[T]8.3. Observation Windows and Population Consistency
To verify Little’s law, estimate L, effective arrival rate, and W for the same population and observation window. Mixing classes or boundaries can yield misleading discrepancies.
8.4. Numerical and Algorithmic Example
If average population is 20 jobs and effective throughput five jobs per second, Little’s law gives four seconds of mean response. Numerator and denominator must refer to the same population; partial traces need care.
8.5. System Boundaries of Little’s Law
Little's law equates long-run mean population to effective throughput multiplied by mean time in the same system. Counting rejected arrivals in λ while measuring only accepted jobs for W and L produces inconsistent results. Jobs present before or remaining after an observation window need appropriate treatment. The law is not an instantaneous identity without suitable stationarity and flow-balance conditions.
L = λₑff·W; Lq = λₑff·Wq8.6. Renewal–Reward Theory
Renewal processes partition a system into appropriate independent identically distributed cycles, allowing long-run reward rates to be computed. Expected reward per cycle divided by expected cycle length yields the time average. Empty-to-empty busy cycles can be renewal points in suitable queue models. Time-dependent forcing or persistent state can violate cycle independence.
8.7. Residual Life and Observation Bias
A job seen in progress at a random time tends to be longer than a randomly sampled job at its start: this is length bias, linked to the second moment. Observing the server busy at a random time and at an arrival epoch are equivalent only under conditions such as PASTA. The observation protocol must specify which epochs are sampled.
E[Sʳ] = E[S²] / (2E[S]) (stationary residual life)8.8. Computational Example and Calculation
At throughput 5 jobs per second and mean sojourn 2 seconds, Little’s law predicts 10 jobs in the stated system boundary, not merely the waiting line.
L=5·2=10 jobs8.9. Independent Verification and Implementation Constraints
When measuring admission and processing stages separately, apply matched λ and W at each boundary.
Unit 9: Queue Disciplines, Priorities, and Admission Control
9.1. Service Disciplines and Fairness
FCFS, priority, shortest-job-first, and round-robin policies trade mean latency against fairness in different ways. Reordering work can affect class-specific waiting even if service demand remains unchanged. Starvation and aging policies require analysis.
9.2. Admission Control and Backpressure
Admission control can bound queue buildup by rejecting work before unlimited accumulation. Backpressure asks upstream systems to reduce production. Drop rates and latency among completed jobs must be evaluated together.
λ[eff] = λ(1−P[block]); U = T[busy]/T[observed]9.3. Class-Specific Latency and Service Share
Measure latency percentiles and actual service share for each priority class. A lower global mean can mask unacceptable starvation of low-priority jobs.
9.4. Numerical and Algorithmic Example
With queue capacity 100, overload can cause rejections. The p99 latency of admitted requests might improve while rejection rate worsens. Reporting latency only for completed requests can misrepresent service quality.
9.5. FCFS, LCFS and Processor Sharing
FCFS serves jobs in arrival order, LCFS favors newer jobs, and processor sharing idealizes concurrent division of service capacity. These disciplines can yield different individual delays and response-time distributions for the same arrivals. Some means coincide under special assumptions, but that does not make disciplines equivalent for every metric. Policy selection must account for fairness and deadline requirements.
W = waiting time + service time9.6. Priority Queues and Starvation
Fixed high-priority classes can cause severe delays or starvation for low-priority jobs. Preemptive and nonpreemptive priority models differ, and whether service may be interrupted must be specified. Aging can reduce starvation by gradually increasing priority, but changes the assumptions of FCFS formulas. Delay and throughput should be reported per job class.
9.7. Finite Buffers, Loss and Admission Control
Finite-capacity queues may reject or drop arrivals when full. Effective throughput counts admitted traffic and need not equal the external offered arrival rate. An M/M/1/K stationary distribution differs from infinite-buffer M/M/1. Admission-control analysis must include loss probability and fairness, not just reduced delay.
λₑff = λ·(1 − Pblock)9.8. Computational Example and Calculation
An offered arrival rate of 10 requests per second with 20 percent blocking produces admitted throughput 8 per second, before accounting for retries.
λₑff = 10·(1−0.2)=89.9. Independent Verification and Implementation Constraints
The mean waiting-time population must be identified as admitted jobs, not all attempted requests.
Unit 10: Queueing Networks and Distributed-System Performance
10.1. Serial and Parallel Service Stages
Queueing networks interconnect multiple service centers that route work between stages. End-to-end latency combines network delay, service, and waiting. Fork–join completion depends on the slowest parallel branch, so independent mean latencies cannot always simply be summed.
10.2. Open and Closed Queueing Networks
Open and closed networks differ in population and arrival assumptions. Product-form results such as Jackson networks require specified Markovian routing and service conditions. Database connection pools and worker counts can form separate service centers.
W[end] = ∑[k=1..K] Wₖ (serial stages of the same job)10.3. Bottlenecks and End-to-End Latency
Inspect utilization and latency distribution for bottlenecks, and counts at each pipeline stage. Ignoring correlations and feedback loops may produce biased predictions.
10.4. Numerical and Algorithmic Example
For the same job traversing two serial stages with mean times 10 and 20 ms, mean total is 30 ms by linearity of expectation. For parallel branches completion is max(T₁,T₂), not the sum, and means alone do not determine its distribution.
10.5. Open Queueing Networks and Traffic Equations
An open queueing network receives external traffic, routes jobs probabilistically among stations and eventually lets them depart. Each node's effective arrival rate includes both external input and routed work from other nodes. Jackson product-form results depend on specific Poisson, exponential-service and routing assumptions. Shared resources or state-dependent routing can break that structure.
Λⱼ = αⱼ + ∑ᵢ Λᵢpᵢⱼ10.6. Serial Stages and End-to-End Delay
For a job visiting sequential queues, end-to-end sojourn time is the sum of realized stage times. Means add linearly, but p95 and p99 percentiles generally do not. Revisits and routing cycles contribute additional delay. Measurements need clear stage boundaries and must separate original arrivals from retry-generated traffic.
10.7. Fork–Join and Parallel-Task Completion
When a job forks into parallel tasks and waits for every branch to finish, completion is governed by the slowest branch. This maximum matters even if branch delays are dependent. Consequently, parallelism can reduce mean task time yet worsen end-to-end tail latency. Standard Jackson product-form assumptions do not directly apply to general fork–join synchronization.
Wfork-join = max(W₁, W₂, …, Wₖ)10.8. Computational Example and Calculation
Station A receives 2 jobs per second externally, and each job visits B with probability 0.5. Without other arrivals or feedback, B receives 1 job per second on average.
Λᴮ = 2·0.5 = 1 job/s10.9. Independent Verification and Implementation Constraints
Node capacity checks use effective routed arrivals, not the external arrival rate copied to each station.
Unit 11: Discrete-Event Simulation and Stochastic Modeling
11.1. Event Calendar and Logical Clock
Discrete-event simulation advances the logical clock from one arrival, service start, or completion event to the next. A priority queue orders events by timestamp with a deterministic tie-breaker. Event-driven simulation may avoid unnecessary fixed time steps.
11.2. Random Seeds and Warm-Up
Random samples use declared seeds and distribution parameters. Warm-up bias, replication count, and autocorrelation affect confidence intervals. A tractable M/M/1 model provides a useful verification target for the simulator.
t[next] = min{t[arrival], t[completion], t[timeout], …}11.3. Simulation Invariants and Confidence Intervals
Maintain invariants such as nonnegative queue length and completed jobs having entered service. Compare sample variance and tail-quantile estimates across independent replications.
11.4. Numerical and Algorithmic Example
If arrival and completion both occur at time 1.0 s, the event calendar must use a documented tie-breaking rule. Inconsistent ordering can change admission decisions. The simulated queue length must remain nonnegative.
11.4.1. Illustrative Pseudocode
events = priority_queue(initial_arrivals)
while events and clock < horizon:
event = events.pop_min()
clock = event.time
update_integrals_of_queue_length()
apply_arrival_or_completion(event)
schedule_followup_events()
check_nonnegative_population()11.5. Event Calendar and Simulation Clock
Discrete-event simulation schedules arrivals, completions and timeouts in a time-ordered priority queue. The simulation clock jumps to the next event rather than iterating over every time tick. Ties among events at equal timestamps may affect queue population, so deterministic ordering is required. Cancellation, draining and rollback events must preserve state invariants.
tₙₑₓₜ = min(tarrival, tcompletion, ttimeout)11.6. Warm-Up Periods and Initialization Bias
Starting simulations with an empty queue can bias early observations toward lower waiting times than stationary behavior. Warm-up deletion aims to reduce initialization bias. The required warm-up depends on run length and mixing behavior. For an unstable queue, choosing a longer warm-up does not repair the model; load and service capacity must be revisited.
11.7. Replications, Confidence Intervals and Randomness
One random seed is not enough for reliable population inference from simulation. Independent replications use different seeds and support confidence intervals for means. Common random numbers can reduce variance when comparing algorithms, provided the experiment is correctly paired. Uncertainty for mean estimates and tail percentiles should be reported separately.
CI(mean) ≈ x̄ ± t·s/√n (independent replication means)11.8. Computational Example and Calculation
With 10 independent simulation replications, a confidence interval for long-run mean delay should be based on replication means rather than pretending individual correlated jobs are independent.
CI = x̄ ± t₉,0.975·s/√1011.9. Independent Verification and Implementation Constraints
Use distinct seeds while recording warm-up and run horizons for every replication.
Unit 12: Queue Stability, Capacity Planning, and Model Validation
12.1. Traffic Intensity and Stability Conditions
Stability concerns preventing unbounded queue growth over the long run. Mean service capacity exceeding arrivals is not universally sufficient for multi-class or feedback-rich networks. Loss, abandonment, and finite buffers require explicit representation.
12.2. Balancing Capacity and Latency Targets
Capacity planning considers peak arrivals, service variability, and tail-delay objectives rather than relying only on average load. Infrastructure costs may be balanced against queueing cost. High CPU utilization alone is not proof of efficient operation.
ρ = λ·E[S]/c; ρ < 1 (necessary for the basic multiserver model)12.3. Parameter Drift and Model Validity
Validate arrival distributions, service demands, rejection rates, and latency distributions against observations. When workload regimes shift, old parameter estimates require revalidation.
12.4. Numerical and Algorithmic Example
In a stable M/M/1 model, mean latency diverges as utilization approaches one. Real systems may bound accumulation through finite buffers or timeouts, but then blocking and loss must be measured.
12.4.1. Illustrative Pseudocode
fit_arrival_and_service_distributions(data)
test_independence_and_stationarity()
simulate_with_multiple_random_seeds()
compare_means_percentiles_rejection_rates()
report_assumptions_and_uncertainty_intervals()12.5. Stability and Traffic Intensity
When long-run offered work persistently exceeds service capacity, queue population tends to grow. The simple M/M/c model requires ρ<1 for stability; general networks need analysis of each bottleneck, routing and workload dependencies. Even below average capacity, long arrival bursts can cause substantial temporary backlogs. Mean stability and bounded-delay guarantees are different concepts.
ρ = λ·E[S]/c; ρ < 1 (basic homogeneous server model)12.6. Backpressure, Admission and Load Shedding
Unbounded queues can accumulate memory during overload. Backpressure slows producers, admission control limits accepted work and load shedding intentionally rejects selected requests. These controls trade latency, loss, fairness and utilization differently. Ignoring retry behavior can create feedback where rejected jobs return as extra arrivals and worsen overload.
12.7. Calibration of Analytical Models Against Measurements
Analytical capacity formulas depend on stated distributional and independence assumptions. Real systems couple CPU, memory, I/O, locks and database resources. Calibration measures arrival and service distributions, p50/p95/p99 latency, queue population and throughput across load levels. Systematic residual patterns indicate model misspecification and should motivate model revision. Public academic reporting does not require exposing confidential production data.
12.8. Computational Example and Calculation
For M/M/1 with μ=10 jobs per second, λ=8 gives mean sojourn 0.5 seconds, while λ=9 gives 1 second. A 12.5-percent arrival increase doubles delay, illustrating nonlinear near-saturation behavior.
W₈=1/(10−8)=0.5 s; W₉=1/(10−9)=1 s12.9. Independent Verification and Implementation Constraints
Before transferring this prediction to a production system, validate service-time independence from load and the arrival-process distribution.
References
- Gross, Donald; Shortle, John F.; Thompson, James M.; Harris, Carl M. Fundamentals of Queueing Theory. 4th ed. Wiley, 2008.
- Ross, Sheldon M. Introduction to Probability Models. 13th ed. Academic Press, 2024.