How many doctors, and when?

Discrete-event simulation, told through an emergency department

Quisim

The staffing problem

An emergency department, open 24 hours

• Patients arrive around the clock, but not evenly
• Each doctor treats one patient at a time
• Too few doctors → long, unsafe waits. Too many → paid idle time

The question: how many doctors per shift, so that waiting stays acceptable at the lowest cost?

Arrivals follow the clock

What we know from the records

Shift	Time	Patients per hour	One patient every … ($1/\lambda$)
Night	00–06	2	30 min
Morning	06–12	6	10 min
Afternoon	12–18	10	6 min
Evening	18–24	4	15 min

A treatment takes 30 minutes on average ($1/\mu = 30$), at any time of day.

The target

For each shift: the fewest doctors such that at least 90% of patients wait no more than 20 minutes.

• Not a statement about the average wait, but about the long tail
• We will get there step by step: first the mean wait, then the 90% target

How to follow along in the free demo

• Open the free demo, show the editor pane (pen icon)
• Examples tab: open a ready-made model, e.g. Erlang C or M/M/k
• Blocks tab: build a small model from the kit; drag from socket to socket for pipes (green) and wires (silver)
• Type new values into the NumberInput fields; the other properties stay fixed in the demo
• Animate (film icon) on to watch, off to get results fast; then Start clean (Ctrl+Enter)
• Monitor tab: the event list, one line per arrival, start and end of treatment

Here 1 model time unit = 1 minute. Full-version users can open Model code under each picture and paste it instead.

One shift as a model

Build it: an afternoon in the department

1. Blocks tab: Arrive → Queue → Perform → Leave, joined by pipes (outlet O → inlet I)
2. RandomExponential V → Arrive D (gap between patients); a second one V → Perform D (treatment time)
3. NumberInput "1/λ" = 6 → first RandomExponential µ; "1/µ" = 30 → second µ; "Doctors" = 6 → Perform C

Model codeOpen in Quisim

NumberInput { id:num0; initValue:6; label:"1/λ"; lower:1 }
Wire { from:num0.measureValue; to:rndA.parameterMean }
RandomExponential { id:rndA; seed:42 }
Wire { from:rndA.measureRandomVariate; to:arrive0.parameterDuration }
Arrive { id:arrive0; label:"Entrance" }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; label:"Waiting room" }
Pipe { from:queue0.outlet; to:perform0.inlet }
NumberInput { id:num1; initValue:30; label:"1/µ"; lower:1 }
Wire { from:num1.measureValue; to:rndB.parameterMean }
RandomExponential { id:rndB; seed:77 }
Wire { from:rndB.measureRandomVariate; to:perform0.parameterDuration }
NumberInput { id:num2; initValue:6; label:"Doctors"; lower:1 }
Wire { from:num2.measureValue; to:perform0.parameterCapacity }
Perform { id:perform0; label:"Treatment" }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; label:"Discharge" }

Build it: watch the waiting room

4. Chart: Queue L (patients waiting) → Chart P[0]
5. Statistics: Queue LW (wait of each patient leaving the waiting room) → Statistics P[0]
6. Animate on, Start clean; open the Monitor tab and follow single patients

Model codeOpen in Quisim

NumberInput { id:num0; initValue:6; label:"1/λ"; lower:1 }
Wire { from:num0.measureValue; to:rndA.parameterMean }
RandomExponential { id:rndA; seed:42 }
Wire { from:rndA.measureRandomVariate; to:arrive0.parameterDuration }
Arrive { id:arrive0; label:"Entrance" }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; label:"Waiting room" }
Pipe { from:queue0.outlet; to:perform0.inlet }
NumberInput { id:num1; initValue:30; label:"1/µ"; lower:1 }
Wire { from:num1.measureValue; to:rndB.parameterMean }
RandomExponential { id:rndB; seed:77 }
Wire { from:rndB.measureRandomVariate; to:perform0.parameterDuration }
NumberInput { id:num2; initValue:6; label:"Doctors"; lower:1 }
Wire { from:num2.measureValue; to:perform0.parameterCapacity }
Perform { id:perform0; label:"Treatment" }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; label:"Discharge" }
Wire { from:queue0.measureLoad; to:chart0.parameters[0] }
Chart { id:chart0; label:"Patients waiting"; height:80; width:300 }
Wire { from:queue0.measureWait; to:stats0.parameters[0] }
Statistics { id:stats0; label:"Mean wait"; height:48; width:128 }

