Model real decisions with linear and integer programming, networks, queues, simulation and decision trees, then present the results so managers act on them.

Optimization & Operations Research

## Why it matters

Queues are everywhere: patients in an emergency department, calls waiting for an agent, trucks at a loading dock, jobs waiting for a server. Customers feel waiting long before they notice anything else, and managers always ask the same questions: how long will people wait, and how many servers do we need? **Queueing theory** answers them with a handful of formulas, and it explains a surprising fact: waits explode as a system gets close to full.

## Describing a queue

**Kendall notation** writes a queue as `A/S/c`:

- `A` is the arrival pattern. `M` (for Markovian, or memoryless) means a Poisson process: arrivals are random and independent, with an average rate `lambda` per unit time, and the gaps between them are exponential.
- `S` is the service time distribution. `M` again means exponential with average rate `mu` (average service time `1/mu`).
- `c` is the number of identical servers.

So **M/M/1** is random arrivals, random service and one server, served first come, first served, with unlimited waiting room. It is the simplest useful model and the baseline for everything else.

The key number is the **utilization** (traffic intensity):

```text
rho = lambda / mu     (fraction of time the server is busy)
```

If `rho >= 1`, work arrives faster than it can be done and the queue grows without limit. The system is **unstable**. A stable queue needs `lambda < mu`.

## The M/M/1 formulas

For a stable M/M/1 queue:

| quantity | meaning | formula |
|---|---|---|
| `L` | average number in the system | `rho / (1 - rho)` |
| `Lq` | average number waiting | `rho^2 / (1 - rho)` |
| `W` | average time in the system | `1 / (mu - lambda)` |
| `Wq` | average time waiting | `rho / (mu - lambda)` |
| `P(wait)` | chance an arrival waits | `rho` |

**Little's law** ties them together and holds for almost any stable queue, not just M/M/1:

```text
L = lambda * W        Lq = lambda * Wq        W = Wq + 1/mu
```

It is a powerful sanity check: if 10 customers arrive per hour and each spends half an hour in the system, there are 5 in the system on average.

## Worked example: a pharmacy counter

Customers arrive at 10 per hour. One pharmacist serves 12 per hour (5 minutes each on average).

```text
rho = 10 / 12 = 0.833
L   = 0.833 / 0.167 = 5 customers
Lq  = 0.833^2 / 0.167 = 4.17 customers
W   = 1 / (12 - 10) = 0.5 hours = 30 minutes
Wq  = 0.833 / 2 = 0.417 hours = 25 minutes
```

A 5-minute service turns into a 25-minute wait. Watch what happens as arrivals rise, with service fixed at 12 per hour:

| arrivals per hour | 6 | 9 | 10 | 11 | 11.5 |
|---|---|---|---|---|---|
| utilization | 50% | 75% | 83% | 92% | 96% |
| average wait (min) | 5 | 15 | 25 | 55 | 115 |

Going from 92% to 96% busy doubles the wait. This **non-linear blow-up** is the most important practical lesson of queueing theory: running people or machines at "100% efficiency" means very long queues.

## More servers: M/M/c and Erlang C

With `c` servers, each serving at rate `mu`, the system is stable when `lambda < c * mu`. Let `a = lambda / mu` (the load in "erlangs") and `rho = a / c`. The **Erlang C** formula gives the chance an arrival must wait:

```text
P(wait) = [a^c / c! / (1 - rho)] / [ sum_{k=0}^{c-1} a^k / k!  +  a^c / c! / (1 - rho) ]
Wq      = P(wait) / (c * mu - lambda)
```

Adding a second pharmacist at the same counter: `P(wait) = 0.245` and `Wq` is about **1.05 minutes**, down from 25. Pooling servers in one shared line is far better than separate lines.

Call centres use Erlang C to staff every half hour. With 100 calls per hour and 12 calls per agent per hour, 9 agents is the minimum for stability; 10 agents give an average wait of about 1.5 minutes and 11 agents about 0.6 minutes. If the target is "average wait under 1 minute", the answer is 11.

## When the formulas do not fit

The formulas assume a steady state, Poisson arrivals and exponential service. Real systems have rush hours, scheduled appointments, customers who give up (abandonment), priorities and service times that are far from exponential. Use the formulas for quick estimates and staffing rules of thumb, and switch to **simulation** (next lesson) when the details matter.

## Common mistakes

- **Mixing time units.** `lambda` and `mu` must use the same unit (per hour, or per minute).
- **Applying formulas to an unstable queue.** If `lambda >= mu`, the "answer" is negative or infinite.
- **Planning for 100% utilization.** Waits grow without limit as `rho` approaches 1.
- **Confusing W and Wq.** Time in system includes service; time in queue does not.

## Quick recap

- M/M/1: Poisson arrivals at `lambda`, exponential service at `mu`, one server; stable only if `lambda < mu`.
- `Wq = rho / (mu - lambda)`, and Little's law `L = lambda * W` links the rest.
- Waits explode as utilization nears 100%.
- Erlang C sizes multi-server systems such as call centres.

## Try it

