Spaces:
Runtime error
Runtime error
File size: 9,328 Bytes
aaef24a f006b5f aaef24a 8eccedb aaef24a 8eccedb aaef24a f006b5f aaef24a f006b5f aaef24a | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 | # /// script
# requires-python = ">=3.13"
# dependencies = [
# "altair",
# "asimpy",
# "marimo",
# "polars==1.24.0",
# ]
# ///
import marimo
__generated_with = "0.22.4"
app = marimo.App(width="medium")
@app.cell(hide_code=True)
def _():
import marimo as mo
import math
import random
import altair as alt
import polars as pl
from asimpy import Environment, Process, Resource
return Environment, Process, Resource, alt, math, mo, pl, random
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
# Basic Ideas in Queueing Theory
## *Arrivals, Servers, and Utilization*
Three concepts underpin every queueing model. The first is *Poisson arrivals*: when customers arrive independently at a constant average rate $\lambda$, gaps between consecutive arrivals follow an *exponential* distribution with mean $1/\lambda$, and the count of arrivals in any window of width $t$ follows a Poisson distribution with mean $\lambda t$. These two descriptions are equivalent: if one holds, the other must too.
The second key idea is that the exponential distribution is *memoryless*. Knowing you have already waited $s$ units gives no information about when the next arrival will come. This property makes the math simple, but means that the exponential distribution isn't a good fit for scenarios where events happen in bursts. Some of the later tutorials will explore models that handle this.
The final concept is *server utilization*. A server completing work at rate $\mu$ has utilization $\rho = \lambda/\mu$, which is the long-run fraction of time it is busy. The system is stable only when $\rho < 1$. When $\rho \geq 1$, arrivals outpace service and the queue grows without bound. Even at exact balance ($\rho = 1$), randomness creates bursts that accumulate faster than the server recovers, so the expected queue length is infinite.
The code below uses a discrete event simulation package called [asimpy](https://asimpy.readthedocs.io/) to model a simple job queue. The technical term for this kind of system is an *[M/M/1 queue](https://en.wikipedia.org/wiki/M/M/1_queue)*: memoryless (Poisson) arrivals, memoryless service times, and one server.
""")
return
@app.cell(hide_code=True)
def _(mo):
sim_time_slider = mo.ui.slider(
start=0,
stop=100_000,
step=1_000,
value=20_000,
label="Simulation time"
)
arrival_rate_slider = mo.ui.slider(
start=1.0,
stop=5.0,
step=0.01,
value=2.0,
label="Arrival rate (λ)"
)
service_rate_slider = mo.ui.slider(
start=1.0,
stop=5.0,
step=0.01,
value=2.0,
label="Service rate (μ)"
)
window_slider = mo.ui.slider(
start=1.0,
stop=5.0,
step=0.01,
value=1.0,
label="Counting window"
)
seed_input = mo.ui.number(
value=192,
step=1,
label="Random seed"
)
run_button = mo.ui.run_button(label="Run simulation")
mo.vstack([sim_time_slider, arrival_rate_slider, service_rate_slider, window_slider, seed_input, run_button])
return (
arrival_rate_slider,
seed_input,
service_rate_slider,
sim_time_slider,
window_slider,
)
@app.cell
def _(
arrival_rate_slider,
random,
seed_input,
service_rate_slider,
sim_time_slider,
window_slider,
):
SIM_TIME = int(sim_time_slider.value)
ARRIVAL_RATE = float(arrival_rate_slider.value)
SERVICE_RATE = float(service_rate_slider.value)
WINDOW = float(window_slider.value)
SEED = int(seed_input.value)
random.seed(SEED)
return ARRIVAL_RATE, SERVICE_RATE, SIM_TIME, WINDOW
@app.cell
def _(Process, random):
class ArrivalSource(Process):
"""Generates arrivals at a Poisson rate and records inter-arrival gaps."""
def init(self, rate, gaps):
self.rate = rate
self.gaps = gaps
async def run(self):
while True:
gap = random.expovariate(self.rate)
await self.timeout(gap)
self.gaps.append(gap)
return (ArrivalSource,)
@app.cell
def _(Process, SERVICE_RATE, random):
class Customer(Process):
def init(self, server, total_service):
self.server = server
self.total_service = total_service
async def run(self):
async with self.server:
svc = random.expovariate(SERVICE_RATE)
await self.timeout(svc)
self.total_service[0] += svc
return (Customer,)
@app.cell
def _(Customer, Process, random):
class Arrivals(Process):
def init(self, rate, server, total_service):
self.rate = rate
self.server = server
self.total_service = total_service
async def run(self):
while True:
await self.timeout(random.expovariate(self.rate))
Customer(self._env, self.server, self.total_service)
return (Arrivals,)
@app.cell
def _(
ARRIVAL_RATE,
ArrivalSource,
Environment,
SIM_TIME,
WINDOW,
alt,
math,
pl,
):
def comparison():
gaps = []
env = Environment()
ArrivalSource(env, ARRIVAL_RATE, gaps)
env.run(until=SIM_TIME)
n = int(SIM_TIME / WINDOW)
counts = [0] * n
t = 0.0
for g in gaps:
t += g
w = int(t / WINDOW)
if w < n:
counts[w] += 1
freq = {}
for c in counts:
freq[c] = freq.get(c, 0) + 1
lam_w = ARRIVAL_RATE * WINDOW
return pl.DataFrame([
{
"k": k,
"observed": freq.get(k, 0) / n,
"theory": (lam_w**k) * math.exp(-lam_w) / math.factorial(k)
}
for k in range(max(counts) + 1)
])
_df = comparison().unpivot(
on=["observed", "theory"], index="k", variable_name="source", value_name="probability"
)
(
alt.Chart(_df).mark_bar(opacity=0.8)
.encode(
x=alt.X("k:O", title=f"Arrivals per window (width {WINDOW})"),
y=alt.Y("probability:Q", title="Probability"),
color=alt.Color("source:N", title="Series"),
xOffset="source:N",
tooltip=["k:O", "source:N", "probability:Q"],
)
.properties(title=f"Arrival Counts: Observed vs. Poisson(λ={ARRIVAL_RATE})")
)
return
@app.cell
def _(Arrivals, Environment, Resource, SERVICE_RATE, SIM_TIME):
def simulate(rho):
env = Environment()
server = Resource(env, capacity=1)
total_service = [0.0]
Arrivals(env, rho * SERVICE_RATE, server, total_service)
env.run(until=SIM_TIME)
busy = total_service[0] / SIM_TIME
return {"rho_target": rho, "rho_observed": round(busy, 4), "idle_frac": round(1.0 - busy, 4)}
return (simulate,)
@app.cell
def _(pl, simulate):
df_sim = pl.DataFrame([simulate(rho) for rho in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95]])
df_sim
return (df_sim,)
@app.cell
def _(alt, df_sim):
_df_plot = df_sim.unpivot(
on=["rho_observed", "idle_frac"],
index="rho_target",
variable_name="metric",
value_name="fraction",
)
chart = (
alt.Chart(_df_plot).mark_line(point=True)
.encode(
x=alt.X("rho_target:Q", title="Target utilization (ρ = λ/μ)"),
y=alt.Y("fraction:Q", title="Fraction of time"),
color=alt.Color("metric:N", title="Metric"),
tooltip=["rho_target:Q", "metric:N", "fraction:Q"],
)
.properties(title="Server Utilization: Busy and Idle Fractions vs. ρ")
)
chart
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
## Understanding the Math
### Why inter-arrival gaps are exponential
Divide $[0, t]$ into $n$ tiny slices of width $\Delta = t/n$. The probability of an arrival in each slice is $\approx \lambda\Delta$; slices are independent. The probability of no arrival across all $n$ slices is $(1 - \lambda t/n)^n \to e^{-\lambda t}$ as $n \to \infty$. This is $P(X > t)$ for the exponential distribution.
### Why mean equals standard deviation for the exponential
If $E[X] = 1/\lambda$, then $E[X^2] = 2/\lambda^2$, so $\text{Var}(X) = 2/\lambda^2 - 1/\lambda^2 = 1/\lambda^2$ and $\text{SD}(X) = 1/\lambda = E[X]$. Equal mean and standard deviation means roughly one-third of gaps are longer than the mean.
### Why the busy fraction equals $\rho$
Over a long run $T$, about $N \approx \lambda T$ customers are served, each occupying the server for mean $1/\mu$. Total service time $\approx N/\mu = \lambda T/\mu = \rho T$. Dividing by $T$ gives busy fraction $= \rho$.
### Why $\rho = 1$ is unstable
At exact balance, the queue length after each service completion performs a symmetric random walk on the non-negative integers. This random walk is *null recurrent*: it returns to zero but with infinite expected return time, so the mean queue length is infinite even when arrivals and service are perfectly matched on average.
""")
return
if __name__ == "__main__":
app.run()
|