A model that plays the day forward

What you just built

• Entities: the patients (the bubbles)
• Queue: the waiting room
• Resource: the doctors (Perform capacity)

What drives it

• Events: arrival, start of treatment, end of treatment
• The clock jumps from event to event: see the Monitor tab

More: How a discrete-event simulation runs

Why the usual tools fall short

The spreadsheet answer

Afternoon: 10 patients/h × 30 min = 5 doctor-hours of work per hour

• Utilisation $\rho = \dfrac{1/\mu}{k \cdot 1/\lambda} = \dfrac{30}{k \cdot 6}$: with $k = 5$ doctors, $\rho = 100\%$
• "Keep up with demand" rule: 5 doctors exactly busy, so 6 to be safe

Predict: five doctors, exactly enough work?

In your model set Doctors = 5. Work arriving = work the doctors can do.

What happens to the waiting room over the run?

Run it: the queue never settles. It drifts upwards and hardly ever empties again; the mean wait comes out in the hundreds of minutes (about 220 in the run below, against about 18 with six doctors) and keeps growing with the run length: at $\rho = 100\%$ there is no slack to recover from a busy hour.

Model codeOpen in Quisim

NumberInput { id:num0; initValue:6; label:"1/λ"; lower:1 }
Wire { from:num0.measureValue; to:rndA.parameterMean }
RandomExponential { id:rndA; seed:42 }
Wire { from:rndA.measureRandomVariate; to:arrive0.parameterDuration }
Arrive { id:arrive0; label:"Entrance" }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; label:"Waiting room" }
Pipe { from:queue0.outlet; to:perform0.inlet }
NumberInput { id:num1; initValue:30; label:"1/µ"; lower:1 }
Wire { from:num1.measureValue; to:rndB.parameterMean }
RandomExponential { id:rndB; seed:77 }
Wire { from:rndB.measureRandomVariate; to:perform0.parameterDuration }
NumberInput { id:num2; initValue:5; label:"Doctors"; lower:1 }
Wire { from:num2.measureValue; to:perform0.parameterCapacity }
Perform { id:perform0; label:"Treatment" }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; label:"Discharge" }
Wire { from:queue0.measureLoad; to:chart0.parameters[0] }
Chart { id:chart0; label:"Patients waiting"; height:80; width:300 }
Wire { from:queue0.measureWait; to:stats0.parameters[0] }
Statistics { id:stats0; label:"Mean wait"; height:48; width:128 }

Predict: six doctors, 83% busy

Back to Doctors = 6. The doctors are idle 17% of the time.

Mean wait in the waiting room: about 2, 5, 20 or 60 minutes?

Run it: the queue now comes back down, but the mean wait is still around 15–20 minutes (three runs gave about 12, 16 and 20; theory: 17.6 min). Each Start clean gives a different value: one run is one possible week, not "the" answer.

Your turn: how many doctors for the afternoon?

In your own model (1/λ = 6, 1/µ = 30), Animate off, Start clean three times per row:

Doctors	ρ	Mean wait, run 1	run 2	run 3	Highest queue in the Chart
5	100%
6	83%
7	71%
8	63%

Question: from which number on do the three runs agree well? Why do they disagree more at 6 than at 8?

Key insight: waiting comes from variability

Waiting is caused by variability, not by the averages: patients come in clusters, some treatments run long, and the busier the team, the longer a surge takes to clear.

• Animate the model with 6 doctors: watch a cluster of arrivals fill the waiting room
• The spreadsheet sees only averages, so it cannot see this