In the **Do** tab, code the M/M/1 metrics and Erlang C, then size a call centre.


Compute utilization, queue length and waiting time for an M/M/1 queue, and size a multi-server system with Erlang C.

Queueing models

Queuing Theory Tutorial - Queues/Lines, Characteristics, Kendall Notation, M/M/1 Queues

Queueing Theory: M/M/1 Queue Worked Example

## Why it matters

Queueing formulas need tidy assumptions. Real service lines have lunch rushes, two kinds of customers, a server who takes breaks and a machine that breaks down. When the maths gets messy, OR analysts **simulate**: they build a small model of the system in code, run it thousands of times with random inputs, and measure what happens. Discrete-event simulation (DES) is used to design hospitals, factories, airports and warehouses, and to test a schedule or layout before anyone spends money on it.

## The core idea

A **discrete-event** simulation jumps from one event to the next instead of ticking through every second. Nothing in the model changes between events.

- **Entities**: the things that flow (customers, patients, parts).
- **Resources**: what entities compete for (servers, beds, machines).
- **Events**: moments when the state changes (an arrival, a service start, a departure).
- **State**: what you track (queue length, which servers are busy).
- **Clock**: the current simulated time.
- **Event list**: future events sorted by time. The simulation repeatedly pops the earliest event, advances the clock to it, updates the state, and schedules any new events it causes.

A one-server line has two event types. An **arrival** joins the queue (or starts service if the server is free) and schedules the next arrival. A **departure** frees the server and starts the next person in line.

## A one-server line in a few lines of Python

For one server, first come first served, you do not even need an explicit event list. Each customer starts service when they arrive or when the server becomes free, whichever is later:

```python
def simulate_line(arrivals, services):
    out = []
    server_free_at = 0.0
    for arrive, service in zip(arrivals, services):
        start = max(arrive, server_free_at)
        finish = start + service
        out.append((start, finish, start - arrive))   # wait = start - arrival
        server_free_at = finish
    return out
```

Check it by hand. Customers arrive at times 0, 1, 2 and 5 and need 2, 2, 1 and 1 minutes:

| customer | arrives | starts | finishes | waits |
|---|---|---|---|---|
| 1 | 0 | 0 | 2 | 0 |
| 2 | 1 | 2 | 4 | 1 |
| 3 | 2 | 4 | 5 | 2 |
| 4 | 5 | 5 | 6 | 0 |

Hand-checking a tiny case like this is **verification**: making sure the code does what you meant.

## Adding randomness

Random inputs come from probability distributions. For an M/M/1 coffee kiosk with 1 arrival per minute and 2 services per minute:

```python
import random

rng = random.Random(42)            # fixed seed: the same run every time
t = 0.0
arrivals, services = [], []
for _ in range(50000):
    t += rng.expovariate(1.0)      # time to the next arrival
    arrivals.append(t)
    services.append(rng.expovariate(2.0))
```

Theory says `Wq = rho / (mu - lambda) = 0.5 / 1 = 0.5` minutes and half of all customers wait. A long simulated run lands close to both numbers. Comparing a simulation with a case you can solve exactly is how you build trust before you add the messy details that the formulas cannot handle: a rush hour, a second barista from 8 to 10, or customers who leave when the line is longer than six.

**Always fix the seed** while developing, so results are reproducible and bugs can be found. Then change the seed to get independent runs.

## Running experiments properly

- **Warm-up.** A simulation that starts empty is not typical. Discard the first stretch of output (the warm-up period) or start from a realistic state.
- **Replications.** One run is one random story. Run many independent replications (different seeds) and report the average with a confidence interval, for example `mean +/- 1.96 * sd / sqrt(n)`.
- **Same random numbers for comparisons.** When comparing two designs, give both the same seeds (common random numbers), so differences come from the design and not from luck.
- **Validation.** Compare the model's output with real data (average waits last month). Verification asks "is the code right?"; validation asks "is the model right?"

## Tools

In Python, **SimPy** is a popular open-source DES library. Processes are generator functions, and resources handle the queueing for you:

```python
import simpy

def customer(env, counter, service_time, waits):
    arrive = env.now
    with counter.request() as req:
        yield req                          # wait for the server
        waits.append(env.now - arrive)
        yield env.timeout(service_time)    # being served

env = simpy.Environment()
counter = simpy.Resource(env, capacity=1)
```

Commercial tools such as Arena, AnyLogic, FlexSim and Simio add drag-and-drop modelling and animation, which helps a lot when you present to managers.

## Common mistakes

- **Trusting one run.** Report several replications and a confidence interval.
- **No warm-up.** Starting empty makes waits look shorter than they really are.
- **Unverified code.** Always test a tiny hand-checked case and a case with a known formula.
- **Wrong rates.** `expovariate(rate)` takes a rate (per minute), not a mean time.
- **Modelling every detail on day one.** Start simple, validate, then add detail.

## Quick recap

- DES jumps between events: arrivals, service starts, departures.
- For one FIFO server, start = max(arrival, server free), wait = start - arrival.
- Fix seeds while developing; use warm-up, replications and confidence intervals for results.
- Verify against hand cases and formulas; validate against real data.

