marimo-team/marimo-learn
3
1# /// script2# requires-python = ">=3.13"3# dependencies = [4# "altair",5# "asimpy",6# "marimo",7# "polars==1.24.0",8# ]9# ///10 11import marimo12 13__generated_with = "0.20.4"14app = marimo.App(width="medium")15 16 17@app.cell(hide_code=True)18def _():19 import marimo as mo20 import random21 import statistics22 23 import altair as alt24 import polars as pl25 26 from asimpy import Environment, Process, Resource27 28 return Environment, Process, Resource, alt, mo, pl, random, statistics29 30 31@app.cell(hide_code=True)32def _(mo):33 mo.md(r"""34 # Queue Formation35 36 ## *Randomness Creates Waiting Even with Spare Capacity*37 38 We now combine arrivals (Poisson at rate $\lambda$) with a server (exponential service at rate $\mu$) into a complete queue. The system is stable because $\rho = \lambda/\mu < 1$, which means that on average, the server handles more work than arrives.39 40 Our first question is, how long is the queue? The surprising answer is that even when the server has plenty of spare capacity, customers wait. The mean number of customers in the system (both waiting and being served) is:41 42 $$L = \frac{\rho}{1 - \rho}$$43 44 The table below gives some representative values:45 46 | $\rho$ | $L$ |47 |:---:|:---:|48 | 0.1 | 0.11 |49 | 0.5 | 1.00 |50 | 0.8 | 4.00 |51 | 0.9 | 9.00 |52 53 When $\rho = 0.5$, half the server's capacity is idle, but there is on average one customer in the system at any moment. That customer either had to wait for a previous customer, or is currently being served. The queue is *never* consistently empty, even at moderate load.54 55 The formula also explains why simple queues are so sensitive to utilization: $L$ blows up as $\rho \to 1$. One way to think about this is that the denominator $(1 - \rho)$ is the spare capacity. As spare capacity vanishes, queue length increases.56 """)57 return58 59 60@app.cell61def _(mo):62 mo.md(r"""63 ## Why Queues Form at All64 65 With deterministic arrivals and service (every customer arrives exactly $1/\lambda$ apart and takes exactly $1/\mu$), a server with $\rho < 1$ would never form a queue: each customer would depart before the next arrived. Randomness changes this. Sometimes three customers arrive close together before the server finishes even one, so the server falls briefly behind. While it recovers, customers wait. These temporary pileups are unavoidable whenever inter-arrival or service times have any variance.66 67 The probability that exactly $n$ customers are in an M/M/1 system at steady state is:68 69 $$P(N = n) = (1 - \rho)\,\rho^n \qquad n = 0, 1, 2, \ldots$$70 71 This is a *geometric distribution* with success probability $1 - \rho$. The formula says the server is idle (i.e., n=0) with probability $1 - \rho$, which is consistent with the utilization result from the previous scenario. Each additional customer in the system is $\rho$ times less likely than the previous count.72 73 This formula $L = \rho/(1-\rho)$ is the foundation of the later M/M/1 nonlinearity scenario, which shows the practical consequences of the $(1-\rho)$ denominator. Every queue-length formula in queueing theory has a similar structure: a traffic factor $\rho$ divided by a spare-capacity factor $(1 - \rho)$, possibly multiplied by a variability correction.74 """)75 return76 77 78@app.cell79def _(mo):80 mo.md(r"""81 ## Implementation82 83 A `Customer` process increments a shared `in_system` counter on arrival and decrements it on departure. A `Monitor` process samples `in_system[0]` every `SAMPLE_INTERVAL` time units. After the simulation, the mean of the samples estimates $L$. The theoretical value $\rho/(1-\rho)$ is computed and compared. By the law of large numbers, this converges to the true steady-state mean as the simulation time grows.84 85 The simulation sweeps $\rho$ from 0.1 to 0.9, confirming the formula at each load level.86 """)87 return88 89 90@app.cell(hide_code=True)91def _(mo):92 sim_time_slider = mo.ui.slider(93 start=0,94 stop=100_000,95 step=1_000,96 value=20_000,97 label="Simulation time",98 )99 100 service_rate_slider = mo.ui.slider(101 start=1.0,102 stop=5.0,103 step=0.01,104 value=2.0,105 label="Service rate",106 )107 108 sample_interval_slider = mo.ui.slider(109 start=1.0,110 stop=5.0,111 step=1.0,112 value=1.0,113 label="Sample interval",114 )115 116 seed_input = mo.ui.number(117 value=192,118 step=1,119 label="Random seed",120 )121 122 run_button = mo.ui.run_button(label="Run simulation")123 124 mo.vstack([125 sim_time_slider,126 service_rate_slider,127 sample_interval_slider,128 seed_input,129 run_button,130 ])131 return (132 sample_interval_slider,133 seed_input,134 service_rate_slider,135 sim_time_slider,136 )137 138 139@app.cell140def _(141 sample_interval_slider,142 seed_input,143 service_rate_slider,144 sim_time_slider,145):146 SIM_TIME = int(sim_time_slider.value)147 SERVICE_RATE = float(service_rate_slider.value)148 SAMPLE_INTERVAL = float(sample_interval_slider.value)149 SEED = int(seed_input.value)150 return SAMPLE_INTERVAL, SEED, SERVICE_RATE, SIM_TIME151 152 153@app.cell154def _(Process, SERVICE_RATE, random):155 class Customer(Process):156 def init(self, server, in_system):157 self.server = server158 self.in_system = in_system159 160 async def run(self):161 self.in_system[0] += 1162 async with self.server:163 await self.timeout(random.expovariate(SERVICE_RATE))164 self.in_system[0] -= 1165 166 return (Customer,)167 168 169@app.cell170def _(Customer, Process, random):171 class Arrivals(Process):172 def init(self, rate, server, in_system):173 self.rate = rate174 self.server = server175 self.in_system = in_system176 177 async def run(self):178 while True:179 await self.timeout(random.expovariate(self.rate))180 Customer(self._env, self.server, self.in_system)181 182 return (Arrivals,)183 184 185@app.cell186def _(Process, SAMPLE_INTERVAL):187 class Monitor(Process):188 """Samples total customers in system at regular intervals."""189 190 def init(self, in_system, samples):191 self.in_system = in_system192 self.samples = samples193 194 async def run(self):195 while True:196 self.samples.append(self.in_system[0])197 await self.timeout(SAMPLE_INTERVAL)198 199 return (Monitor,)200 201 202@app.cell203def _(204 Arrivals,205 Environment,206 Monitor,207 Resource,208 SERVICE_RATE,209 SIM_TIME,210 statistics,211):212 def simulate(rho):213 arrival_rate = rho * SERVICE_RATE214 env = Environment()215 server = Resource(env, capacity=1)216 in_system = [0]217 samples = []218 Arrivals(env, arrival_rate, server, in_system)219 Monitor(env, in_system, samples)220 env.run(until=SIM_TIME)221 sim_L = statistics.mean(samples)222 theory_L = rho / (1.0 - rho)223 return {224 "rho": rho,225 "sim_L": round(sim_L, 4),226 "theory_L": round(theory_L, 4),227 "error_pct": round(100.0 * (sim_L - theory_L) / theory_L, 2),228 }229 230 return (simulate,)231 232 233@app.cell234def _(SEED, pl, random, simulate):235 def sweep():236 rows = [simulate(rho) for rho in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]]237 return pl.DataFrame(rows)238 239 random.seed(SEED)240 df = sweep()241 df242 return (df,)243 244 245@app.cell246def _(alt, df):247 df_plot = df.unpivot(248 on=["sim_L", "theory_L"],249 index="rho",250 variable_name="source",251 value_name="L",252 )253 chart = (254 alt.Chart(df_plot)255 .mark_line(point=True)256 .encode(257 x=alt.X("rho:Q", title="Utilization (ρ)"),258 y=alt.Y("L:Q", title="Mean customers in system (L)"),259 color=alt.Color("source:N", title="Source"),260 tooltip=["rho:Q", "source:N", "L:Q"],261 )262 .properties(title="Queue Formation: Simulated vs. Theoretical L = ρ/(1−ρ)")263 )264 chart265 return266 267 268@app.cell269def _(mo):270 mo.md(r"""271 ## Understanding the Math272 273 ### Why is the queue length geometric?274 275 The M/M/1 queue can be analyzed as a random walk on the non-negative integers. When the server is busy, the queue grows by 1 with each arrival (with rate $\lambda$) and shrinks by 1 with each service completion (with rate $\mu$). The ratio $\lambda/\mu = \rho$ is the probability that the queue grows rather than shrinks at the next event. Under steady state, the probability of being at level $n$ is proportional to $\rho^n$ — because reaching level $n$ requires $n$ consecutive "up" steps. Normalizing so the probabilities sum to 1 gives $(1-\rho)\rho^n$.276 277 ### Deriving $L = \rho/(1-\rho)$ from the geometric distribution278 279 Given $P(N = n) = (1 - \rho)\rho^n$, the mean is:280 281 $$L = E[N] = \sum_{n=0}^{\infty} n \cdot (1-\rho)\rho^n = (1-\rho) \sum_{n=0}^{\infty} n\rho^n$$282 283 A result from basic calculus is that the geometric series $\sum_{n=0}^{\infty} \rho^n = 1/(1-\rho)$. Differentiating both sides with respect to $\rho$:284 285 $$\sum_{n=0}^{\infty} n\rho^{n-1} = \frac{1}{(1-\rho)^2}$$286 287 Multiply both sides by $\rho$:288 289 $$\sum_{n=0}^{\infty} n\rho^n = \frac{\rho}{(1-\rho)^2}$$290 291 Substituting back:292 293 $$L = (1-\rho) \cdot \frac{\rho}{(1-\rho)^2} = \frac{\rho}{1-\rho}$$294 295 ### Checking the formula at the boundaries296 297 When $\rho \to 0$: there are almost no arrivals, so $L \to 0$, i.e., the server is nearly always idle.298 299 When $\rho \to 1$: the spare capacity is $(1-\rho) \to 0$, so $L \to \infty$, i.e., the queue grows without bound.300 301 Both limits match physical intuition. Note that the formula is exact (not an approximation) for an M/M/1 queues in steady state.302 """)303 return304 305 306if __name__ == "__main__":307 app.run()308