The queueing formula: Erlang C

One shift is an M/M/k queue

Random arrivals (M), random treatment times (M), $k$ doctors: open Examples → Erlang C.

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 Function block computes the textbook formulas; the blocks next to it simulate the same queue.

What the formula says

• Offered load $a = \lambda/\mu$, utilisation $\rho = a/k$; stable only if $\rho < 1$
• $P(\text{wait})$: share of patients who find all $k$ doctors busy (Erlang C)
• Mean wait $W_q = \dfrac{P(\text{wait})}{k\mu - \lambda}$

Erlang C gives exact answers in seconds, but only for one steady queue with exponential times.

Predict: the afternoon with six doctors

In Erlang C set 1/λ = 6, 1/µ = 30, k servers = 6.

Which share of patients must wait at all? 10%, 30% or 60%?

Run it: $P(\text{wait}) = 0.59$, $\rho = 0.83$, $W_q = 17.6$ min. The simulated mean wait (Statistics, 20 replications) comes out at about 19 ± 5 min: the formula's 17.6 lies inside the interval.

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:6; label:"1/λ\ninterarrival\ntime"; lower:1; x:22; y:48; height:72; width:96 }
NumberInput { id: number1; initValue:30; label:"1/µ\nservice time"; x:22; y:164; lower:1; height:58; width:96 }
NumberInput { id: number2; initValue:6; 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] }

Your turn: the afternoon in numbers

Erlang C example, 1/λ = 6, 1/µ = 30, Animate off. Set k servers and Start clean:

k	ρ	P(wait)	Wq formula	Wq simulated
5
6
7
8

1. What is the smallest stable k? What do the formula boxes and the Chart show at k = 5?
2. What is the smallest k with a mean wait below 5 minutes?

Check: the afternoon

k	ρ	P(wait)	$W_q$ (formula)
5	100%	formula boxes show nan; simulated Wq in the hundreds of minutes, the queue chart climbs	∞
6	83%	0.59	17.6 min
7	71%	0.32	4.9 min
8	63%	0.17	1.7 min

One extra doctor (6 → 7) cuts the mean wait by more than two thirds.

Careful: the mean is not our target

The Erlang C example reports the mean wait and $P(\text{wait})$. Our target is the share of patients seen within 20 minutes, which the example does not show.

For M/M/k there is one more formula: $P(W_q \le 20) = 1 - P(\text{wait})\, e^{-(k\mu - \lambda)\cdot 20}$

Afternoon: $k = 6$ → 70%, $k = 7$ → 91.5%, $k = 8$ → 97.7%

Your turn: night, morning, evening

Group work, one shift per group. Erlang C example, 1/µ = 30:

Group	Shift	1/λ	k to try
A	Night	30	1, 2, 3, 4
B	Morning	10	3, 4, 5, 6
C	Evening	15	2, 3, 4, 5

Record ρ, P(wait) and Wq for each k. Then compute $P(W_q \le 20)$ with the formula above. Report: the smallest stable k, and the smallest k that reaches 90% within 20 minutes.

Check: all four shifts

Shift	1/λ	Smallest stable k	P(wait) at that k	Smallest k with ≥ 90% within 20 min
Night	30	2	0.33	3 (97.6%)
Morning	10	4	0.51	5 (93.8%)
Afternoon	6	6	0.59	7 (91.5%)
Evening	15	3	0.44	4 (95.4%)

Key insight: every shift needs one doctor more than "just stable", and stable is already one more than the spreadsheet's workload.

And beyond the formula?

• Formulas exist only for simple, idealised systems: one queue, a load that never changes, a long-run steady state
• A real day changes its load every few hours; triage, emergency priority, several stations: no formula

→ We need to replay the department, with all its randomness. More: Why simulate?

The full model: measuring the target

The emergency department as blocks

Each patient's wait is recorded; the target is met when 90% of the waits are 20 minutes or less.

Measuring the service level

"At least 90% wait ≤ 20 min" is the same as "the 90% quantile of the wait is ≤ 20 min".