## Try it

In the **Do** tab, simulate a service line and check it against the M/M/1 formula.


Simulate a single-server service line with random arrivals and check the results against queueing theory.

Discrete-event simulation

Understanding Discrete Event Simulation, Part 1: What Is Discrete Event Simulation

Simulating a Queue: Basic Discrete Event Simulation

## Why it matters

Some decisions are made once, under real uncertainty: build a big plant or a small one, launch now or run a pilot first, bid on a contract or not. There is no long-run average to optimize, and nobody knows what demand will be. **Decision analysis** gives a disciplined way to lay out the choices, the uncertain events and the payoffs, and to compare options fairly. A decision tree is also one of the best tools you have for explaining a recommendation to managers, because they can see every assumption.

## Building blocks of a decision tree

- **Decision nodes** (drawn as squares): you choose one branch.
- **Chance nodes** (drawn as circles): nature chooses a branch, with probabilities that sum to 1.
- **End points**: the payoff of that path, usually profit or cost.

The tree is drawn left to right in the order things happen, and solved **right to left** (rolling back):

- At a chance node, take the **expected value**: the sum of probability times value.
- At a decision node, take the **best** branch and note which one it is.

The expected value of an option is often called its **EMV** (expected monetary value).

## Worked example: plant size

A company must decide how to serve a new market. Demand will be strong with probability 0.6 or weak with probability 0.4. Profits are in $k.

- **Build large**: 200 if demand is strong, -60 if weak.
- **Build small**: if demand is strong, the company can then **expand** (profit 140) or **stay** small (profit 90). If demand is weak, profit is 20.
- **Wait**: profit 0.

```text
[decide] -- large -- (demand) -- 0.6 strong --> 200
         |                   \-- 0.4 weak   --> -60
         |-- small -- (demand) -- 0.6 strong --> [decide] -- expand --> 140
         |                   |                            \-- stay   -->  90
         |                   \-- 0.4 weak   --> 20
         \-- wait --> 0
```

Roll back from the right.

```text
small, strong demand: best of expand (140) and stay (90) = 140
EMV(small) = 0.6 * 140 + 0.4 * 20 = 84 + 8 = 92
EMV(large) = 0.6 * 200 + 0.4 * (-60) = 120 - 24 = 96
EMV(wait)  = 0
```

The EMV-best choice is **build large** (96), with small close behind (92). The plan is a full strategy: "build small, and expand if demand is strong" is how the small option is valued.

## Sensitivity: how sure must we be?

Let `p` be the probability of strong demand.

```text
EMV(large) = 200p - 60(1 - p) = 260p - 60
EMV(small) = 140p + 20(1 - p) = 120p + 20
```

They are equal when `140p = 80`, so `p = 0.571`. Our estimate of 0.6 is only just above the break-even point. If the true probability is 0.5, small wins (80 versus 70). That is worth telling the decision maker: "the choice hinges on whether you believe strong demand is more than about 57% likely."

## The value of information

Suppose a market study could tell us the demand in advance. With **perfect information** we would pick the best option for each outcome:

```text
strong: max(200, 140, 0) = 200
weak:   max(-60, 20, 0)  = 20
EV with perfect information = 0.6 * 200 + 0.4 * 20 = 128
EVPI = 128 - 96 = 32
```

The **expected value of perfect information** (EVPI) is 32 ($k): never pay more than that for any study, because no real study is perfect. A real, imperfect study is valued with Bayes' rule (the expected value of sample information), using what you learned about conditional probabilities.

## Beyond expected value

EMV assumes the decision maker is **risk neutral**. The large plant loses 60 with probability 0.4; a small firm might not survive that. Two ways to handle risk:

- Show the **risk profile**: the distribution of outcomes for each option, not just the average.
- Use a **utility function** that values losses more heavily than gains, and roll back expected utility instead of expected money.

Probabilities come from data, models (like the simulations in the last lesson) or structured expert judgement. Write down where each one came from.

## Common mistakes

- **Rolling back in the wrong direction.** Always solve from the end points back to the first decision.
- **Probabilities that do not sum to 1** at a chance node.
- **Double counting costs**, for example subtracting the plant cost at the end point and again on the branch.
- **Ignoring later decisions.** The option to expand is what makes the small plant attractive.
- **Reporting only the EMV.** Show the break-even probabilities and the worst case too.

## Quick recap

- Decision trees use decision nodes (choose the best) and chance nodes (take expected values).
- Roll back from right to left to find the best strategy and its EMV.
- Sensitivity analysis finds the break-even probabilities.
- EVPI caps what information is worth; risk attitude may justify a lower-EMV option.

## Try it

In the **Do** tab, write a recursive function that rolls back a decision tree and picks the best option.


Build a decision tree, roll it back with expected values, and measure how much information is worth.

Decision analysis

Decision Analysis 3: Decision Trees

Decision Tree Analysis - Intro and Example with Expected Monetary Value

Estimate waits with queueing formulas, simulate a service line, and compare risky options with decision trees.