Discrete-Event Simulation by Building

Week 2 — Replications and confidence intervals

Quisim

Week 2: How much can we trust one run?

This week

Goal

Find out how good a simulation number is: run it again, run it many times, and put an interval around the answer.

• Rerun M/M/k: same model, different answers
• Build an M/M/1 with seed knobs and make replications by hand
• Let the Statistics example do it for you: mean ± radius
• Test the simulation against the Erlang C formula
• The price of precision: how many runs are enough?

How to follow along in the free demo

• Show the editor pane (pen icon, Ctrl+E): below it are the Examples and Blocks tabs
• Examples tab: click M/M/k, Statistics or Erlang C to load it
• Blocks tab: on an empty canvas, click a block to add it; drag from socket to socket to draw a pipe (green) or a wire (silver)
• Type new values into the NumberInput fields; press Enter to keep them
• Start clean (Ctrl+Enter) runs from time 0; Animate (film icon) to watch, off for fast runs
• Monitor tab: the event list, while animating or paused

The demo cannot open or paste code. Full-version users can copy each slide's Model code instead.

Recap: the M/M/1 example from week 1

Model codeOpen in Quisim

NumberInput { id: number0; initValue:125; label:"1/λ"; lower:1 }
NumberInput { id: number1; initValue:100; label:"1/µ"; lower:1 }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
RandomExponential { id:rnd0 }
RandomExponential { id:rnd1 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration }
Arrive { id:arrive0 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Perform { id:perform0 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0 }
Wire { from: queue0.measureLoad; to: chart0.parameters[0] }
Statistics { id: statistics0; height:48; width:128; label:"Sojourn\ntime" }
Wire { from: leave0.measureSojourn; to: statistics0.parameters[0] }
Chart { id: chart0; height:66; width:384 }

$\rho = \lambda/\mu = 100/125 = 0.8$, and theory says $W = \dfrac{1}{\mu-\lambda} = 500$ ms. Week 1's runs rarely showed 500. Today we find out why.

Same model, different answers

Reopen the M/M/k example

Model codeOpen in Quisim

NumberInput { id: number0; initValue:55; label:"1/λ"; lower:1 }
NumberInput { id: number1; initValue:500; label:"1/µ"; lower:1 }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
RandomExponential { id:rnd0 }
RandomExponential { id:rnd1 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration }
Arrive { id:arrive0 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0 }
Pipe { from:queue0.outlet; to:perform0.inlet }
NumberInput { id: number2; initValue:10; label:"Capacity"; lower:1 }
Wire { from: number2.measureValue; to: perform0.parameterCapacity }
Perform { id:perform0 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0 }
Wire { from: queue0.measureLoad; to: chart0.parameters[0] }
Wire { from: perform0.measureLoad; to: chart0.parameters[1] }
Statistics { id: statistics0; height:48; width:128; label:"Sojourn\ntime" }
Wire { from: leave0.measureSojourn; to: statistics0.parameters[0] }
Chart { id: chart0; numParameters:2; height:66; width:384 }

10 servers, $\rho = \dfrac{500}{55 \cdot 10} \approx 0.91$. Examples tab → M/M/k → Start clean. Note the Sojourn time.

Predict: press Start clean again

Question

Nothing in the model changed. Will the mean sojourn time be the same as before?

Run it: a different number every time. In our test runs a single 10 s run gave anything from about 450 to over 1 100 ms, on average about 650; Erlang C says $W \approx 883$ ms.

• The random blocks have no fixed seed: each Start clean draws a new history from a clock-based seed
• The runs also start empty, so early customers never wait: most runs land below 883

Your turn: three runs of M/M/k

Is the 11th server worth it? (5 min)

In the M/M/k example, press Start clean three times with Capacity 10, then three times with Capacity 11. Record the Sojourn time.

Capacity	run 1	run 2	run 3	lowest	highest
10
11

Do the two ranges overlap? Could you tell your manager "11 servers cut the time in system by X ms"?

Key insight: a run is a sample

One simulation run is one sample of a random experiment, like watching a call centre on one single Monday.

• Nobody sizes a call centre from one Monday
• To learn about the mean, we need many samples and a little statistics

Be the Steering block: replications by hand

Build it: M/M/1 (1/3)

Model codeOpen in Quisim

NumberInput { id: number0; initValue:125; label:"1/λ" }
NumberInput { id: number1; initValue:100; label:"1/µ" }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
RandomExponential { id:rnd0 }
RandomExponential { id:rnd1 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration }
Arrive { id:arrive0 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Perform { id:perform0 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0 }

New model (Ctrl+N). Pipes: Arrive O → Queue I → Perform I, Perform O → Leave I. Wires: RandomExponential V → Arrive D, a second RandomExponential V → Perform D, NumberInput N → each µ socket. Type 125 and 100.

Build it: add a Statistics and a Chart (2/3)

Model codeOpen in Quisim

NumberInput { id: number0; initValue:125; label:"1/λ" }
NumberInput { id: number1; initValue:100; label:"1/µ" }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
RandomExponential { id:rnd0 }
RandomExponential { id:rnd1 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration }
Arrive { id:arrive0 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Perform { id:perform0 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0 }
Wire { from: queue0.measureLoad; to: chart0.parameters[0] }
Statistics { id: statistics0 }
Wire { from: leave0.measureSojourn; to: statistics0.parameters[0] }
Chart { id: chart0 }

Wires: Leave LS (sojourn time) → Statistics P[0]; Queue L (queue length) → Chart P[0]. Start clean: your canvas now does what the M/M/1 example does.

Build it: two seed knobs (3/3)

Model codeOpen in Quisim

NumberInput { id: number0; initValue:125; label:"1/λ" }
NumberInput { id: number1; initValue:100; label:"1/µ" }
NumberInput { id: number2; initValue:1; label:"seed A" }
NumberInput { id: number3; initValue:2; label:"seed S" }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
Wire { from: number2.measureValue; to: rnd0.parameterSeed }
Wire { from: number3.measureValue; to: rnd1.parameterSeed }
RandomExponential { id:rnd0 }
RandomExponential { id:rnd1 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration }
Arrive { id:arrive0 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Perform { id:perform0 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0 }
Wire { from: queue0.measureLoad; to: chart0.parameters[0] }
Statistics { id: statistics0 }
Wire { from: leave0.measureSojourn; to: statistics0.parameters[0] }
Chart { id: chart0 }

Two more NumberInputs: N → S (seed) of the arrival RandomExponential, N → S of the service one. Type 1 and 2. Seed 0 means "no fixed seed": a new history every run. Demo blocks have no labels: remember which knob is which.

Predict: same seeds, run again

Question

Seeds 1 and 2, Start clean twice. Then seeds 3 and 4. What do you expect each time?

Run it: with seeds 1 and 2, both runs agree to the last digit, Chart included. Seeds 3 and 4 give a new history and a new mean.

• Same seed, same "random" numbers: results are reproducible, a bug can be found again
• Use different seeds for the two streams: the same seed twice would make every service time exactly 0.8 × an interarrival time (try 1 and 1: nobody ever waits, sojourn about 80)

Your turn: five replications by hand

You are the Steering block (10 min)

Run your model with five seed pairs and fill in the mean Sojourn per run.

seeds	1, 2	3, 4	5, 6	7, 8	9, 10
mean sojourn $\hat\theta_i$

Compute the average $\hat\Theta$, the standard deviation $S$, and $\delta = 2.78 \cdot S/\sqrt{5}$. Is 500 inside $\hat\Theta \pm \delta$?

Expect the average of a single run anywhere from about 150 to 900 ms, most below 500, averaging about 400: each run starts empty.

From five runs to one interval

Your table holds $R = 5$ independent estimates $\hat\theta_i$. Their average is roughly normal (Central Limit Theorem), so

$$\hat\Theta \pm \delta, \qquad \delta = t_{R-1,\,0.975}\cdot\frac{S_R}{\sqrt{R}}, \qquad S_R^2 = \frac{1}{R-1}\sum_{i=1}^{R}\big(\hat\theta_i-\hat\Theta\big)^2$$

• $S_R$: how much your five runs disagree
• $t_{R-1,\,0.975}$: wide for small $R$ (2.78 at $R=5$, 2.09 at $R=20$, → 1.96)
• $R = 1$ gives no interval: one number has no spread

Replications in Quisim: the Statistics example

The Statistics example

Model codeOpen in Quisim

NumberInput { id: number1; initValue:100; label:"1/µ\nservice time"; lower:1; height:58; width:96 }
NumberInput { id: number0; initValue:125; label:"1/λ\ninterarrival time"; lower:1; height:58; width:112 }
RandomExponential { id:rnd0 }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from: number2.measureValue; to: perform0.parameterCapacity }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration }
NumberInput { id: numberinput1; initValue:10000; label:"Duration"; width:112; lower:1; upper:1000000 }
NumberInput { id: numberinput0; initValue:20; label:"Replications"; width:112; lower:1; upper:100 }
NumberInput { id: number2; initValue:1; label:"k servers"; lower:1 }
Wire { from: numberinput1.measureValue; to: steering0.parameterDuration }
Steering { id: steering0; duration:50000 }
Wire { from: numberinput0.measureValue; to: steering0.parameterReplications }
RandomExponential { id:rnd1 }
Arrive { id:arrive0 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Perform { id:perform0 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0 }
Statistics { id: statistics0; timeAverage:true; numParameters:2; warmupWidth:5000; label:"Queue length\nService load"; height:192; width:152; x:646; y:245 }
Wire { from: queue0.measureLoad; to: statistics0.parameters[0] }
Wire { from: perform0.measureLoad; to: statistics0.parameters[1] }
Statistics { id: statistics2; warmupWidth:5000; label:"Sojourn time"; height:192; width:152; x:646; y:12 }
Statistics { id: statistics1; numParameters:2; warmupWidth:5000; label:"Waiting time\nService time"; height:192; width:152; x:470; y:245 }
Wire { from: queue0.measureWait; to: statistics1.parameters[0] }
Wire { from: rnd1.measureRandomVariate; to: statistics1.parameters[1] }
Wire { from: leave0.measureSojourn; to: statistics2.parameters[0] }

Your M/M/1 plus a Steering block: Replications $R = 20$ runs of Duration 10 000 ms, each with fresh random numbers. Each Statistics block drops the first 5 000 ms (warm-up) and shows mean ± radius $\delta$ in its last row.

Predict: Replications = 1

Question

Set Replications to 1 and Start clean. What will the Sojourn time block show in its last row?

Run it: one replication

Model codeOpen in Quisim

NumberInput { id: number1; initValue:100; label:"1/µ\nservice time"; lower:1; height:58; width:96 }
NumberInput { id: number0; initValue:125; label:"1/λ\ninterarrival time"; lower:1; height:58; width:112 }
RandomExponential { id:rnd0 }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from: number2.measureValue; to: perform0.parameterCapacity }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration }
NumberInput { id: numberinput1; initValue:10000; label:"Duration"; width:112; lower:1; upper:1000000 }
NumberInput { id: numberinput0; initValue:1; label:"Replications"; width:112; lower:1; upper:100 }
NumberInput { id: number2; initValue:1; label:"k servers"; lower:1 }
Wire { from: numberinput1.measureValue; to: steering0.parameterDuration }
Steering { id: steering0; duration:50000 }
Wire { from: numberinput0.measureValue; to: steering0.parameterReplications }
RandomExponential { id:rnd1 }
Arrive { id:arrive0 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Perform { id:perform0 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0 }
Statistics { id: statistics0; timeAverage:true; numParameters:2; warmupWidth:5000; label:"Queue length\nService load"; height:192; width:152; x:646; y:245 }
Wire { from: queue0.measureLoad; to: statistics0.parameters[0] }
Wire { from: perform0.measureLoad; to: statistics0.parameters[1] }
Statistics { id: statistics2; warmupWidth:5000; label:"Sojourn time"; height:192; width:152; x:646; y:12 }
Statistics { id: statistics1; numParameters:2; warmupWidth:5000; label:"Waiting time\nService time"; height:192; width:152; x:470; y:245 }
Wire { from: queue0.measureWait; to: statistics1.parameters[0] }
Wire { from: rnd1.measureRandomVariate; to: statistics1.parameters[1] }
Wire { from: leave0.measureSojourn; to: statistics2.parameters[0] }

Run it: a mean, but no ± radius. One run has no spread to estimate. Back to 20 and the radius returns.

Your turn: read all four statistics

Check the model against what you know (5 min)

Statistics example, default values ($R = 20$). Record mean ± $\delta$:

Sojourn $W$	Waiting $W_q$	Service time	Queue length $L_q$	Service load

1. Is Service time close to $1/\mu = 100$? Is Service load close to $\rho = 0.8$?
2. Little's law says $L_q = \lambda W_q = W_q / 125$. Does it hold for your numbers?

Key insight: replications do your table for you

The Steering block runs the model $R$ times with fresh seeds; the Statistics block turns the $R$ means into $\hat\Theta \pm \delta$, exactly the calculation you did by hand.

Simulation meets theory: Erlang C

The Erlang C example

Model codeOpen in Quisim

NumberInput { id: numberinput1; initValue:10000; label:"Duration"; width:96; lower:1; upper:1000000; x:22; y:464 }
NumberInput { id: numberinput0; initValue:20; label:"Replications"; width:96; lower:1; upper:100; x:22; y:524 }
Wire { from: numberinput0.measureValue; to: steering0.parameterReplications }
Wire { from: numberinput1.measureValue; to: steering0.parameterDuration }
Steering { id: steering0; replications:20; duration:100000; x:161; y:464 }

NumberInput { id: number0; initValue:125; label:"1/λ\ninterarrival\ntime"; lower:1; x:22; y:48; height:72; width:96 }
NumberInput { id: number1; initValue:100; label:"1/µ\nservice time"; x:22; y:164; lower:1; height:58; width:96 }
NumberInput { id: number2; initValue:1; x:22; y:284; label:"k servers"; lower:1; width:96 }
Wire { from: number0.measureValue; to: function0.parameters[0] }
Wire { from: number1.measureValue; to: function0.parameters[1] }
Wire { from: number2.measureValue; to: function0.parameters[2] }
Function { id: function0; code:"const a = P[1] / P[0], ρ = a / P[2];\nif (ρ >= 0 && ρ < 1) { // Jagerman method\n  let t = 1, s = 1, n = 1;\n  while (n < P[2]) { t *= a / n; s += t; n += 1; }\n  t *= a / n;\n  const Pidle = 1 / (s + t / (1 - ρ)), Pwait = t / (1 - ρ) * Pidle;\n  const Lq = Pwait * ρ / (1 - ρ), Wq = Lq * P[0];\n  M = [ρ, Pwait, Lq, Lq + a, Wq, Wq + P[1], Pidle];\n}\n\n\n\n"; numMeasures:7; numParameters:3; label:"Erlang C"; height:48; width:96; x:162; y:284 }
Number { id: num0; label:"ρ utilization"; width: 156; x:313; y:164 }
Wire { from: function0.measures[0]; to: num0.parameterValue }
Number { id: num1; label:"P(wait) Erlang C"; width: 156; x:313; y:224 }
Wire { from: function0.measures[1]; to: num1.parameterValue }

Number { id: num2; label:"Lq mean queue length"; width: 156; x:313; y:284 }
Wire { from: function0.measures[2]; to: num2.parameterValue }

Number { id: num3; label:"L mean in system"; width: 156; x:313; y:344 }
Wire { from: function0.measures[3]; to: num3.parameterValue }

Number { id: num4; label:"Wq mean wait time"; width: 156; x:313; y:404 }
Wire { from: function0.measures[4]; to: num4.parameterValue }

Number { id: num5; label:"W mean sojourn time"; width: 156; x:313; y:464 }
Wire { from: function0.measures[5]; to: num5.parameterValue }

Number { id: num6; label:"P(idle) idle probability"; width: 156; x:313; y:524 }
Wire { from: function0.measures[6]; to: num6.parameterValue }

Wire { from: number0.measureValue; to: rnd0.parameterMean }
RandomExponential { id:rnd0; x:162; y:48; width:96 }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
RandomExponential { id:rnd1; x:162; y:164; width:96 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration; bend:-90; incline:100; debend:-100; decline:100 }
Arrive { id:arrive0; x:312; y:48 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; x:421; y:48 }
Wire { from: queue0.measureLoad; to: chart0.parameters[0] }
Chart { id: chart0; x:502; y:387; label:"Queue length"; height:185; width:421 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Wire { from: number2.measureValue; to: perform0.parameterCapacity; bend:-20; incline:20; decline:-10 }
Perform { id:perform0; x:529; y:48 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; x:631; y:48 }
Statistics { id: statistics1; numParameters:2; warmupWidth:5000; x:503; y:166; label:"Wq Waiting time\n1/µ Service time" }
Wire { from: queue0.measureWait; to: statistics1.parameters[0] }
Wire { from: rnd1.measureRandomVariate; to: statistics1.parameters[1] }
Statistics { id: statistics0; timeAverage:true; numParameters:2; warmupWidth:5000; x:719; y:166; label:"Lq Queue length\nLs=L-Lq Service load" }
Wire { from: queue0.measureLoad; to: statistics0.parameters[0] }
Wire { from: perform0.measureLoad; to: statistics0.parameters[1] }
Statistics { id: statistics2; warmupWidth:5000; x:719; y:28; label:"W Sojourn time"; height:96 }
Wire { from: leave0.measureSojourn; to: statistics2.parameters[0] }

The same simulation, plus a Function block that computes the Erlang C formula for the same inputs; the Number blocks in the middle show its results.

Erlang C in one slide

For M/M/$k$ with offered load $a = \lambda/\mu$ and $\rho = a/k < 1$:

$$P(\text{wait}) = C(k,a) = \frac{\dfrac{a^k}{k!}\,\dfrac{1}{1-\rho}}{\displaystyle\sum_{n=0}^{k-1}\frac{a^n}{n!} + \frac{a^k}{k!}\,\frac{1}{1-\rho}}, \qquad W_q = \frac{C(k,a)}{k\mu-\lambda}, \qquad W = W_q + \frac{1}{\mu}$$

• Assumes Poisson arrivals, exponential service, $k$ servers, one FIFO queue, steady state
• Example defaults ($k = 1$): $\rho = 0.8$, $P(\text{wait}) = 0.8$, $W_q = 400$, $W = 500$ ms: read them off the Number blocks

Predict: is 500 inside the interval?

Question

Erlang C example, default values, Start clean. The formula says $W = 500$. Will 500 lie inside the W Sojourn time interval $\hat\Theta \pm \delta$?

Run it: often not. A typical result is about 420 ± 130, and in our tests about 1 run in 3 missed 500.

• A 95% interval should miss only 1 time in 20. Something is off: keep that thought

Your turn: the class experiment

Everyone runs the same model (5 min)

Erlang C example, defaults ($R = 20$). Every Start clean takes a new seed from the clock, so everyone gets a different interval. Do three runs:

run	mean $\hat\Theta$	radius $\delta$	500 inside?
1
2
3

Show of hands: how many of the class's intervals missed 500? Expected for 95%: about 1 in 20.

What "95%" means

The 95% describes the procedure

If the model matches the formula, about 95 of 100 such intervals cover the true mean. Your one interval either covers it or not.

• It does not say 95% of customers wait in this range: single waits vary far more
• So if a quarter of the class misses, either the class was unlucky, or the model does not measure what the formula describes

Predict: will more replications fix it?

Question

Set Replications to 100. Will the interval now cover 500 more often, or less often?

Run it: 100 replications

Model codeOpen in Quisim

NumberInput { id: numberinput1; initValue:10000; label:"Duration"; width:96; lower:1; upper:1000000; x:22; y:464 }
NumberInput { id: numberinput0; initValue:100; label:"Replications"; width:96; lower:1; upper:100; x:22; y:524 }
Wire { from: numberinput0.measureValue; to: steering0.parameterReplications }
Wire { from: numberinput1.measureValue; to: steering0.parameterDuration }
Steering { id: steering0; replications:20; duration:100000; x:161; y:464 }

NumberInput { id: number0; initValue:125; label:"1/λ\ninterarrival\ntime"; lower:1; x:22; y:48; height:72; width:96 }
NumberInput { id: number1; initValue:100; label:"1/µ\nservice time"; x:22; y:164; lower:1; height:58; width:96 }
NumberInput { id: number2; initValue:1; x:22; y:284; label:"k servers"; lower:1; width:96 }
Wire { from: number0.measureValue; to: function0.parameters[0] }
Wire { from: number1.measureValue; to: function0.parameters[1] }
Wire { from: number2.measureValue; to: function0.parameters[2] }
Function { id: function0; code:"const a = P[1] / P[0], ρ = a / P[2];\nif (ρ >= 0 && ρ < 1) { // Jagerman method\n  let t = 1, s = 1, n = 1;\n  while (n < P[2]) { t *= a / n; s += t; n += 1; }\n  t *= a / n;\n  const Pidle = 1 / (s + t / (1 - ρ)), Pwait = t / (1 - ρ) * Pidle;\n  const Lq = Pwait * ρ / (1 - ρ), Wq = Lq * P[0];\n  M = [ρ, Pwait, Lq, Lq + a, Wq, Wq + P[1], Pidle];\n}\n\n\n\n"; numMeasures:7; numParameters:3; label:"Erlang C"; height:48; width:96; x:162; y:284 }
Number { id: num0; label:"ρ utilization"; width: 156; x:313; y:164 }
Wire { from: function0.measures[0]; to: num0.parameterValue }
Number { id: num1; label:"P(wait) Erlang C"; width: 156; x:313; y:224 }
Wire { from: function0.measures[1]; to: num1.parameterValue }

Number { id: num2; label:"Lq mean queue length"; width: 156; x:313; y:284 }
Wire { from: function0.measures[2]; to: num2.parameterValue }

Number { id: num3; label:"L mean in system"; width: 156; x:313; y:344 }
Wire { from: function0.measures[3]; to: num3.parameterValue }

Number { id: num4; label:"Wq mean wait time"; width: 156; x:313; y:404 }
Wire { from: function0.measures[4]; to: num4.parameterValue }

Number { id: num5; label:"W mean sojourn time"; width: 156; x:313; y:464 }
Wire { from: function0.measures[5]; to: num5.parameterValue }

Number { id: num6; label:"P(idle) idle probability"; width: 156; x:313; y:524 }
Wire { from: function0.measures[6]; to: num6.parameterValue }

Wire { from: number0.measureValue; to: rnd0.parameterMean }
RandomExponential { id:rnd0; x:162; y:48; width:96 }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
RandomExponential { id:rnd1; x:162; y:164; width:96 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration; bend:-90; incline:100; debend:-100; decline:100 }
Arrive { id:arrive0; x:312; y:48 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; x:421; y:48 }
Wire { from: queue0.measureLoad; to: chart0.parameters[0] }
Chart { id: chart0; x:502; y:387; label:"Queue length"; height:185; width:421 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Wire { from: number2.measureValue; to: perform0.parameterCapacity; bend:-20; incline:20; decline:-10 }
Perform { id:perform0; x:529; y:48 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; x:631; y:48 }
Statistics { id: statistics1; numParameters:2; warmupWidth:5000; x:503; y:166; label:"Wq Waiting time\n1/µ Service time" }
Wire { from: queue0.measureWait; to: statistics1.parameters[0] }
Wire { from: rnd1.measureRandomVariate; to: statistics1.parameters[1] }
Statistics { id: statistics0; timeAverage:true; numParameters:2; warmupWidth:5000; x:719; y:166; label:"Lq Queue length\nLs=L-Lq Service load" }
Wire { from: queue0.measureLoad; to: statistics0.parameters[0] }
Wire { from: perform0.measureLoad; to: statistics0.parameters[1] }
Statistics { id: statistics2; warmupWidth:5000; x:719; y:28; label:"W Sojourn time"; height:96 }
Wire { from: leave0.measureSojourn; to: statistics2.parameters[0] }

Run it: the interval shrinks to about ± 60, but stays centred near 420: now it misses 500 more often than not, in our tests about 2 runs in 3.

Why the simulation sits below 500

• Each run starts with an empty queue, and only 10 000 ms long
• The warm-up cut drops the first 5 000 ms, but at $\rho = 0.8$ the queue needs many thousands of ms to "forget" the empty start
• So every replication measures a slightly too quiet system: a bias, the same in every run

More $R$ makes it worse

Replications shrink the noise, so the interval tightens around the wrong value.

Your turn: longer runs

Duration 10 000 against 100 000 (10 min)

Erlang C example, $R = 20$, Animate off. Two runs per setting:

Duration	mean $\hat\Theta$	radius $\delta$	500 inside?
10 000
10 000
100 000
100 000

What moved: the mean, the radius, or both?

Run it: a long run lands on theory

Model codeOpen in Quisim

NumberInput { id: numberinput1; initValue:100000; label:"Duration"; width:96; lower:1; upper:1000000; x:22; y:464 }
NumberInput { id: numberinput0; initValue:20; label:"Replications"; width:96; lower:1; upper:100; x:22; y:524 }
Wire { from: numberinput0.measureValue; to: steering0.parameterReplications }
Wire { from: numberinput1.measureValue; to: steering0.parameterDuration }
Steering { id: steering0; replications:20; duration:100000; x:161; y:464 }

NumberInput { id: number0; initValue:125; label:"1/λ\ninterarrival\ntime"; lower:1; x:22; y:48; height:72; width:96 }
NumberInput { id: number1; initValue:100; label:"1/µ\nservice time"; x:22; y:164; lower:1; height:58; width:96 }
NumberInput { id: number2; initValue:1; x:22; y:284; label:"k servers"; lower:1; width:96 }
Wire { from: number0.measureValue; to: function0.parameters[0] }
Wire { from: number1.measureValue; to: function0.parameters[1] }
Wire { from: number2.measureValue; to: function0.parameters[2] }
Function { id: function0; code:"const a = P[1] / P[0], ρ = a / P[2];\nif (ρ >= 0 && ρ < 1) { // Jagerman method\n  let t = 1, s = 1, n = 1;\n  while (n < P[2]) { t *= a / n; s += t; n += 1; }\n  t *= a / n;\n  const Pidle = 1 / (s + t / (1 - ρ)), Pwait = t / (1 - ρ) * Pidle;\n  const Lq = Pwait * ρ / (1 - ρ), Wq = Lq * P[0];\n  M = [ρ, Pwait, Lq, Lq + a, Wq, Wq + P[1], Pidle];\n}\n\n\n\n"; numMeasures:7; numParameters:3; label:"Erlang C"; height:48; width:96; x:162; y:284 }
Number { id: num0; label:"ρ utilization"; width: 156; x:313; y:164 }
Wire { from: function0.measures[0]; to: num0.parameterValue }
Number { id: num1; label:"P(wait) Erlang C"; width: 156; x:313; y:224 }
Wire { from: function0.measures[1]; to: num1.parameterValue }

Number { id: num2; label:"Lq mean queue length"; width: 156; x:313; y:284 }
Wire { from: function0.measures[2]; to: num2.parameterValue }

Number { id: num3; label:"L mean in system"; width: 156; x:313; y:344 }
Wire { from: function0.measures[3]; to: num3.parameterValue }

Number { id: num4; label:"Wq mean wait time"; width: 156; x:313; y:404 }
Wire { from: function0.measures[4]; to: num4.parameterValue }

Number { id: num5; label:"W mean sojourn time"; width: 156; x:313; y:464 }
Wire { from: function0.measures[5]; to: num5.parameterValue }

Number { id: num6; label:"P(idle) idle probability"; width: 156; x:313; y:524 }
Wire { from: function0.measures[6]; to: num6.parameterValue }

Wire { from: number0.measureValue; to: rnd0.parameterMean }
RandomExponential { id:rnd0; x:162; y:48; width:96 }
Wire { from: number1.measureValue; to: rnd1.parameterMean }
RandomExponential { id:rnd1; x:162; y:164; width:96 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration; bend:-90; incline:100; debend:-100; decline:100 }
Arrive { id:arrive0; x:312; y:48 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; x:421; y:48 }
Wire { from: queue0.measureLoad; to: chart0.parameters[0] }
Chart { id: chart0; x:502; y:387; label:"Queue length"; height:185; width:421 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Wire { from: number2.measureValue; to: perform0.parameterCapacity; bend:-20; incline:20; decline:-10 }
Perform { id:perform0; x:529; y:48 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; x:631; y:48 }
Statistics { id: statistics1; numParameters:2; warmupWidth:5000; x:503; y:166; label:"Wq Waiting time\n1/µ Service time" }
Wire { from: queue0.measureWait; to: statistics1.parameters[0] }
Wire { from: rnd1.measureRandomVariate; to: statistics1.parameters[1] }
Statistics { id: statistics0; timeAverage:true; numParameters:2; warmupWidth:5000; x:719; y:166; label:"Lq Queue length\nLs=L-Lq Service load" }
Wire { from: queue0.measureLoad; to: statistics0.parameters[0] }
Wire { from: perform0.measureLoad; to: statistics0.parameters[1] }
Statistics { id: statistics2; warmupWidth:5000; x:719; y:28; label:"W Sojourn time"; height:96 }
Wire { from: leave0.measureSojourn; to: statistics2.parameters[0] }

Duration 100 000: typically about 500 ± 65, and 500 is inside about 9 times in 10. The empty start now weighs little.

Key insight: noise and bias

Replications shrink the noise; only a longer run (or a proper warm-up) removes the bias of the empty start.

• This is how a simulation is verified: rebuild a case with a known answer and check that the interval covers it

Your turn: M/M/2 at the same utilization

Two slow servers or one fast one? (5 min)

Erlang C example, Duration 100 000, $R = 20$. Set $1/\mu = 200$, $k = 2$: still $\rho = 0.8$.

	$k$	$1/\mu$	$P(\text{wait})$ formula	$W$ formula	$W$ simulated ± $\delta$
one fast	1	100
two slow	2	200

Does each formula value lie inside its interval? Which set-up is better for the customer, and why?

The price of precision

Predict: how many runs for half the width?

Question

At $R = 20$ the interval is $\pm\delta$. How many replications give $\pm\delta/2$?

$$\delta = t\cdot\frac{S}{\sqrt{R}} \quad\Longrightarrow\quad \frac{\delta}{2} \text{ needs about } 4R = 80 \text{ replications}$$

$S$ is a property of the model, not of $R$: it stays roughly the same. Only $\sqrt{R}$ works for you.

Your turn: CI width against $R$

Measure the $1/\sqrt{R}$ law (10 min)

Erlang C example, Duration 100 000, Animate off. Record W Sojourn time:

$R$	5	10	20	50	100
mean $\hat\Theta$
radius $\delta$
$\delta\cdot\sqrt{R}$

Sketch $\delta$ against $R$. What happens to the last row?

Run it: what the table shows

Typical values from our runs (Duration 100 000, default inputs):

$R$	5	10	20	50	100
radius $\delta$	≈ 160	≈ 100	≈ 68	≈ 43	≈ 31
$\delta\cdot\sqrt{R}$	≈ 360	≈ 320	≈ 300	≈ 310	≈ 310

• $\delta$ falls fast at first, then flattens; $\delta\sqrt{R}$ stays roughly constant
• A bit higher at $R = 5$: there $t_{4,\,0.975} = 2.78$ instead of about 2

Key insight: precision costs $\sqrt{R}$

The interval shrinks like $1/\sqrt{R}$: half the width costs four times the runs, one more digit costs a hundred times.

• A limit of all sampling methods, not of Quisim

Build it: less variability, fewer runs

Model codeOpen in Quisim

NumberInput { id: number0; initValue:125; label:"1/λ" }
NumberInput { id: number1; initValue:100; label:"D" }
NumberInput { id: number2; initValue:1; label:"seed A" }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
Wire { from: number2.measureValue; to: rnd0.parameterSeed }
Wire { from: number1.measureValue; to: perform0.parameterDuration }
RandomExponential { id:rnd0 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Arrive { id:arrive0 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0 }
Pipe { from:queue0.outlet; to:perform0.inlet }
Perform { id:perform0 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0 }
Wire { from: queue0.measureLoad; to: chart0.parameters[0] }
Statistics { id: statistics0 }
Wire { from: leave0.measureSojourn; to: statistics0.parameters[0] }
Chart { id: chart0 }

Back to your own model. Select the service RandomExponential and its seed knob and delete them (Del); wire the 100 NumberInput straight into Perform D. Now it is M/D/1: theory $W = 300$ ms.

Your turn: how many runs does M/D/1 need?

Same five seeds, less spread (5 min)

Run your M/D/1 with seeds 1, 3, 5, 7, 9 and fill in the table, as you did for M/M/1:

seed	1	3	5	7	9	$S$
mean sojourn

Compare $S$ with your M/M/1 table. How many M/M/1 runs give the same $\delta$ as five M/D/1 runs? (Hint: $R \propto S^2$.)

Expect $S$ about half of M/M/1's: M/M/1 needs about four times the runs.

Your turn: is the difference real?

One fast server or two slow ones? (10 min)

Erlang C example, Duration 100 000. Formula: $W = 500$ ($k = 1$, $1/\mu = 100$) against $W \approx 556$ ($k = 2$, $1/\mu = 200$). Run both with $R = 5$, $20$, $100$:

$R$	$k=1$: $\hat\Theta \pm \delta$	$k=2$: $\hat\Theta \pm \delta$	intervals overlap?
5
20
100

From which $R$ on can the simulation alone tell the two apart?

Key insight: precise enough for the decision

Choose $R$ so that the intervals of the options you compare stop overlapping; a small difference needs many runs.

• Plan a study: a few pilot runs, estimate $S$, then pick $R$ from $\delta = t\,S/\sqrt{R}$
• Other levers: longer runs, a proper warm-up, less variability in the system itself

Wrap-up

What we learned this week

• A run is one sample: a fixed seed replays it, a new seed draws another history
• Replications turn many samples into $\hat\Theta \pm \delta$, $\delta = t\cdot S/\sqrt{R}$
• "95%" describes the procedure: about 1 interval in 20 misses
• Replications fix noise, not bias: short runs from an empty start sit below theory
• With long enough runs, simulation agrees with Erlang C
• Precision is expensive: half the width costs four times the runs

Take-home

M/M/2 at 80% utilization

Erlang C example: $1/\lambda = 125$, $1/\mu = 200$, $k = 2$, Duration 100 000.

1. Run with $R = 5$ and note the radius $\delta$ of W Sojourn time
2. Predict how many replications give an interval half as wide; then find it by trying
3. Is the formula value ($W \approx 556$) inside your intervals?

$R$	5	(your prediction)	(tries)
$\delta$

Bonus: why is the answer to 2. less than $4 \times 5 = 20$? (Hint: $t_{R-1,\,0.975}$.)

Next week

Previous: Week 1: From blocking to queues to multi-server systems

Next: Week 3: Beyond simple queues — priorities and triage

More: Reading the confidence interval · Steering block · Statistics block · Course outline

Preview

Erlang C assumes one FIFO line. In an emergency department, critical patients go first. What does that do to everyone else's wait?

For the lecturer

Running this session

• Timing: recap and M/M/k reruns 15 min, build and hand replications 25 min, Statistics example 10 min, Erlang C and bias 25 min, price of precision 15 min
• Students work in their own browser in the free demo; the slides' Model code boxes are only for full-version users
• Ask for a written prediction before every run; collect the M/M/k ranges and the class experiment's misses on the board
• If time is short, drop the M/D/1 build and the "is the difference real?" task: the take-home covers the same idea

Notes on the models

• All example models are unseeded: every Start clean differs. Only the hand-built model has seed knobs (seed 0 = clock)
• Demo kit used: Arrive, Queue, Perform, Leave, 2 RandomExponential, NumberInputs, Statistics, Chart. The demo kit has no Steering block: replications need the Statistics or Erlang C example
• Erlang C defaults (Duration 10 000, warm-up 5 000) give a biased mean near 420 (Quisim test runs, 300 replications: 422): about 1 in 3 intervals miss at $R = 20$, about 2 in 3 at $R = 100$. That is the point of the "bias" slides
• At Duration 100 000 the bias nearly vanishes; the mean is about 500 and coverage is about 90–94% in our runs, a little under 95% because the averages of single runs are skewed
• Replications is capped at 100 in both examples; turn Animate off for $R = 100$