Model codeOpen in Quisim

RandomExponential { id:rndA; mean:6; seed:42; x:22; y:48 }
Wire { from:rndA.measureRandomVariate; to:arrive0.parameterDuration }
Arrive { id:arrive0; label:"Entrance"; x:162; y:48 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; label:"Waiting room"; x:295; y:48 }
Pipe { from:queue0.outlet; to:perform0.inlet }
RandomExponential { id:rndB; mean:30; seed:77; x:287; y:210 }
Wire { from:rndB.measureRandomVariate; to:perform0.parameterDuration }
Number { id:cap0; value:6; x:440; y:210 }
Wire { from:cap0.measureValue; to:perform0.parameterCapacity }
Perform { id:perform0; label:"Treatment"; x:433; y:48 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; label:"Discharge"; x:570; y:48 }
Wire { from:queue0.measureWait; to:stats0.parameters[0] }
Statistics { id:stats0; x:640; y:100; height:120; width:200; label:"90% wait at most"; warmupWidth:500; quantile:0.9 }
Steering { id:steer0; duration:10500; replications:50; x:22; y:210 }

Statistics with quantile 0.9 on the waits: with 6 doctors, 90% of patients wait at most ≈ 54 ± 7 min (50 replications), far above 20. The quantile and Steering settings need the full version.

Randomness, on purpose

• Arrival gaps and treatment times are drawn at random around the recorded averages
• One run is one possible course of events, so the model runs 50 times (replications)
• The first 500 minutes start from an empty waiting room and are discarded (warm-up)

Result: not one number, but an estimate with a confidence interval: 54 ± 7 min. The M/M/k formula gives 53 min for the 90% quantile (and only 70% of patients within 20 min).

Predict: how many doctors for 90%?

With 6 doctors, 90% of afternoon patients wait up to about 54 min. How many doctors does the simulation need to bring this to 20 minutes or less?

Run it: set the Number block to 7 → about 18 ± 1.5 min, below 20 (formula: 17.6 min, i.e. 91.5% within 20 min). One doctor more cuts the 90% wait to a third.

Scenarios and the grid search

Asking "what if?" systematically

• Scenario = one shift × one number of doctors: 4 × 10 = 40 scenarios
• Steering runs 40 editions in one go; a Function turns the edition number $e$ into one scenario
• For each scenario: the 90% wait with a confidence interval (batch means)

// e = 0 … 39: shift = e / 10 (rounded down), doctors = e % 10 + 1
M = [[30, 10, 6, 15][Math.floor(P[0] / 10)], P[0] % 10 + 1]

The grid search model

Model codeOpen in Quisim

Steering { id:steer0; duration:10500; editions:40 }
Wire { from:steer0.measureEdition; to:grid0.parameters[0] }
Function { id:grid0; label:"Scenario"; numMeasures:2; code:"M = [[30, 10, 6, 15][Math.floor(P[0] / 10)], P[0] % 10 + 1]" }
Wire { from:grid0.measures[0]; to:rndA.parameterMean }
Wire { from:grid0.measures[1]; to:perform0.parameterCapacity }
RandomExponential { id:rndA; seed:42 }
Wire { from:rndA.measureRandomVariate; to:arrive0.parameterDuration }
Arrive { id:arrive0; label:"Entrance" }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; label:"Waiting room" }
Pipe { from:queue0.outlet; to:perform0.inlet }
RandomExponential { id:rndB; mean:30; seed:77 }
Wire { from:rndB.measureRandomVariate; to:perform0.parameterDuration }
Perform { id:perform0; label:"Treatment" }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; label:"Discharge" }
Wire { from:grid0.measures[0]; to:stats0.parameters[0] }
Wire { from:grid0.measures[1]; to:stats0.parameters[1] }
Wire { from:queue0.measureWait; to:stats0.parameters[2] }
Statistics { id:stats0; label:"90% wait at most"; warmupWidth:500; numParameters:3; numBatches:10; quantile:0.9; height:120; width:220 }

Function M[0] → arrival mean, M[1] → doctors; both also go to Statistics to label each edition. 1/λ = 6: 6 doctors ≈ 42 ± 22 min, 7 doctors ≈ 18 ± 11 min, 8 doctors ≈ 7 min.

Doctors needed per shift

Shift	Time	"Keep up" rule	Simulation (≥ 90% in 20 min)
Night	00–06	2	3
Morning	06–12	4	5
Afternoon	12–18	6	7
Evening	18–24	3	4
Total		15	19 doctor-shifts / day

The simulation confirms what your Erlang C groups found.

What the decision rests on

• The busiest shift needs the most doctors, but not in proportion to patients
• Uncertainty is largest exactly at the threshold, where the 90% wait crosses 20 minutes
• Far below or far above the threshold the answer is clear-cut

Key insight: the plan is the starting point for scheduling. Breaks, handovers, absences and part-time contracts come on top.

A whole day: shift changes

The plan 3/5/7/4 in one model

Model codeOpen in Quisim

Steering { id:steer0; duration:4320; replications:30 }
Arrive { id:shift0; label:"Shift change"; beginAtZero:false; duration:360 }
Pipe { from:shift0.outlet; to:shiftEnd0.inlet }
Leave { id:shiftEnd0 }
Wire { from:shift0.measureNumPassed; to:plan0.parameters[0] }
Function { id:plan0; label:"Staffing plan"; numMeasures:6; code:"const s = P[0] % 4;\nM = [[30, 10, 6, 15][s], [3, 5, 7, 4][s], s == 0, s == 1, s == 2, s == 3]" }
Wire { from:plan0.measures[0]; to:rndA.parameterMean }
Wire { from:plan0.measures[1]; to:perform0.parameterCapacity }
RandomExponential { id:rndA; seed:42 }
Wire { from:rndA.measureRandomVariate; to:arrive0.parameterDuration }
Arrive { id:arrive0; label:"Entrance" }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; label:"Waiting room" }
Pipe { from:queue0.outlet; to:perform0.inlet }
RandomExponential { id:rndB; mean:30; seed:77 }
Wire { from:rndB.measureRandomVariate; to:perform0.parameterDuration }
Perform { id:perform0; label:"Treatment" }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; label:"Discharge" }
Wire { from:queue0.measureWait; to:stAll.parameters[0] }
Statistics { id:stAll; quantile:0.9; label:"Whole day"; height:100; width:150 }
Wire { from:queue0.measureWait; to:st0.parameters[0] }
Wire { from:plan0.measures[2]; to:st0.parameterActive }
Statistics { id:st0; quantile:0.9; label:"Night"; height:100; width:150 }
Wire { from:queue0.measureWait; to:st1.parameters[0] }
Wire { from:plan0.measures[3]; to:st1.parameterActive }
Statistics { id:st1; quantile:0.9; label:"Morning"; height:100; width:150 }
Wire { from:queue0.measureWait; to:st2.parameters[0] }
Wire { from:plan0.measures[4]; to:st2.parameterActive }
Statistics { id:st2; quantile:0.9; label:"Afternoon"; height:100; width:150 }
Wire { from:queue0.measureWait; to:st3.parameters[0] }
Wire { from:plan0.measures[5]; to:st3.parameterActive }
Statistics { id:st3; quantile:0.9; label:"Evening"; height:100; width:150 }

Every 360 minutes a "Shift change" event switches arrival mean and number of doctors. M[2…5] switch on the Statistics block of the current shift only (by the time treatment starts). 30 runs of 3 days each.

Predict: does the plan hold over a whole day?

Each shift reaches 90% within 20 minutes on its own. Joined into one day: which shift's patients fare worst?

Run it: over the whole day 90% wait at most about 15 min: the target holds on average. But the evening misses it badly, about 30 ± 12 min (night ≈ 6, morning ≈ 7, afternoon ≈ 12): at 18:00 the team drops from 7 to 4, while the afternoon queue is still there.

Your turn: fix the evening

Full version (in the demo: discuss on paper). In the Staffing plan Function change [3, 5, 7, 4]:

Plan	90% wait, whole day	Night	Evening	Doctor-shifts
3, 5, 7, 4				19
3, 5, 7, 5				20
2, 5, 7, 5				19

Question: is one more evening doctor worth it, or can the night give one away? What would a handover overlap (afternoon team stays 1 hour longer) need in the model?

Key insight: shifts are not independent

A queue carries over from one shift into the next: staffing each shift "on its own" misses the moment the team gets smaller.

More: Time-dependent arrival rates

What the model leaves out

Still missing

• Treatment times assumed exponential; arrivals random within a shift
• No triage, no emergency priority, no breaks

The case study extends the model: longer, more realistic treatment times (Gamma, mean 50 min) and a prioritised emergency stream.

The extended model: emergency priority

Model codeOpen in Quisim

Number { id: number0; value:15; x:22; y:48 }
Wire { from: number0.measureValue; to: rnd0.parameterMean }
RandomExponential { id:rnd0; seed:88; x:162; y:48 }
Arrive { id:arrive0; label:"Entrance"; x:310; y:48 }
Pipe { from:arrive0.outlet; to:queue0.inlet }
Queue { id:queue0; label:"Waiting room"; x:451; y:48 }
Number { id: number1; value:25; x:20; y:191 }
Wire { from: number1.measureValue; to: rnd1.parameterScale }
Number { id: number2; x:91; y:191; value:2; height:48; width:64 }
Wire { from: number2.measureValue; to: rnd1.parameterShape; bend:-6; incline:-18; debend:-18; decline:-19 }
RandomGamma { id:rnd1; seed:143; x:162; y:191 }
Wire { from:rnd0.measureRandomVariate; to:arrive0.parameterDuration }
Wire { from:rnd1.measureRandomVariate; to:perform0.parameterDuration }
Perform { id:perform0; capacity:4; label:"Treatment"; x:630; y:185 }
Pipe { from:perform0.outlet; to:leave0.inlet }
Leave { id:leave0; label:"Exit"; x:718; y:185 }
Note { id: note0; x:123; y:191; label:"Treatment time" }
Note { id: note1; x:123; y:49; label:"Inter-arrival time" }
Wire { from: queue0.measureWait; to: chart0.parameters[0] }
Chart { id: chart0; x:586; y:7; height:61; width:185; label:"Waiting time" }
Wire { from: queue0.measureWait; to: statistics0.parameters[0] }
Statistics { id: statistics0; x:586; y:96; height:35; width:184; warmupWidth:10000 }
Number { id: number3; value:180; x:22; y:345 }
RandomExponential { id:rnd2; seed:1231; x:162; y:345 }
Arrive { id:arrive1; label:"Emergency entrance"; x:310; y:345 }
Queue { id:queue1; label:"Emergency waiting"; x:451; y:345 }
Function { id: function0; x:540; y:354; code:"M[0] = P[0]>0 ? 1 : 0" }
Wire { from: queue1.measureLoad; to: function0.parameters[0] }
Pipe { from: queue0.outlet; to: select0.inlets[0] }
Wire { from: function0.measures[0]; to: select0.parameters[0] }
Select { id: select0; x:540; y:185 }
Pipe { from: select0.outlet; to: perform0.inlet }
Pipe { from: queue1.outlet; to: select0.inlets[1] }
Note { id: note2; x:123; y:342; label:"Emergency\ninter-arrival time" }
Wire { from: number3.measureValue; to: rnd2.parameterMean }
Pipe { from:arrive1.outlet; to:queue1.inlet }
Wire { from:rnd2.measureRandomVariate; to:arrive1.parameterDuration }

Select takes the emergency queue whenever it is not empty. Extended model

Can we trust it?

Four questions to ask of any simulation

• Data: where do arrival and treatment figures come from? How old are they?
• Validation: does the model reproduce today's waiting times with today's staffing?
• Confidence intervals: is the difference between two options bigger than the noise?
• Warm-up: the model starts with an empty waiting room, so the first stretch is discarded

More: Model testing · Output analysis

Your turn: validate before you trust

Erlang C example, afternoon values (1/λ = 6, 1/µ = 30, k = 6). Set Replications to 1, 5, 20 and Start clean:

Replications	Wq simulated	Width of the interval
1		(none: one run has no ±)
5
20

Question: how many replications until the simulated Wq reliably agrees with the formula's 17.6 min?

Where else this applies

Same pattern, different setting

Anything with arrivals, queues and limited capacity:

• Call centres and service desks: agents per hour
• Check-in, security, ticket counters: open lanes per time slot
• Production lines: buffers, machines, shifts (extrusion case study)
• Logistics and warehouses: pickers, docks, stock levels

Next steps

• Read the full case, with every model and exercise: Hospital case study
• Learn the building blocks: Introduction to Quisim (slides) and the Quisim handbook
• Make arrivals follow the clock within one model: Time-dependent arrival rates

Appendix

A1 Data and parameters

• Time unit: minutes
• Arrivals: exponential gaps, mean $\mu_A$ = 30 / 10 / 6 / 15 min (night / morning / afternoon / evening)
• Treatment: exponential, mean $\mu_B = 30$ min
• Workload rule: $c_{\min} = \lceil \mu_B / \mu_A \rceil$ → 2 / 4 / 6 / 3
• Utilisation $\rho = \mu_B / (c \cdot \mu_A)$, e.g. afternoon with $c = 6$: $30 / 36 = 83\%$

Data in the case study

A2 The base model, block by block

• RandomExponential (mean $\mu_A$) → Arrive "Entrance"
• Queue "Waiting room" → Perform "Treatment"
• RandomExponential (mean 30) → Perform duration; Number → Perform capacity $c$
• Perform → Leave "Discharge"
• Queue measureWait → Statistics with quantile: 0.9 (the wait that 90% of patients do not exceed)
• Not via a Function M[0] = P[0]<=20 ? 1 : 0: a Function only re-evaluates when its input changes, so runs of zero waits are dropped and the share comes out too low (≈ 52% instead of 70%)

Step 2 in the case study

A3 Run settings

• Single scenario: Steering duration 10,500, 50 replications
• Statistics warmupWidth: 500 (first 500 min discarded)
• Afternoon, $c = 6$: 90% wait ≈ 54 ± 7 min; $c = 7$: ≈ 18 ± 1.5 min

A4 Grid search over editions

// Function "Scenario", P[0] = Steering EN (edition 0 … 39)
M = [[30, 10, 6, 15][Math.floor(P[0] / 10)],  // inter-arrival mean
     P[0] % 10 + 1]                            // number of doctors

• M[0] → RandomExponential mean; M[1] → Perform capacity
• M[0], M[1] also go to Statistics P[0], P[1] to label each edition
• The case study does the same with an Explorer block (for … yield)

A5 Grid search statistics

• Steering editions: 40, duration 10,500
• Statistics numParameters: 3, numBatches: 10, warmupWidth: 500
• Confidence interval per edition by batch means: 10 batches of 1,000 min each
• quantile: 0.9 on P[2] = Queue measureWait

Step 4 in the case study

A6 Erlang C (M/M/k) formulas

• $a = \lambda/\mu$, $\rho = a/k < 1$
• $P(\text{wait}) = \dfrac{\frac{a^k}{k!}\,\frac{1}{1-\rho}}{\sum_{n=0}^{k-1}\frac{a^n}{n!} + \frac{a^k}{k!}\,\frac{1}{1-\rho}}$
• $W_q = \dfrac{P(\text{wait})}{k\mu - \lambda}$, $\quad P(W_q \le t) = 1 - P(\text{wait})\,e^{-(k\mu-\lambda)\,t}$
• The Erlang C example's Function block evaluates the first two (Jagerman's recursion)