10. Applied OR Models

The earlier chapters built general tools: linear and integer programming, network algorithms, dynamic programming and stochastic models. This chapter applies them to four problems that every factory, warehouse, construction site and service operation faces: how much stock to hold, when each activity of a project must happen, in what order jobs should go through machines, and how much demand to expect. These models are among the oldest and most used in OR. Their closed-form answers and simple rules run inside ERP systems, project-management software and shop-floor schedulers, so knowing where they come from, and when their assumptions fail, is directly useful in practice.

<a id="101-inventory-theory"></a>

10.1 Inventory Theory

Inventory is a buffer between supply and demand. Holding it costs money, and running out of it costs money too. This subchapter starts with the deterministic economic order quantity and two extensions (quantity discounts and finite production rate). It then adds uncertainty with reorder points and safety stock, the single-period newsvendor model, and the $(s, S)$ policies that are optimal for periodic review with fixed ordering costs. The first three topics share one item: annual demand $D = 1200$ units, ordering cost $K = 100$ per order and holding cost $h = 6$ per unit per year, the EOQ example of Section 5.1.

Economic Order Quantity Model

Why it matters. Every replenishment decision balances two opposing costs. Ordering often in small lots keeps little stock on the shelf but pays the fixed cost of an order (paperwork, transport, receiving, a machine setup) many times. Ordering rarely in large lots saves those fixed costs but ties up capital and space in stock. The economic order quantity (EOQ), published by Harris in 1913, finds the best balance. It is the starting point of all inventory theory, and its square-root law still sets lot sizes in many ERP systems.

Theory. The EOQ model makes strong assumptions:

• demand is known and constant at $D$ units per year,
• each order costs a fixed $K$, regardless of its size,
• holding one unit for one year costs $h$ (often a rate $i$ times the unit value $c$: capital, storage, insurance, obsolescence),
• orders arrive all at once, after a known lead time, and shortages are not allowed.

Under these assumptions the optimal policy orders the same quantity $Q$ every time the stock reaches the level that covers the lead time. Inventory follows a sawtooth:

inventory
  Q +\        |\        |\
    | \       | \       | \
    |  \      |  \      |  \
Q/2 +...\.....|...\.....|...\.....   average inventory = Q/2
    |    \    |    \    |    \
  R +-----\---|-----\---|-----\---   reorder point: an order is placed here
    |      \  |      \  |      \
  0 +-------\-|-------\-|-------\--> time
     <-T=Q/D->
          <-L->   lead time: the order placed at R arrives just as stock hits 0

Each cycle lasts $T = Q/D$ years, so there are $D/Q$ orders per year, and the average stock is $Q/2$. The annual cost, excluding the purchase cost $cD$ that does not depend on $Q$, is

$$ C(Q) = \underbrace{\frac{K D}{Q}}{\text{ordering}} + \underbrace{\frac{h Q}{2}}{\text{holding}}. $$

$C$ is convex for $Q > 0$. Setting $C'(Q) = -KD/Q^2 + h/2 = 0$ gives the EOQ

$$ Q^* = \sqrt{\frac{2 K D}{h}}, \qquad C(Q^*) = \sqrt{2 K D h}. $$

At the optimum the two costs are equal, $KD/Q^* = hQ^*/2$. This balance is a useful sanity check on any lot size.

Three properties make the EOQ more useful than its assumptions suggest:

• Square-root law. Doubling demand increases the lot size by only $\sqrt{2} \approx 1.41$. Consolidating the demand of several warehouses into one therefore reduces total cycle stock.
• Flat cost curve. Dividing $C(Q)$ by $C(Q^)$ gives $$ \frac{C(Q)}{C(Q^)} = \frac{1}{2}\left(\frac{Q}{Q^} + \frac{Q^}{Q}\right). $$ A lot 50% too large costs only about 8% more. Errors in $K$, $D$ or $h$ are damped twice, once by the square root and once by the flat bottom. Lot sizes can be rounded to pallets or truckloads at little cost.
• Reorder point. With a constant lead time $L$, order when the stock covers the lead-time demand: $R = D L$. If $L$ is longer than a cycle, compare $R$ with the inventory position (stock on hand plus stock on order), not with the stock on hand.

An analogy is grocery shopping. Each trip costs time (the fixed cost $K$), and a full fridge ties up money and spoils food (holding $h$). People who live far from the store shop rarely and buy a lot, which is the square-root law in action.

When several items are ordered together, or ordering intervals must fit a calendar, power-of-two policies restrict each cycle to a base period times $2^k$. Roundy (1985) showed that they cost at most 6% more than the unrestricted optimum, thanks to the flat cost curve.

Example: lot size and its sensitivity. The code computes the EOQ for the running item and shows how slowly the cost rises around it.

import numpy as np

K, D, h = 100.0, 1200.0, 6.0      # cost per order, demand per year, holding cost per unit per year

def cost(Q):
    return K * D / Q + h * Q / 2  # annual ordering cost + annual holding cost

Q_star = np.sqrt(2 * K * D / h)   # square-root formula
print(f"Q* = {Q_star:.0f}, orders/yr = {D / Q_star:.1f}, "
      f"cycle = {Q_star / D * 52:.1f} weeks, C(Q*) = {cost(Q_star):.0f}")

for Q in (100, 150, 200, 250, 300, 400):
    print(f"Q = {Q:3d}: ordering {K * D / Q:5.0f} + holding {h * Q / 2:5.0f}"
          f" = {cost(Q):5.0f}  ({cost(Q) / cost(Q_star) - 1:+.1%})")

# Deterministic lead time of 2 weeks: reorder when the stock covers the lead time
L = 2 / 52
print(f"reorder point R = D * L = {D * L:.0f} units")
# Q* = 200, orders/yr = 6.0, cycle = 8.7 weeks, C(Q*) = 1200
# Q = 100: ordering  1200 + holding   300 =  1500  (+25.0%)
# Q = 150: ordering   800 + holding   450 =  1250  (+4.2%)
# Q = 200: ordering   600 + holding   600 =  1200  (+0.0%)   <- the two costs are equal
# Q = 250: ordering   480 + holding   750 =  1230  (+2.5%)
# Q = 300: ordering   400 + holding   900 =  1300  (+8.3%)
# Q = 400: ordering   300 + holding  1200 =  1500  (+25.0%)
# reorder point R = D * L = 46 units

The optimal policy orders 200 units six times a year, whenever the stock falls to 46 units. Lots between 150 and 250 all cost within 4.2% of the optimum, and halving or doubling the lot costs 25% more. The curve is asymmetric: lots that are too small hurt more than lots that are too large by the same amount.

Key takeaways

• The EOQ $Q^* = \sqrt{2KD/h}$ balances fixed ordering costs against holding costs. At the optimum the two are equal.
• The cost curve is flat near $Q^*$, so the EOQ is robust to errors in its inputs and to rounding.
• Lot size grows with the square root of demand, which is why consolidating demand reduces cycle stock.
• The reorder point covers lead-time demand. Compare it with the inventory position when orders overlap.

Practice

1. (Beginner) A shop sells 5,000 units per year. Each order costs 50 and holding costs 2 per unit per year. Compute the EOQ, the number of orders per year, the cycle length in weeks and the annual cost. Check that ordering and holding costs are equal at the optimum.
2. (Intermediate) The ordering cost $K$ in the running example was estimated as 100, but the true value is 200. The firm keeps ordering 200 units. Using the ratio formula, by what percentage does its cost exceed the cost of the true optimum? Explain why the error is so small.
3. (Interview) "A supplier offers to cut our ordering cost in half by automating purchase orders. How much will our total inventory-related cost fall?" Hint: $C^* = \sqrt{2KDh}$ depends on $\sqrt{K}$. Also think about what the lower $K$ does to the lot size and to warehouse space.

EOQ with Quantity Discounts

Why it matters. Suppliers rarely charge one price for every order size. They give discounts for larger orders to share their own economies of scale. A discount can justify ordering far more than the EOQ, but it also raises holding costs. Buyers need a quick, correct way to decide whether a price break is worth it.

Theory. Once the price depends on $Q$, the purchase cost $c(Q),D$ must be added to the objective because it is no longer constant:

$$ TC(Q) = c(Q), D + \frac{K D}{Q} + \frac{i, c(Q), Q}{2}, $$

where the holding cost $h = i,c$ is a rate $i$ times the price paid. Two discount schemes are common:

Scheme	How it works	Cost curve
All-units	the price for the whole order drops once $Q$ reaches a break point $q_k$	a separate curve per tier, with downward jumps at the breaks
Incremental	only the units beyond $q_k$ get the lower price, like tax brackets	continuous, with kinks at the breaks

For all-units discounts, the optimum is found tier by tier:

1. For each tier $k$ with price $c_k$ and range $[q_k, q_{k+1})$, compute $\text{EOQ}_k = \sqrt{2KD/(i c_k)}$.
2. Clamp it to the range: if $\text{EOQ}_k < q_k$, use $q_k$ (the smallest quantity that earns the discount); if it is above the range, use the upper end, though that tier is then beaten by the next cheaper tier.
3. Evaluate $TC$ at each candidate and keep the cheapest.

The break points are often optimal, because the cost curves are convex within each tier and the only reason to leave a tier's own EOQ is to reach a cheaper price. For incremental discounts the purchase cost per unit becomes an average that falls with $Q$. The candidates are then the tier EOQs computed with this average cost, and break points are never optimal.

The analogy is a bulk pack at the supermarket. The larger pack is cheaper per unit, but it is only a good deal if you can store it and use it before it spoils.

Example: three price tiers. The running item costs 24 per unit for orders below 300, 23.50 from 300 to 999, and 23.00 from 1,000. The holding rate is 25% per year, so $h = 6$ at the base price.

import numpy as np

K, D, i = 100.0, 1200.0, 0.25           # order cost, annual demand, holding rate (25% of price per year)
tiers = [(0, 24.0), (300, 23.5), (1000, 23.0)]   # (minimum order quantity, unit price), all-units

best = None
for k, (q_min, c) in enumerate(tiers):
    q_max = tiers[k + 1][0] - 1 if k + 1 < len(tiers) else np.inf
    h = i * c                            # holding cost depends on the price paid
    eoq = np.sqrt(2 * K * D / h)
    Q = min(max(eoq, q_min), q_max)      # push the EOQ into the tier's range
    total = c * D + K * D / Q + h * Q / 2   # purchase cost now matters: it differs by tier
    print(f"price {c:5.2f}: EOQ = {eoq:5.0f} -> Q = {Q:5.0f}, "
          f"purchase {c * D:6.0f} + order {K * D / Q:4.0f} + hold {h * Q / 2:5.0f} = {total:6.0f}")
    if best is None or total < best[1]:
        best = (Q, total)
print(f"best: order {best[0]:.0f} units, annual cost {best[1]:.0f}")
# price 24.00: EOQ =   200 -> Q =   200, purchase  28800 + order  600 + hold   600 =  30000
# price 23.50: EOQ =   202 -> Q =   300, purchase  28200 + order  400 + hold   881 =  29481
# price 23.00: EOQ =   204 -> Q =  1000, purchase  27600 + order  120 + hold  2875 =  30595
# best: order 300 units, annual cost 29481

Moving from 200 to 300 units saves 600 per year on purchases and 200 on orders, and costs only 281 in extra holding. The largest discount is not worth it: reaching 1,000 units saves another 600 on purchases but adds about 2,000 in holding cost. Notice that the EOQ barely changes across tiers (200 to 204). The discount matters through the break points, not through the formula.

Key takeaways

• With price breaks, the purchase cost $c(Q) D$ enters the objective and the holding cost $i,c(Q)$ changes by tier.
• For all-units discounts, evaluate each tier's EOQ clamped to its range. The optimum is often exactly at a break point.
• Bigger discounts are not always better: holding costs grow linearly with the lot, while the savings from the price are fixed per unit.
• Incremental discounts lead to a smooth problem whose optimum is never at a break point.

Practice

1. (Beginner) Repeat the example with a holding rate of 10% per year instead of 25%. Which tier wins now, and why does a lower holding rate favor larger orders?
2. (Intermediate) Modify the code for incremental discounts: the first 299 units of every order cost 24, units 300 to 999 cost 23.50 and units from 1,000 on cost 23.00. Write the average unit cost as a function of $Q$ and find the optimal order size numerically.
3. (Interview) "A supplier offers 2% off if we order a full quarter's demand at once. Our holding rate is 20% per year. Should we accept?" Hint: compare the annual saving $0.02,cD$ with the extra holding cost of moving from the EOQ to $D/4$. Also consider obsolescence risk, cash flow and storage space, which the formula leaves out.

Production Lot Size Model

Why it matters. When a firm makes the item itself, the "order" is a production run. The units do not arrive all at once: they flow off the line at a finite rate while customers keep taking units away. Stock therefore builds up more slowly and peaks lower than in the EOQ model. This changes the best batch size. The model also leads to a central question in factories that make many products on one line: how often should the line switch between them?

Theory. The economic production quantity (EPQ) model keeps the EOQ assumptions except that a run of $Q$ units is produced at rate $P > D$ units per year. A setup costs $K$. During the run, which lasts $t_p = Q/P$, inventory grows at rate $P - D$. Afterwards it falls at rate $D$:

inventory
                 peak = Q (1 - D/P)
         /\                  /\
        /  \                /  \           up:   P - D   (producing and selling)
       /    \              /    \          down: D       (selling only)
      /      \            /      \
  ---+--------\----------+--------\------> time
     |<-t_p->|
     |<------- T = Q/D ->|

The peak is $Q(1 - D/P)$ and the average stock is half of it, so the annual cost is

$$ C(Q) = \frac{K D}{Q} + \frac{h Q}{2}\left(1 - \frac{D}{P}\right), \qquad Q^* = \sqrt{\frac{2 K D}{h,(1 - D/P)}}, \qquad C^* = \sqrt{2 K D h \left(1 - \frac{D}{P}\right)}. $$

The EPQ is the EOQ with an effective holding cost $h(1 - D/P)$. Two limits check the formula:

• $P \to \infty$: production is instantaneous, and the EOQ is recovered.
• $P \to D$: the line runs all the time just to keep up, inventory never builds up, and the lot becomes infinite. The line simply never switches off.

The ratio $\rho = D/P$ is the share of line time this product needs. When several products share one line, their runs must not overlap. That is the economic lot scheduling problem (ELSP), which is NP-hard in general. The classic heuristic is the common cycle: the line makes every product once per cycle of length $T$. Lot sizes are $Q_i = D_i T$, and the cost

$$ C(T) = \frac{\sum_i K_i}{T} + \frac{T}{2}\sum_i h_i D_i (1 - \rho_i) \quad\Longrightarrow\quad T^* = \sqrt{\frac{2 \sum_i K_i}{\sum_i h_i D_i (1 - \rho_i)}} $$

is again an EOQ-type trade-off. It is feasible only if $\sum_i \rho_i < 1$ and the cycle leaves room for the setup times, $\sum_i s_i \le T,(1 - \sum_i \rho_i)$.

Example: effect of the production rate, and a three-product rotation. The running item is now made in-house. The code varies the production rate, then plans a common cycle for three products sharing the line.

import numpy as np

K, D, h = 100.0, 1200.0, 6.0          # setup cost, demand per year, holding cost per unit per year

def epq(P):
    rho = D / P                       # fraction of time the line is busy with this product
    Q = np.sqrt(2 * K * D / (h * (1 - rho)))
    return Q, Q * (1 - rho), np.sqrt(2 * K * D * h * (1 - rho))   # lot, peak stock, annual cost

for P in (np.inf, 9600, 4800, 2400, 1500):
    Q, peak, C = epq(P)
    print(f"P = {P:>5.0f}/yr: lot = {Q:4.0f}, run = {Q / P * 52:4.1f} wk, "
          f"cycle = {Q / D * 52:4.1f} wk, peak stock = {peak:3.0f}, cost = {C:5.0f}")

# Three products on one line, rotated in a common cycle T (the line makes each once per cycle)
Ks = np.array([100.0, 80.0, 150.0]); Ds = np.array([1200.0, 2000.0, 600.0])
hs = np.array([6.0, 2.0, 10.0]);     Ps = np.array([4800.0, 8000.0, 3000.0])
rho = Ds / Ps
print(f"utilization = {rho.sum():.2f}")                    # must be < 1
T = np.sqrt(2 * Ks.sum() / np.sum(hs * Ds * (1 - rho)))    # optimal common cycle (years)
print(f"common cycle T* = {T * 52:.1f} weeks, lots = {np.round(Ds * T)}, "
      f"cost = {Ks.sum() / T + np.sum(hs * Ds * (1 - rho)) * T / 2:.0f}")
indep = sum(np.sqrt(2 * k * d * hh * (1 - r)) for k, d, hh, r in zip(Ks, Ds, hs, rho))
print(f"independent EPQs (ignoring clashes) = {indep:.0f}")
# P =   inf/yr: lot =  200, run =  0.0 wk, cycle =  8.7 wk, peak stock = 200, cost =  1200   <- EOQ
# P =  9600/yr: lot =  214, run =  1.2 wk, cycle =  9.3 wk, peak stock = 187, cost =  1122
# P =  4800/yr: lot =  231, run =  2.5 wk, cycle = 10.0 wk, peak stock = 173, cost =  1039
# P =  2400/yr: lot =  283, run =  6.1 wk, cycle = 12.3 wk, peak stock = 141, cost =   849
# P =  1500/yr: lot =  447, run = 15.5 wk, cycle = 19.4 wk, peak stock =  89, cost =   537
# utilization = 0.70
# common cycle T* = 11.6 weeks, lots = [268. 447. 134.], cost = 2952
# independent EPQs (ignoring clashes) = 2932

As the line slows down, lots get larger but the peak stock gets smaller, because more of each run is sold while it is being produced. At $P = 1500$ the line spends 80% of its time on this product. For the three-product line, the common cycle costs only 20 more per year (0.7%) than the three independent EPQs. Those independent EPQs have cycles of different lengths, so their runs would collide on the line. The common cycle gives up almost nothing and guarantees a feasible rotation, provided the 30% of free line time covers the setups.

Key takeaways

• With a finite production rate $P$, the EPQ is the EOQ with an effective holding cost $h(1 - D/P)$, so lots are larger and peak stock is lower.
• $D/P$ is the share of line time the product needs. As it approaches 1, the line should simply run continuously.
• Several products on one line create the economic lot scheduling problem. A common rotation cycle is a simple, near-optimal, always-consistent heuristic.
• Setup times consume capacity: a feasible cycle must leave time for them after production.

Practice

1. (Beginner) A bottling line can fill 60,000 bottles per month. Demand is 20,000 per month, a setup costs 900 and holding one bottle for a month costs 0.02. Compute the EPQ, the run length, the cycle length and the peak inventory.
2. (Intermediate) In the three-product example, every setup takes the same time $s$. What is the longest $s$ (in working days, 5 per week) for which the common cycle $T^*$ is still feasible? If setups take twice that long, find the smallest feasible $T$ and the resulting cost increase.
3. (Interview) "Lean manufacturing says to cut setup times so you can run small batches. How does that relate to the EPQ formula?" Hint: shorter setups usually mean a lower setup cost $K$. Follow what $Q^* \propto \sqrt{K}$ does to stock and lead times, and what freed capacity does to the feasibility of a rotation.

Reorder Point and Safety Stock

Why it matters. In reality, demand during the lead time is uncertain. If the reorder point covers only the average lead-time demand, the stock runs out in about half of the cycles. Safety stock is extra stock held to cover demand above the average during the lead time. Setting it is one of the most frequent quantitative decisions in supply chains, and errors in it show up directly as stock-outs or as excess capital.

Theory. In a continuous review $(Q, R)$ policy, the inventory position is watched continuously, and an order of $Q$ units is placed when it drops to the reorder point $R$. The order quantity $Q$ usually comes from the EOQ. A stock-out can only happen during the lead time, so only the demand during lead time (DDLT) matters for $R$:

$$ R = \mu_L + \underbrace{z,\sigma_L}_{\text{safety stock}}. $$

With i.i.d. daily demand of mean $d$ and standard deviation $\sigma_d$, and a fixed lead time of $L$ days, $\mu_L = d L$ and $\sigma_L = \sigma_d \sqrt{L}$. The square root comes from adding independent variances, so safety stock grows more slowly than the lead time. If the lead time is itself random, with mean $\bar L$ and standard deviation $s_L$, the DDLT is a random sum (Section 9.1) and

$$ \sigma_L = \sqrt{\bar L, \sigma_d^2 + d^2 s_L^2}. $$

The second term is often the larger one. An unreliable supplier can require more safety stock than volatile demand.

The safety factor $z$ is set from a service target. The two common targets measure different things:

Measure	Definition	Formula (normal DDLT)
Cycle service level (CSL, $\alpha$)	probability that a replenishment cycle has no stock-out	$\alpha = \Phi(z)$, so $z = \Phi^{-1}(\alpha)$
Fill rate ($\beta$)	fraction of demand served directly from stock	$\beta = 1 - \sigma_L, G(z) / Q$

Here $G(z) = \varphi(z) - z,(1 - \Phi(z))$ is the standard normal loss function, the expected shortfall $E[(Z - z)^+]$. The fill rate depends on $Q$: with large orders there are few cycles per year, and therefore few opportunities to run out. A 90% CSL can still give a 99% fill rate. Customers feel the fill rate, so quoting the CSL to them can be misleading.

The model assumes the DDLT is roughly normal, which is reasonable when the lead time covers many days of demand. For slow-moving items (a few units per lead time), a Poisson or negative binomial distribution should be used instead. In practice, $\sigma_d$ should be the standard deviation of the forecast errors (Section 10.4) rather than of demand itself, because forecasts already explain part of the variation.

Example: service targets and an unreliable lead time. Daily demand has mean 20 and standard deviation 6, the lead time is 9 days and orders are 400 units. The code computes reorder points for several service levels, then makes the lead time random (uniform on 6 to 12 days) and checks the result by simulation.

import numpy as np
from scipy.stats import norm

d, sd = 20.0, 6.0                    # daily demand: mean and standard deviation
L, Q = 9, 400                        # lead time (days), order quantity

def G(z):                            # standard normal loss function E[(Z - z)^+]
    return norm.pdf(z) - z * (1 - norm.cdf(z))

mu_L, sigma_L = d * L, sd * np.sqrt(L)   # demand during the lead time
print(f"lead-time demand: mean {mu_L:.0f}, sd {sigma_L:.0f}")
for csl in (0.80, 0.90, 0.95, 0.99):
    z = norm.ppf(csl)
    fill = 1 - sigma_L * G(z) / Q    # fraction of demand served from stock
    print(f"CSL {csl:.0%}: z = {z:.2f}, safety stock = {z * sigma_L:5.1f}, "
          f"R = {mu_L + z * sigma_L:5.1f}, fill rate = {fill:.4f}")

# Random lead time: uniform on {6, ..., 12} days (mean 9, sd 2)
sL = 2.0
sigma_rand = np.sqrt(L * sd**2 + d**2 * sL**2)   # random-sum variance (Section 9.1)
R = mu_L + norm.ppf(0.95) * sigma_rand
print(f"random lead time: sd {sigma_rand:.1f}, safety stock = {R - mu_L:.1f}, R = {R:.1f}")

# Monte Carlo check of the cycle service level P(lead-time demand <= R)
rng = np.random.default_rng(seed=7)
n = 200_000
lt = rng.integers(6, 13, size=n)                      # one lead time per replenishment cycle
ddlt = np.array([rng.normal(d, sd, k).sum() for k in lt])
print(f"simulated CSL with R = {R:.1f}: {np.mean(ddlt <= R):.4f}")
# lead-time demand: mean 180, sd 18
# CSL 80%: z = 0.84, safety stock =  15.1, R = 195.1, fill rate = 0.9950
# CSL 90%: z = 1.28, safety stock =  23.1, R = 203.1, fill rate = 0.9979
# CSL 95%: z = 1.64, safety stock =  29.6, R = 209.6, fill rate = 0.9991
# CSL 99%: z = 2.33, safety stock =  41.9, R = 221.9, fill rate = 0.9998
# random lead time: sd 43.9, safety stock = 72.1, R = 252.1
# simulated CSL with R = 252.1: 0.9521

Going from 95% to 99% CSL needs 41% more safety stock, a typical sign of diminishing returns at high service levels. The fill rates are all above 99%, because each cycle sells 400 units and a stock-out, when it happens, is short. A lead time that varies by only ±3 days multiplies the safety stock by 2.4 (from 29.6 to 72.1 units), even though average demand and the average lead time are unchanged. The simulation confirms the 95% target (95.2%), so the normal approximation works well here.

Key takeaways

• The reorder point is mean lead-time demand plus safety stock $z,\sigma_L$. Only demand during the lead time can cause a stock-out.
• $\sigma_L$ grows with $\sqrt{L}$ for a fixed lead time. Lead-time variability adds $d^2 s_L^2$ to the variance and often dominates.
• Cycle service level and fill rate are different targets. The fill rate also depends on the order quantity.
• Use the standard deviation of forecast errors, not of raw demand, and a discrete distribution for slow movers.

Practice

1. (Beginner) Weekly demand is normal with mean 150 and standard deviation 40, and the lead time is 4 weeks. Compute the reorder point for a 97.5% cycle service level.
2. (Intermediate) In the example, find the reorder point that gives a fill rate of exactly 99.5% with $Q = 400$, by solving $\sigma_L G(z) = 0.005,Q$ for $z$ numerically. Which CSL does it correspond to? Repeat with $Q = 100$.
3. (Interview) "We moved production from a local supplier (lead time 1 week, very reliable) to an overseas one (6 weeks, plus or minus 2). The unit cost fell 8%. What happens to our inventory?" Hint: compute safety stock under both, add pipeline stock ($d \bar L$, owned in transit), and compare the extra holding cost with the purchase saving.

Newsvendor Model

Why it matters. Some products can be sold only during one short period: newspapers, fresh food, fashion collections, seasonal toys, airline seats, conference catering. There is a single ordering decision before demand is known. Ordering too little loses profitable sales, and ordering too much leaves goods that must be scrapped or marked down. Section 1.1 showed that ordering the mean demand is wrong. The newsvendor model gives the exact answer in closed form, and the same logic sizes capacity, overbooking and safety stock.

Theory. Let $D$ be random demand with cdf $F$, $p$ the selling price, $c$ the unit cost and $v < c$ the salvage value of a leftover unit. Define the two costs of being wrong by one unit:

• underage cost $c_u = p - c$: the margin lost on a sale that could not be made,
• overage cost $c_o = c - v$: the loss on a unit bought but not sold.

Marginal analysis gives the optimum. Consider whether to add one more unit to an order of $Q$. It will be sold if $D > Q$, which happens with probability $1 - F(Q)$, and it will be left over otherwise. The extra unit is worth adding while its expected gain beats its expected loss:

$$ c_u,(1 - F(Q)) \ge c_o,F(Q) \quad\Longleftrightarrow\quad F(Q) \le \frac{c_u}{c_u + c_o}. $$

The optimal order therefore satisfies

$$ F(Q^*) = \frac{c_u}{c_u + c_o}, \qquad \text{the } \textbf{critical ratio}. $$

It is the critical-ratio quantile of demand. For discrete demand, take the smallest $Q$ with $F(Q) \ge c_u/(c_u + c_o)$. With an empirical distribution, that is a quantile of past demand.

The intuition is in the ratio. If a lost sale costs four times as much as a leftover, it pays to be right 80% of the time: order the 80th percentile. Cheap leftovers and high margins push the order above the mean, while expensive leftovers and thin margins push it below.

For normal demand $N(\mu, \sigma^2)$, $Q^* = \mu + z\sigma$ with $z = \Phi^{-1}(c_u/(c_u + c_o))$. The loss function $G$ of the previous topic gives the expected lost sales $\sigma G(z)$, the expected leftovers $Q - \mu + \sigma G(z)$ and the expected profit

$$ E[\text{profit}] = p,(\mu - \sigma G(z)) + v,(Q - \mu + \sigma G(z)) - c,Q. $$

The newsvendor is also the simplest two-stage stochastic program (Section 9.4), and the critical-ratio formula appears in many disguises:

Application	Underage	Overage
Base-stock level per period	backorder cost $b$	holding cost $h$
Airline overbooking (how many extra seats to sell)	an empty seat's revenue	compensation to a bumped passenger
Capacity of a plant built before demand is known	lost margin	idle capacity cost

Conversely, a service-level target $P(D \le Q) \ge 1 - \alpha$, as in a chance constraint, is a newsvendor whose critical ratio has been set directly to $1 - \alpha$ instead of being derived from costs.

Example: the kiosk of Section 1.1, then normal demand with returns.

import numpy as np
from scipy.stats import norm

# (1) The kiosk of Section 1.1: price 10, cost 2, no salvage, demand uniform on [50, 150]
p, c, v = 10.0, 2.0, 0.0
cu, co = p - c, c - v                     # underage: lost margin; overage: loss on a leftover
cr = cu / (cu + co)
print(f"uniform: critical ratio = {cr:.2f}, Q* = {50 + cr * 100:.0f}")

# (2) Normal demand (mean 100, sd 30) and leftovers returned for 1 each
mu, sigma, v = 100.0, 30.0, 1.0
cu, co = p - c, c - v
cr = cu / (cu + co)
z = norm.ppf(cr)
Q = mu + z * sigma
G = norm.pdf(z) - z * (1 - norm.cdf(z))   # standard normal loss function
lost = sigma * G                          # expected unmet demand E[(D - Q)^+]
left = Q - mu + lost                      # expected leftovers E[(Q - D)^+]
profit = p * (mu - lost) + v * left - c * Q
print(f"normal: CR = {cr:.3f}, z = {z:.2f}, Q* = {Q:.1f}, "
      f"lost = {lost:.2f}, leftovers = {left:.2f}, E[profit] = {profit:.1f}")

rng = np.random.default_rng(seed=3)
D = np.maximum(rng.normal(mu, sigma, 1_000_000), 0)
def sim(q):
    return np.mean(p * np.minimum(q, D) + v * np.maximum(q - D, 0) - c * q)
for q in (100, 120, round(Q), 150):
    print(f"  order {q:3d}: simulated E[profit] = {sim(q):.1f}")

# (3) Discrete demand from 20 past days: Q* = smallest q with F(q) >= CR
history = np.array([83, 112, 95, 130, 101, 77, 144, 99, 108, 121,
                    90, 117, 86, 135, 104, 93, 126, 111, 70, 98])
q_emp = np.sort(history)[int(np.ceil(cr * len(history))) - 1]
print(f"empirical: CR = {cr:.3f} -> order {q_emp} (the {int(np.ceil(cr * 20))}th smallest of 20)")
# uniform: critical ratio = 0.80, Q* = 130
# normal: CR = 0.889, z = 1.22, Q* = 136.6, lost = 1.61, leftovers = 38.23, E[profit] = 748.9
#   order 100: simulated E[profit] = 692.4
#   order 120: simulated E[profit] = 739.3
#   order 137: simulated E[profit] = 749.0      <- matches the closed form
#   order 150: simulated E[profit] = 744.8
# empirical: CR = 0.889 -> order 130 (the 18th smallest of 20)

The closed form reproduces the search result of Section 1.1 (order 130) without simulating anything. With a return value of 1, leftovers are cheaper, so the critical ratio rises to 0.889 and the kiosk orders 137 papers, 37 more than the mean demand. On average it is left with 38 unsold papers, and that is optimal. The profit curve is flat near the optimum: ordering 120 or 150 loses less than 2%. Ordering the mean, 100, loses 7.6%. The empirical version needs no distributional assumption, but with only 20 days of history its upper quantiles are imprecise.

Key takeaways

• Order the critical-ratio quantile of demand: $F(Q^*) = c_u / (c_u + c_o)$, with $c_u$ the lost margin and $c_o$ the loss per leftover.
• The optimum exceeds the mean when underage costs more than overage, and vice versa. Ordering the mean is optimal only when the two costs are equal (for symmetric demand).
• For normal demand, the loss function gives expected lost sales, leftovers and profit in closed form.
• The same quantile logic sets base-stock levels, overbooking limits, capacity and chance constraints.

Practice

1. (Beginner) A bakery sells croissants for 2.50. They cost 0.80 to make, and leftovers go to a food bank (value 0). Daily demand is normal with mean 200 and standard deviation 40. How many should it bake? How does the answer change if leftovers can be sold the next day at 1.00?
2. (Intermediate) Prove the critical-ratio formula formally: write $E[\text{profit}(Q)]$ for a continuous demand density, differentiate it using $\frac{d}{dQ} E[\min(Q, D)] = 1 - F(Q)$, and show that the second derivative is negative.
3. (Interview) "A fashion retailer orders once per season with a 60% gross margin, and unsold items are sold at 70% off. Its demand forecast has a large error. Should the order be above or below the forecast?" Hint: compute $c_u$ and $c_o$ from the margin and the markdown, then the critical ratio. Discuss how the forecast's uncertainty, not only its mean, changes the order.

(s, S) Policies

Why it matters. Many items are reviewed periodically (every day, every week) instead of continuously, and every order has a fixed cost. With a fixed cost, ordering every period up to a fixed level wastes money on small orders. Ordering only when stock is low, but then ordering a lot, should be better. The $(s, S)$ policy formalizes this: it is the provably optimal policy for a broad class of inventory problems, and its "min/max" form is built into most ERP systems.

Theory. At the start of each period, the inventory position $x$ is reviewed:

$$ \text{if } x \le s: \text{ order } S - x \text{ (up to } S\text{)}, \qquad \text{otherwise: do not order.} $$

$s$ is the reorder level (min) and $S$ the order-up-to level (max). Two special cases are familiar:

• $K = 0$ (no fixed cost): ordering every period is free, so the optimal policy is a base-stock policy, $s = S - 1$. Its level $S$ is the newsvendor quantile of demand over one period plus the lead time, with critical ratio $b/(b + h)$ for backorder cost $b$ and holding cost $h$.
• Continuous review: the analog is the $(Q, R)$ policy of the previous topic, with $s = R$ and orders of size about $Q$.

Optimality. Consider a multi-period problem with a fixed cost $K$ plus a linear cost per unit ordered, a convex expected holding and shortage cost per period, and backlogged demand. Scarf (1960) proved that an $(s, S)$ policy is optimal in every period. The proof is a dynamic program (Section 6.2). The value function is not convex because of the fixed cost, but it is $K$-convex:

$$ K + g(y + a) \ge g(y) + a,\frac{g(y) - g(y - b)}{b} \qquad \text{for all } a \ge 0,\ b > 0. $$

$K$-convexity is the property that makes "do nothing above $s$, jump to $S$ below it" optimal: below $s$, the saving from moving to $S$ is worth the fixed cost; above $s$ it is not.

Computing $s$ and $S$. Exact values come from dynamic programming or from a Markov chain evaluation of each candidate pair. Practical approximations are common:

• $s$ from the base-stock or reorder-point logic (cover demand over the review period plus the lead time at the target service level),
• $S = s + Q$ with $Q$ the EOQ,
• or the "power approximation" of Ehrhardt, a regression fitted to exact solutions.

The analogy is refueling a car. You do not stop for gas every time the tank drops a little (fixed cost: the detour). You wait until the gauge is low ($s$), then fill up completely ($S$).

Example: exact evaluation of $(s, S)$ policies. Weekly demand is Poisson with mean 10, an order costs $K = 64$, holding costs 1 per unit per week and a backorder costs 9 per unit per week. Orders arrive immediately (zero lead time). With a given $(s, S)$, the post-order stock level $y \in {s+1, \dots, S}$ is a Markov chain (Section 9.2), so the long-run cost of each policy can be computed exactly from the stationary distribution.

import numpy as np
from scipy.stats import poisson

lam, K, h, b = 10.0, 64.0, 1.0, 9.0       # Poisson demand per period, order cost, holding, backorder
d = np.arange(0, 61)                      # demand support (P(D > 60) is negligible)
pmf = poisson.pmf(d, lam)

def L(y):                                 # expected holding + backorder cost after ordering up to y
    return np.sum(pmf * (h * np.maximum(y - d, 0) + b * np.maximum(d - y, 0)))

def avg_cost(s, S):
    """Exact long-run cost per period: Markov chain on the post-order level y in {s+1..S}."""
    levels = np.arange(s + 1, S + 1)
    n = len(levels)
    P = np.zeros((n, n))
    for i, y in enumerate(levels):
        for j, y2 in enumerate(levels):
            if y2 <= y:
                P[i, j] += poisson.pmf(y - y2, lam)      # no order: y2 = y - D
        P[i, -1] += poisson.sf(y - s - 1, lam)           # D >= y - s: drop to <= s, order up to S
    A = np.vstack([P.T - np.eye(n), np.ones(n)])         # solve pi P = pi, sum(pi) = 1
    pi = np.linalg.lstsq(A, np.r_[np.zeros(n), 1.0], rcond=None)[0]
    order_prob = sum(pi[i] * poisson.sf(y - s - 1, lam) for i, y in enumerate(levels))
    return sum(pi[i] * L(y) for i, y in enumerate(levels)) + K * order_prob, order_prob

best = min(((s, S) for s in range(0, 20) for S in range(s + 1, s + 60)),
           key=lambda p: avg_cost(*p)[0])
c, f = avg_cost(*best)
print(f"optimal (s, S) = {best}: cost/period = {c:.2f}, orders in {f:.0%} of periods")

S_bs = int(poisson.ppf(b / (b + h), lam))  # base stock: newsvendor fractile, order every period
c, f = avg_cost(S_bs - 1, S_bs)
print(f"base stock S = {S_bs}: cost/period = {c:.2f}, orders in {f:.0%} of periods")

Q = round(np.sqrt(2 * K * lam / h))        # EOQ heuristic: s from the base stock, S = s + EOQ
c, f = avg_cost(S_bs - 1, S_bs - 1 + Q)
print(f"heuristic (s, S) = ({S_bs - 1}, {S_bs - 1 + Q}): cost/period = {c:.2f}, orders in {f:.0%} of periods")
# optimal (s, S) = (6, 40): cost/period = 35.02, orders in 26% of periods
# base stock S = 14: cost/period = 69.87, orders in 100% of periods
# heuristic (s, S) = (13, 49): cost/period = 39.54, orders in 24% of periods

The optimal policy orders about once every four weeks, when stock falls to 6 or less, and then fills up to 40. The base-stock policy, which would be optimal without the fixed cost, orders every week and costs twice as much, almost all of it in order costs. The simple heuristic ($s$ from the newsvendor fractile, $S = s + \text{EOQ}$) is within 13% of the optimum. Its $s$ is too high. With large orders, most periods start well above $s$, and only the last period before an order is exposed to shortages. Lowering $s$ (and $S$ with it) saves holding cost in every period at the price of a few backorders in that last period.

Key takeaways

• An $(s, S)$ policy orders up to $S$ whenever the inventory position is at or below $s$.
• With a fixed order cost, it is optimal for periodic review with backlogging (Scarf, via $K$-convexity). Without a fixed cost it reduces to a base-stock policy.
• Exact parameters come from dynamic programming or Markov chain evaluation. The EOQ-plus-reorder-point heuristic is a good first approximation.
• $S - s$ plays the role of the order quantity and $s$ the role of the reorder point.

Practice

1. (Beginner) With $s = 20$ and $S = 60$, trace five weeks of a periodic-review system starting at 45 units, for demands 12, 18, 9, 25, 14. Write down when orders are placed and their sizes.
2. (Intermediate) Rerun the example with $K = 16$ and with $K = 256$. How do $s$, $S$ and $S - s$ move? Compare $S - s$ with the EOQ in each case.
3. (Interview) "Our ERP uses min/max planning with min = 2 weeks of demand and max = 6 weeks for every item. What would you change?" Hint: these parameters ignore order costs, demand variability, lead times and service targets. Think about which items gain most from item-specific $(s, S)$ values, and how you would compute them at scale.

<a id="102-project-scheduling"></a>

10.2 Project Scheduling

A project is a one-off set of activities linked by precedence relations: a building, a product launch, a software release, a plant shutdown for maintenance. This subchapter represents projects as networks, computes the shortest schedule and its critical activities with CPM, adds uncertain durations with PERT, trades time against money by crashing, and finally respects limited resources. All topics use the same eight-activity software release:

Activity	Description	Duration (days)	Predecessors
A	requirements	3	-
B	design	5	A
C	backend	8	B
D	frontend	6	B
E	database setup	4	A
F	integration	3	C, D, E
G	user documentation	2	D
H	testing and release	4	F, G

Project Networks

Why it matters. Before computing anything, a project must be written down as activities, durations and "what must finish before what". A drawing of these relations shows which work can run in parallel and which must wait. It also exposes mistakes such as circular dependencies, which make a schedule impossible. Every scheduling algorithm of this subchapter works on this network, so getting its representation right comes first.

Theory. A project consists of activities with durations and precedence relations: "$i \to j$" means $j$ cannot start until $i$ has finished. Precedence relations must form a directed acyclic graph (DAG, Section 3.1). A cycle would mean that an activity waits for itself. Two drawing conventions exist:

Convention	Nodes	Arcs	Notes
Activity-on-node (AON)	activities	precedence relations	no dummies needed; used by modern software (precedence diagramming method)
Activity-on-arrow (AOA)	events (moments in time)	activities	the original PERT/CPM notation; needs dummy activities of zero duration

The AON network of the running example, with a start and an end node added:

                              +--> C (8) -------+
                              |                 v
start --> A (3) --+--> B (5) -+--> D (6) --+--> F (3) --> H (4) --> end
                  |                        |     ^         ^
                  |                        +-----|-> G (2) +
                  |                              |
                  +--> E (4) --------------------+

AOA needs dummies whenever two activities share some, but not all, of their predecessors. Suppose R requires P and Q, while S requires only Q. Activity P ends at event 2 and Q at event 3, and a dummy arc carries Q's completion into event 2:

        P              R
  (1) -----> (2) -----------> (4)
   |          ^
   |  Q       | dummy (duration 0)
   +-------> (3) -----------> (5)
                     S

In AON this is simply four nodes and three arcs. This is why AON dominates in practice.

Real projects also use generalized precedence relations with time lags. The standard one is finish-to-start (FS), "$j$ starts at least $\ell$ days after $i$ finishes". The others are start-to-start (SS), finish-to-finish (FF) and start-to-finish (SF). For example, "testing may start 2 days after coding starts" is SS with lag 2. All of them are linear inequalities between start times, of the form $S_j \ge S_i + \text{lag}$, so the methods below extend to them directly.

Every algorithm processes activities in a topological order: an order in which each activity comes after all its predecessors. Kahn's algorithm builds one in $O(n + m)$ time by repeatedly removing an activity with no unprocessed predecessors. If it gets stuck before placing every activity, the remaining ones contain a cycle. The algorithm thus doubles as a data check.

Example: topological order and cycle detection. The project data below are reused by every later topic of this subchapter.

from collections import deque

# Project data: activity -> (duration in days, predecessors)
project = {
    "A": (3, []),               # requirements
    "B": (5, ["A"]),            # design
    "C": (8, ["B"]),            # backend
    "D": (6, ["B"]),            # frontend
    "E": (4, ["A"]),            # database setup
    "F": (3, ["C", "D", "E"]),  # integration
    "G": (2, ["D"]),            # user documentation
    "H": (4, ["F", "G"]),       # testing and release
}

def topological_order(proj):
    """Kahn's algorithm: repeatedly take an activity whose predecessors are all placed."""
    indeg = {a: len(pred) for a, (_, pred) in proj.items()}
    succ = {a: [] for a in proj}
    for a, (_, pred) in proj.items():
        for p in pred:
            succ[p].append(a)
    ready = deque(a for a in proj if indeg[a] == 0)
    order = []
    while ready:
        a = ready.popleft()
        order.append(a)
        for s in succ[a]:
            indeg[s] -= 1
            if indeg[s] == 0:               # all predecessors of s are now placed
                ready.append(s)
    if len(order) < len(proj):
        raise ValueError(f"precedence cycle among {sorted(set(proj) - set(order))}")
    return order, succ

order, succ = topological_order(project)
print("topological order:", " ".join(order))
print("successors:", {a: s for a, s in succ.items() if s})

bad = dict(project, A=(3, ["H"]))          # data-entry error: A waits for H
try:
    topological_order(bad)
except ValueError as e:
    print("error:", e)
# topological order: A B E C D F G H
# successors: {'A': ['B', 'E'], 'B': ['C', 'D'], 'C': ['F'], 'D': ['F', 'G'], 'E': ['F'], 'F': ['H'], 'G': ['H']}
# error: precedence cycle among ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H']

One wrong predecessor entry ("A waits for H") closes a loop through the whole project, so no activity can ever start. Kahn's algorithm finds no activity without predecessors and reports all of them. Project-management tools run exactly this kind of check whenever a dependency is added.

Key takeaways

• A project is a DAG of activities and precedence relations. A cycle makes the project impossible to schedule.
• Activity-on-node is the standard representation. Activity-on-arrow needs dummy activities to express partial overlaps of predecessors.
• Generalized precedence relations (FS, SS, FF, SF with lags) are all linear inequalities between start times.
• A topological order, computed in linear time by Kahn's algorithm, drives every scheduling pass and detects cycles.

Practice

1. (Beginner) Draw the AOA network of the running example. Which dummy activities are needed, and why?
2. (Intermediate) Extend topological_order to support a lag on each relation, stored as (predecessor, lag) pairs, and compute the earliest start of every activity when integration (F) must start at least 1 day after the database setup (E) ends.
3. (Interview) "How many different topological orders does a project network have, and why might the order matter to a scheduler even though CPM's result does not depend on it?" Hint: count them for a small example, and think about the resource-constrained heuristics later in this subchapter, which schedule activities one by one in list order.

Critical Path Method

Why it matters. Two questions come up in every project review: when can we finish, and which activities must not slip? The critical path method (CPM), developed at DuPont in 1957, answers both with two passes over the network. It also tells how much each non-critical activity can be delayed without delaying the project, which is the room a manager has to move people around.

Theory. With durations $d_j$, CPM computes four times for every activity:

• Forward pass (in topological order): the earliest start $ES_j = \max_{i \to j} EF_i$ (0 if $j$ has no predecessors) and the earliest finish $EF_j = ES_j + d_j$. The project duration is $T = \max_j EF_j$.
• Backward pass (in reverse order): the latest finish $LF_j = \min_{j \to k} LS_k$ ($T$ if $j$ has no successors) and the latest start $LS_j = LF_j - d_j$. These are the latest times that still allow finishing by $T$.

Two kinds of slack follow:

• Total float $TF_j = LS_j - ES_j$: how long $j$ can be delayed without delaying the project. It is shared along a path: if one activity uses it, its successors lose it.
• Free float $FF_j = \min_{j \to k} ES_k - EF_j$: how long $j$ can be delayed without delaying any successor's earliest start. It is never larger than the total float.

Activities with zero total float are critical. They form at least one critical path, a longest path from start to end. Any delay on a critical activity delays the whole project by the same amount. CPM is just the longest-path algorithm for DAGs (Section 3.2). Longest paths are NP-hard in general graphs but linear-time in DAGs, because time only moves forward.

CPM is also a linear program:

$$ \min\ t_{\text{end}} \quad \text{s.t.} \quad t_j \ge t_i + d_i \ \ \forall (i \to j), \qquad t \ge 0. $$

Its dual sends one unit of flow from start to end along the arcs and maximizes the total duration collected, so it picks a longest path. By complementary slackness, the dual variables are 1 on the critical arcs. They are shadow prices: shortening a critical activity by one day shortens the project by one day, until another path becomes critical.

Example: forward and backward passes.

# Reuses project and topological_order() from the Project Networks topic.
order, succ = topological_order(project)
dur = {a: d for a, (d, _) in project.items()}
pred = {a: p for a, (_, p) in project.items()}

ES, EF = {}, {}
for a in order:                                   # forward pass: earliest times
    ES[a] = max((EF[p] for p in pred[a]), default=0)
    EF[a] = ES[a] + dur[a]
T = max(EF.values())                              # project duration

LS, LF = {}, {}
for a in reversed(order):                         # backward pass: latest times
    LF[a] = min((LS[s] for s in succ[a]), default=T)
    LS[a] = LF[a] - dur[a]

print(f"{'act':>3} {'dur':>3} {'ES':>3} {'EF':>3} {'LS':>3} {'LF':>3} {'TF':>3} {'FF':>3}")
for a in order:
    TF = LS[a] - ES[a]                            # total float: delay that keeps T
    FF = min((ES[s] for s in succ[a]), default=T) - EF[a]   # free float: delay that keeps every ES
    print(f"{a:>3} {dur[a]:3d} {ES[a]:3d} {EF[a]:3d} {LS[a]:3d} {LF[a]:3d} {TF:3d} {FF:3d}"
          + ("   critical" if TF == 0 else ""))
print("project duration:", T)
# act dur  ES  EF  LS  LF  TF  FF
#   A   3   0   3   0   3   0   0   critical
#   B   5   3   8   3   8   0   0   critical
#   E   4   3   7  12  16   9   9
#   C   8   8  16   8  16   0   0   critical
#   D   6   8  14  10  16   2   0
#   F   3  16  19  16  19   0   0   critical
#   G   2  14  16  17  19   3   3
#   H   4  19  23  19  23   0   0   critical
# project duration: 23

The release takes 23 days, and the critical path is A-B-C-F-H, driven by the backend. The same schedule as a Gantt chart, with total float shown as dots:

day        0    5    10   15   20   25
           |    |    |    |    |    |
A critical ###
B critical    #####
C critical         ########
D                  ######..
E             ####.........
F critical                 ###
G                        ##...
H critical                    ####

The frontend D has 2 days of total float but 0 free float. Delaying it by a day costs the project nothing, but it pushes the documentation G later, so G's float shrinks from 3 to 2 days. The database setup E can slip by 9 days without affecting anything. Adding developers to E would be wasted. To finish earlier, the backend C, or another critical activity, must get shorter.

Key takeaways

• The forward pass gives earliest times and the project duration. The backward pass gives latest times.
• Total float is the delay an activity can absorb without delaying the project, and free float the delay that keeps every successor's earliest start.
• Critical activities have zero total float and form a longest path. Only shortening them shortens the project.
• CPM is longest path in a DAG, and the dual of its LP identifies the critical arcs.

Practice

1. (Beginner) Recompute the CPM table by hand if the frontend D takes 9 days instead of 6. What is the new critical path and duration?
2. (Intermediate) Solve the CPM linear program of the example with scipy.optimize.linprog and print the dual values (res.ineqlin.marginals). Check that the nonzero ones are exactly the arcs of the critical path.
3. (Interview) "A project has 2,000 activities and three critical paths. The sponsor asks you to cut two weeks. Where do you start?" Hint: activities common to all critical paths shorten all of them at once. Consider float on near-critical paths, which become critical once you start cutting, and the option of changing the logic (overlapping activities) rather than the durations.

PERT and Activity Time Estimation

Why it matters. CPM treats durations as known, but nobody knows exactly how long a new backend will take to build. The Program Evaluation and Review Technique (PERT), created in 1958 for the Polaris submarine missile program, adds uncertainty to activity times. It turns "we will finish in 23 days" into "we have a 60% chance of finishing within 25 days". That is the kind of statement a deadline negotiation needs.

Theory. For each activity, experts give three estimates: optimistic $a$, most likely $m$ and pessimistic $b$. PERT models the duration as a beta distribution on $[a, b]$ with mode $m$, and uses the approximations

$$ \mu = \frac{a + 4m + b}{6}, \qquad \sigma = \frac{b - a}{6}. $$

The weight 4 on the mode and the "range equals six standard deviations" rule are conventions, not exact properties of every beta distribution. Asking for three numbers forces estimators to think about risk, which is the real value of the method.

PERT then makes two assumptions:

1. The critical path computed with mean durations stays critical.
2. Activity times are independent, so the path duration $T$ is a sum of independent variables. By the central limit theorem it is approximately normal with mean $\sum \mu_j$ and variance $\sum \sigma_j^2$ over the critical path: $$ P(T \le t) \approx \Phi!\left(\frac{t - \sum_{j \in CP} \mu_j}{\sqrt{\sum_{j \in CP} \sigma_j^2}}\right). $$

The first assumption is the weak point. The project ends at the maximum over all paths, and $E[\max(X, Y)] \ge \max(E[X], E[Y])$. Whenever a near-critical path has high variance, it is sometimes the longest one, and PERT is optimistic about both the mean and the probability of meeting the deadline. This is the merge-event bias: at every node where paths join, delays on any incoming path propagate. Monte Carlo simulation (Chapter 12) removes the bias. It samples all durations, runs the forward pass, and repeats. It also gives the criticality index of each activity, the fraction of runs in which it is critical.

Example: analytic PERT versus simulation. The frontend estimate is very uncertain (4, 6, 16 days), so the non-critical path through D sometimes becomes the longest.

import numpy as np
from scipy.stats import norm
# Reuses project and topological_order() from the Project Networks topic.

# Three-point estimates (optimistic a, most likely m, pessimistic b) in days
est = {"A": (2, 3, 5), "B": (4, 5, 8), "C": (6, 8, 14), "D": (4, 6, 16),
       "E": (3, 4, 6), "F": (2, 3, 4), "G": (1, 2, 4), "H": (3, 4, 7)}
mean = {k: (a + 4 * m + b) / 6 for k, (a, m, b) in est.items()}
var = {k: ((b - a) / 6) ** 2 for k, (a, m, b) in est.items()}

order, _ = topological_order(project)
pred = {a: p for a, (_, p) in project.items()}

def finish_times(d):                       # forward pass, vectorized over simulation runs
    EF = {}
    for a in order:
        start = np.zeros_like(d[a])
        for p in pred[a]:
            start = np.maximum(start, EF[p])   # an activity starts when its last predecessor ends
        EF[a] = start + d[a]
    return EF

path = ["A", "B", "C", "F", "H"]           # critical path using the means
mu = sum(mean[a] for a in path)
sd = np.sqrt(sum(var[a] for a in path))
print(f"PERT: critical path mean = {mu:.2f}, sd = {sd:.2f}")
for deadline in (25, 27, 29):
    print(f"  P(T <= {deadline}) = {norm.cdf((deadline - mu) / sd):.3f}")

# Monte Carlo with Beta-PERT activity times: every path, not just one
rng = np.random.default_rng(seed=11)
n = 200_000
def beta_pert(a, m, b, size):
    alpha = 1 + 4 * (m - a) / (b - a)      # Beta-PERT shape parameters
    beta = 1 + 4 * (b - m) / (b - a)
    return a + (b - a) * rng.beta(alpha, beta, size)
d = {k: beta_pert(*est[k], n) for k in est}
EF = finish_times(d)
T = EF["H"]
print(f"simulation: mean = {T.mean():.2f}, sd = {T.std():.2f}")
for deadline in (25, 27, 29):
    print(f"  P(T <= {deadline}) = {np.mean(T <= deadline):.3f}")

# Criticality index: how often C (backend) vs D (frontend) drives the start of F
crit_C = np.mean(EF["C"] >= np.maximum(EF["D"], EF["E"]))
print(f"backend C is on the critical path in {crit_C:.0%} of runs, frontend D in {1 - crit_C:.0%}")
# PERT: critical path mean = 24.50, sd = 1.74
#   P(T <= 25) = 0.613
#   P(T <= 27) = 0.925
#   P(T <= 29) = 0.995
# simulation: mean = 25.01, sd = 1.92
#   P(T <= 25) = 0.515
#   P(T <= 27) = 0.846
#   P(T <= 29) = 0.976
# backend C is on the critical path in 72% of runs, frontend D in 28%

Analytic PERT promises a 92.5% chance of finishing within 27 days. The simulation, which lets the frontend path take over when it runs long, gives 84.6%. The frontend is critical in 28% of the runs although its mean path (23.2 days) is more than a day shorter than the backend path (24.5 days). A manager who watches only the "critical path" would ignore the riskiest activity of the project. With deadlines that matter, simulate.

Key takeaways

• PERT estimates each duration from three points: mean $(a + 4m + b)/6$ and standard deviation $(b - a)/6$.
• Summing means and variances along the critical path gives a normal approximation of the project duration and of the probability of meeting a deadline.
• PERT ignores near-critical paths, so it is optimistic (merge-event bias). The error is largest when parallel paths have similar lengths or high variance.
• Monte Carlo simulation gives unbiased completion probabilities and criticality indices for every activity.

Practice

1. (Beginner) An activity has estimates $a = 4$, $m = 7$, $b = 16$ days. Compute its PERT mean and standard deviation. Why is the mean larger than the most likely value?
2. (Intermediate) Find, with the analytic PERT formula and with the simulation, the deadline that gives a 90% chance of on-time completion. How many days of extra buffer does the merge-event bias add?
3. (Interview) "Our project-management tool reports 95% confidence for the launch date, yet our projects are almost always late. What could be wrong?" Hint: merge-event bias, correlated durations (the same team, the same optimistic estimator, shared risks), estimates that are anchored on the most likely value, and Parkinson's law (work expands to fill the time available).

Time-Cost Trade-Offs (Crashing)

Why it matters. Projects can often be accelerated by spending more money: overtime, extra staff, faster equipment, express shipping. Meanwhile every day of project duration costs overhead, and late completion may carry penalties. Crashing finds which activities to accelerate, and by how much, to reach a deadline at minimum cost, or to minimize total project cost.

Theory. Each activity $j$ has a normal duration $d_j$ and a crash duration $d_j - u_j$, the shortest possible. Shortening it by $y_j \in [0, u_j]$ days costs $c_j y_j$, where the cost slope $c_j$ is the extra cost per day saved. With start times $t_j$, the cheapest way to finish by a deadline $T$ is the LP

$$ \begin{aligned} \min \quad & \sum_j c_j, y_j \ \text{s.t.} \quad & t_i + d_i - y_i \le t_j && \forall (i \to j) \ & t_n + d_n - y_n \le T && \text{(the last activity } n\text{)} \ & 0 \le y_j \le u_j, \quad t_j \ge 0. \end{aligned} $$

Solving it for every $T$ gives the project time-cost curve. It is convex and piecewise linear: each day saved costs at least as much as the previous one. Adding the indirect cost (overhead, penalties) per day, $k,T$, and minimizing the sum gives the optimal project duration.

The classical manual procedure explains what the LP does. Repeatedly shorten the project by one day by crashing the cheapest set of activities that cuts every critical path. When two paths are critical, that set may be one activity common to both, or one activity on each path. Fulkerson (1961) and Kelley (1961) made this exact with a min-cut computation on the critical subnetwork. Each step's cost is the capacity of a minimum cut (Section 3.4). The greedy one-activity-at-a-time version can be fooled: sometimes it pays to lengthen back an activity that was crashed earlier. The LP handles this automatically.

Example: the time-cost curve of the release. Crash limits and cost slopes per day:

Activity	A	B	C	D	E	F	G	H
Max days saved $u_j$	1	2	3	2	1	1	1	1
Cost per day $c_j$	400	300	200	250	150	500	100	350

Overhead costs 450 per project day.

import numpy as np
from scipy.optimize import linprog
# Reuses project from the Project Networks topic.

acts = list(project)
n = len(acts)
idx = {a: i for i, a in enumerate(acts)}
dur = np.array([project[a][0] for a in acts], dtype=float)
max_crash = np.array([1, 2, 3, 2, 1, 1, 1, 1], dtype=float)           # days each activity can lose
slope = np.array([400, 300, 200, 250, 150, 500, 100, 350], dtype=float)  # cost per day crashed
indirect = 450.0                                                         # overhead per project day

# Variables: [t_A..t_H (start times), y_A..y_H (days crashed)]
A_ub, b_ub = [], []
for a, (_, preds) in project.items():
    for p in preds:                           # t_p + dur_p - y_p <= t_a
        row = np.zeros(2 * n)
        row[idx[p]], row[n + idx[p]], row[idx[a]] = 1, -1, -1
        A_ub.append(row); b_ub.append(-dur[idx[p]])
finish = np.zeros(2 * n)                      # t_H + dur_H - y_H <= deadline
finish[idx["H"]], finish[n + idx["H"]] = 1, -1
bounds = [(0, None)] * n + [(0, m) for m in max_crash]
c = np.r_[np.zeros(n), slope]                 # minimize the total crashing cost

print(" T  crash cost  indirect   total   crashed days")
for T in range(23, 13, -1):
    res = linprog(c, A_ub=np.array(A_ub + [finish]), b_ub=b_ub + [T - dur[idx["H"]]],
                  bounds=bounds)
    if res.status != 0:
        print(f"{T:2d}  infeasible: shortest possible duration reached")
        break
    y = res.x[n:]
    plan = " ".join(f"{a}{y[i]:.0f}" for i, a in enumerate(acts) if y[i] > 1e-6)
    print(f"{T:2d}  {res.fun:10.0f}  {indirect * T:8.0f}  {res.fun + indirect * T:6.0f}   {plan}")
#  T  crash cost  indirect   total   crashed days
# 23           0     10350   10350
# 22         200      9900   10100   C1
# 21         400      9450    9850   C2          <- path A-B-D-F-H is now critical too
# 20         700      9000    9700   B1 C2
# 19        1000      8550    9550   B2 C2
# 18        1350      8100    9450   B2 C2 H1
# 17        1750      7650    9400   A1 B2 C2 H1
# 16        2200      7200    9400   A1 B2 C3 D1 H1
# 15        2700      6750    9450   A1 B2 C3 D1 F1 H1
# 14  infeasible: shortest possible duration reached

The first two days come from the backend C, the cheapest critical activity at 200 per day. At 21 days the frontend path A-B-D-F-H becomes critical as well. Crashing C alone no longer helps, so the LP switches to the shared activities B (300), H (350) and A (400) rather than C and D together (450). Once those are used up, it crashes C and D together. Each day saved costs more than the previous one, as convexity requires. The total cost is lowest at 16 or 17 days (9,400), where the marginal crash cost (450) equals the daily overhead. Going below that loses money, and 15 days is the shortest possible duration.

Key takeaways

• Crashing shortens activities at a cost per day. The cheapest way to meet a deadline is an LP in start times and crash amounts.
• The project time-cost curve is convex and piecewise linear. Adding indirect costs gives the optimal duration, where marginal crash cost equals daily overhead.
• Every step must cut all critical paths at once, so shared activities, or combinations across parallel paths, are crashed.
• Crash only critical activities: shortening anything else costs money and saves no time.

Practice

1. (Beginner) Using the table above, which activity would you crash first to save one day, and which would you never crash, whatever the deadline? Explain using total float.
2. (Intermediate) Add the overhead to the LP as a cost on the project end time, making $T$ a variable, and solve a single LP for the optimal duration. Then add a late penalty of 1,000 per day beyond day 18 and solve again.
3. (Interview) "Adding more developers to a late software project makes it later (Brooks's law). How does that fit the crashing model?" Hint: the model assumes linear, independent crash costs. Think about onboarding time, communication overhead that grows with team size, and whether the activity can actually be split.

Resource-Constrained Project Scheduling

Why it matters. CPM assumes unlimited resources: any number of activities can run in parallel. In reality, the same four developers, one crane or one test lab must be shared. The CPM schedule then asks for more people than exist on some days, and the true duration is longer. The resource-constrained project scheduling problem (RCPSP) is the most studied problem in project scheduling. Its heuristics are inside every major project-management tool's "level resources" button.

Theory. Each activity $j$ needs $r_{jk}$ units of renewable resource $k$ (people, machines) while it runs, and resource $k$ has a capacity of $R_k$ units per period. The RCPSP chooses start times $S_j$ that respect precedence and, at every time $t$,

$$ \sum_{j ,:, S_j \le t < S_j + d_j} r_{jk} \le R_k \qquad \forall k, $$

while minimizing the makespan. The RCPSP is strongly NP-hard: the job shop problem (Section 10.3) is the special case where every resource is a machine of capacity 1. Exact methods (branch and bound, time-indexed MIPs with variables $x_{jt} = 1$ if $j$ starts at $t$, constraint programming) solve instances with 30 activities routinely and many with 60. Larger instances use heuristics.

The workhorse heuristic is the schedule generation scheme (SGS) driven by a priority rule:

• Serial SGS. Take the activities one at a time, in a precedence-feasible priority list, and start each at the earliest time that respects its predecessors and leaves enough resources. It produces active schedules, in which no activity can start earlier without delaying another. For regular objectives such as the makespan, some active schedule is optimal, and every active schedule comes from some list. Searching over lists is therefore enough, which is what genetic algorithms and other metaheuristics for the RCPSP do (Chapter 11).
• Parallel SGS. Move forward in time and, at each decision point, start as many eligible activities as resources allow, in priority order. It produces non-delay schedules, which never leave resources idle when some activity could use them. Non-delay schedules may exclude every optimal schedule.

Common priority rules:

Rule	Priority to the activity with	Idea
LFT	smallest CPM latest finish time	urgency with respect to the deadline
MinSlack	smallest total float	protect the critical path
SPT	shortest duration	finish many activities quickly
GRD	greatest resource demand ($d_j \times r_j$)	place the hardest activities first

Resource leveling is a related problem: the deadline is fixed, and the goal is to smooth the resource profile (fewer peaks, less hiring and firing) by moving activities within their float.

Example: four developers. The activities need A 2, B 3, C 2, D 2, E 2, F 3, G 1 and H 2 developers, and only 4 are available. The CPM schedule needs 5 people on days 3 to 6, when the design B and the database setup E run together. The code runs the serial SGS with four priority rules, and finds the optimum by trying every precedence-feasible list (there are 23).

from itertools import permutations
# Reuses project from the Project Networks topic.

demand = {"A": 2, "B": 3, "C": 2, "D": 2, "E": 2, "F": 3, "G": 1, "H": 2}   # developers needed
CAP = 4                                                                      # developers available
dur = {a: d for a, (d, _) in project.items()}
pred = {a: p for a, (_, p) in project.items()}

def serial_sgs(activity_list):
    """Schedule each activity, in list order, at the earliest precedence- and resource-feasible time."""
    usage = [0] * 100                                   # developers in use on each day
    start = {}
    for a in activity_list:
        t = max((start[p] + dur[p] for p in pred[a]), default=0)
        while any(usage[u] + demand[a] > CAP for u in range(t, t + dur[a])):
            t += 1                                      # shift right until it fits
        for u in range(t, t + dur[a]):
            usage[u] += demand[a]
        start[a] = t
    return start, max(start[a] + dur[a] for a in start)

def priority_list(key):
    """Precedence-feasible list: among eligible activities, pick the one with the smallest key."""
    done, lst = set(), []
    while len(lst) < len(project):
        eligible = [a for a in project if a not in done and all(p in done for p in pred[a])]
        a = min(eligible, key=key)
        lst.append(a); done.add(a)
    return lst

# CPM latest finish times and total floats (from the CPM topic)
LF = {"A": 3, "B": 8, "C": 16, "D": 16, "E": 16, "F": 19, "G": 19, "H": 23}
TF = {"A": 0, "B": 0, "C": 0, "D": 2, "E": 9, "F": 0, "G": 3, "H": 0}
rules = {"LFT (latest finish)": lambda a: LF[a],
         "MinSlack":            lambda a: TF[a],
         "SPT (shortest first)": lambda a: dur[a],
         "GRD (most work first)": lambda a: -dur[a] * demand[a]}
for name, key in rules.items():
    lst = priority_list(key)
    start, ms = serial_sgs(lst)
    print(f"{name:22} list {''.join(lst)} -> makespan {ms}")

# Exact optimum for this tiny instance: some activity list always yields an optimal schedule
best = min((serial_sgs(p) for p in permutations(project)
            if all(p.index(q) < p.index(a) for a in p for q in pred[a])), key=lambda r: r[1])
print("optimal makespan:", best[1], "starts:", best[0])

start, ms = serial_sgs(priority_list(rules["LFT (latest finish)"]))
for a in sorted(start, key=start.get):                 # text Gantt chart of the LFT schedule
    print(f"{a} ({demand[a]}) |" + " " * start[a] + "#" * dur[a])
print("devs used |" + "".join(str(sum(demand[a] for a in start if start[a] <= t < start[a] + dur[a]))
                              for t in range(ms)))
# LFT (latest finish)    list ABCDEFGH -> makespan 25
# MinSlack               list ABCDGEFH -> makespan 27
# SPT (shortest first)   list AEBDGCFH -> makespan 27
# GRD (most work first)  list ABCDEFGH -> makespan 25
# optimal makespan: 25 starts: {'A': 0, 'B': 3, 'C': 8, 'D': 8, 'E': 14, 'F': 18, 'G': 16, 'H': 21}
# A (2) |###
# B (3) |   #####
# C (2) |        ########
# D (2) |        ######
# E (2) |              ####
# G (1) |                ##
# F (3) |                  ###
# H (2) |                     ####
# devs used |2223333344444444333332222

With four developers the release takes 25 days instead of 23. The database setup E cannot run beside the design B (3 + 2 > 4) or beside C and D together, so it waits until the frontend finishes on day 14. That is 2 days past its CPM latest start, so the project slips by 2 days. LFT and GRD find the optimum. MinSlack schedules the documentation G before E, so E has to wait until day 16. SPT starts the short activity E first and pushes the critical design B back by 4 days. Both lose 2 more days. No single rule is best on every instance, which is why practical tools try several rules or improve the list with local search. Note also that the critical path has changed: E, which had 9 days of float under CPM, is now critical. With scarce resources, criticality comes from resource conflicts as well as precedence (the critical chain idea).

Key takeaways

• The RCPSP adds resource capacities to CPM. It is strongly NP-hard and generalizes the job shop problem.
• The serial schedule generation scheme turns any priority list into an active schedule, and some list always gives an optimal one.
• Priority rules (LFT, MinSlack, SPT, GRD) are fast but none dominates. Use several, or improve lists with metaheuristics.
• Resource conflicts can make activities with large CPM float critical. Float computed without resources is optimistic.

Practice

1. (Beginner) Draw the developer profile of the CPM schedule of the example day by day, and mark the days where it exceeds 4.
2. (Intermediate) Implement the parallel SGS: at each time $t$, start eligible activities in priority order while resources allow, then move $t$ to the next activity completion. Compare its makespans with the serial SGS for the four rules.
3. (Interview) "Your project tool's resource leveling pushed the end date out by three weeks. The sponsor wants it back. What are your options?" Hint: more resources (and which ones: the bottleneck), shortening resource-heavy activities, changing the logic (overlaps), changing who does what (resource flexibility), and checking whether the leveling heuristic itself is far from optimal.

<a id="103-machine-and-job-scheduling"></a>

10.3 Machine and Job Scheduling

Project scheduling deals with one project and its activities. Machine scheduling deals with many jobs competing for a few machines: orders in a factory, patients in an operating theater, tasks on CPU cores, aircraft on a runway. This subchapter introduces the standard classification of these problems and solves the single-machine case with simple sorting rules. It then moves to several machines in series (flow shops and Johnson's algorithm), jobs with individual routes (job shops), and identical machines in parallel. A recurring theme is the boundary between problems that a sorting rule solves exactly and NP-hard ones where heuristics must be used (Section 4.4).

Single-Machine Scheduling

Why it matters. A single machine is the simplest scheduling setting, and it is also common in practice: one bottleneck machine, one surgeon, one CPU core, one loading dock. Its results are the building blocks for harder problems, because complex shops are often scheduled by decomposing them into single-machine subproblems at the bottleneck. It also shows the most important proof technique in scheduling, the exchange argument.

Theory. Scheduling problems are classified with Graham's three-field notation $\alpha \mid \beta \mid \gamma$:

Field	Meaning	Examples
$\alpha$: machine environment	how many machines and how they are arranged	$1$ single machine, $Pm$ identical parallel, $Qm$ uniform, $Rm$ unrelated, $Fm$ flow shop, $Jm$ job shop
$\beta$: job characteristics	constraints on the jobs	$r_j$ release dates, $prec$ precedence, $pmtn$ preemption allowed, $s_{jk}$ setup times, $prmu$ permutation
$\gamma$: objective	what is minimized	$C_{\max}$ makespan, $\sum w_j C_j$, $L_{\max}$, $\sum T_j$, $\sum U_j$

Each job $j$ has a processing time $p_j$, and possibly a weight $w_j$ and a due date $d_j$. A schedule determines the completion time $C_j$, from which follow the lateness $L_j = C_j - d_j$ (negative when early), the tardiness $T_j = \max(L_j, 0)$ and the tardy indicator $U_j = 1$ if $T_j > 0$.

On one machine without release dates, no idle time is useful and $C_{\max} = \sum_j p_j$ for every sequence. The interesting objectives depend on the order. The weighted sum of completion times $\sum w_j C_j$ measures weighted average throughput time, or inventory cost when $w_j$ is a holding cost. Smith's rule (1956), also called WSPT (weighted shortest processing time), is optimal for $1 ,||, \sum w_j C_j$: sequence jobs in increasing order of $p_j / w_j$.

Exchange argument. Take any sequence in which job $i$ is immediately followed by job $j$, and swap them. Only the completion times of $i$ and $j$ change. If $t$ is the time at which the pair starts, then

$$ \text{($i$ first)}: \ w_i (t + p_i) + w_j (t + p_i + p_j), \qquad \text{($j$ first)}: \ w_j (t + p_j) + w_i (t + p_j + p_i). $$

The difference ($i$ first minus $j$ first) is $p_i w_j - p_j w_i$. Putting $i$ first is better exactly when $p_i / w_i < p_j / w_j$. Any sequence that violates WSPT has an adjacent pair that can be swapped to improve it, so only the WSPT order can be optimal. The analogy is a supermarket queue that lets customers with few items, or in a hurry, go first: it minimizes total waiting.

Small changes to the problem change its complexity dramatically:

Problem	Complexity	Solution
$1 ,\vert\vert, \sum w_j C_j$	polynomial	WSPT (Smith's rule)
$1 ,\vert\vert, L_{\max}$	polynomial	EDD (next topic)
$1 ,\vert\vert, \sum U_j$	polynomial	Moore-Hodgson (next topic)
$1 ,\vert, prec ,\vert, L_{\max}$	polynomial	Lawler's algorithm (schedule backwards)
$1 ,\vert, r_j, pmtn ,\vert, \sum C_j$	polynomial	SRPT (shortest remaining processing time)
$1 ,\vert\vert, \sum T_j$	NP-hard (ordinary sense)	dynamic programming, pseudo-polynomial
$1 ,\vert\vert, \sum w_j T_j$	strongly NP-hard	branch and bound, heuristics (e.g. ATC rule)
$1 ,\vert, r_j ,\vert, L_{\max}$	strongly NP-hard	branch and bound

Example: WSPT versus brute force. Six jobs with processing times $p$, weights $w$ and due dates $d$:

from itertools import permutations

# job: (processing time p, weight w, due date d)
jobs = {1: (4, 2, 9), 2: (2, 1, 6), 3: (6, 4, 14), 4: (3, 3, 8), 5: (5, 1, 24), 6: (1, 2, 3)}

def evaluate(seq):
    t, wC, wT = 0, 0, 0
    for j in seq:
        p, w, d = jobs[j]
        t += p                                # completion time C_j
        wC += w * t
        wT += w * max(t - d, 0)               # weighted tardiness
    return wC, wT

wspt = sorted(jobs, key=lambda j: jobs[j][0] / jobs[j][1])   # Smith's rule: increasing p/w
spt = sorted(jobs, key=lambda j: jobs[j][0])
print("WSPT", wspt, "sum wC = %d, sum wT = %d" % evaluate(wspt))
print("SPT ", spt, "sum wC = %d, sum wT = %d" % evaluate(spt))

all_seqs = list(permutations(jobs))                          # 720 sequences
best_wC = min(all_seqs, key=lambda s: evaluate(s)[0])
best_wT = min(all_seqs, key=lambda s: evaluate(s)[1])
print("best for sum wC:", list(best_wC), evaluate(best_wC)[0])
print("best for sum wT:", list(best_wT), evaluate(best_wT)[1])

# Exchange argument: swapping adjacent jobs i, j changes sum wC by p_i w_j - p_j w_i
i, j = 5, 3                                     # a long, light job placed before a heavy one
(pi, wi, _), (pj, wj, _) = jobs[i], jobs[j]
print(f"i={i} before j={j} costs {pi * wj - pj * wi:+d} versus j before i")
# WSPT [6, 4, 3, 1, 2, 5] sum wC = 119, sum wT = 20
# SPT  [6, 2, 4, 1, 5, 3] sum wC = 142, sum wT = 30
# best for sum wC: [6, 4, 3, 1, 2, 5] 119      <- WSPT is optimal
# best for sum wT: [2, 6, 4, 1, 3, 5] 10       <- but not for weighted tardiness
# i=5 before j=3 costs +14 versus j before i

WSPT matches the best of all 720 sequences for $\sum w_j C_j$. Ignoring weights (plain SPT) costs 19% more, because SPT runs the heavy job 3 last. For weighted tardiness WSPT is twice the optimum, since it ignores due dates. That objective is strongly NP-hard, and practical schedulers use rules that combine both, such as the apparent tardiness cost (ATC) rule, or local search.

Key takeaways

• Graham's notation $\alpha \mid \beta \mid \gamma$ classifies scheduling problems by machines, job constraints and objective.
• On one machine, WSPT (increasing $p_j / w_j$) minimizes the weighted sum of completion times.
• The exchange argument proves optimality of sorting rules: any violating sequence has an adjacent pair whose swap helps.
• Small changes (release dates, tardiness, weights) turn easy problems into NP-hard ones. Check the complexity before searching for an exact rule.

Practice

1. (Beginner) Four jobs have $(p_j, w_j)$ = (3, 1), (5, 4), (2, 2), (6, 3). Find the WSPT sequence and its $\sum w_j C_j$.
2. (Intermediate) Prove with an exchange argument that SPT minimizes $\sum C_j$ on one machine. Then show that it also minimizes the average number of jobs waiting (WIP) over the schedule, using Little's law (Section 9.3).
3. (Interview) "An operating system scheduler always runs the shortest job first. What could go wrong?" Hint: processing times are not known in advance, new short jobs keep arriving (starvation of long jobs), and fairness or deadlines matter. Look up how real schedulers approximate SRPT.

Shortest Processing Time and EDD Rules

Why it matters. Dispatching rules, which decide "which job next?" by sorting, run most shop floors, help desks and computing clusters. They need no optimization software and can be applied in real time. Two rules matter most: SPT (shortest processing time first) for speed and EDD (earliest due date first) for punctuality. Knowing exactly what each one optimizes, and what it sacrifices, is basic scheduling literacy.

Theory. On one machine, with all jobs available at time 0:

• SPT minimizes $\sum C_j$, the mean flow time. Since $\sum L_j = \sum C_j - \sum d_j$, it also minimizes mean lateness. By Little's law, it minimizes the average number of jobs in the system (work in process). The exchange argument of the previous topic, with $w_j = 1$, proves it.
• EDD (Jackson's rule, 1955) minimizes the maximum lateness $L_{\max}$, and therefore also the maximum tardiness $T_{\max}$. Exchange argument: if $d_i > d_j$ and $i$ is right before $j$, swapping them does not increase the larger of their two latenesses, because after the swap $j$ finishes earlier, and $i$ finishes when $j$ used to finish but has a later due date. If some sequence meets every due date, EDD finds one.
• Moore-Hodgson (1968) minimizes the number of tardy jobs $\sum U_j$. Add jobs in EDD order. Whenever the current job is late, remove the longest job scheduled so far and put it at the end. Removing the longest job frees the most time.

Each rule optimizes one criterion and can be poor on the others:

Rule	Optimal for	Weakness
SPT	mean flow time, mean lateness, mean WIP	long jobs wait indefinitely; ignores due dates
EDD	maximum lateness and tardiness	can make many jobs slightly late; ignores job length
Moore-Hodgson	number of tardy jobs	the jobs it gives up on can be very late
FIFO (arrival order)	nothing in particular	perceived as fair; easy to explain

Minimizing total tardiness $\sum T_j$ is NP-hard. SPT works well when due dates are loose (everything will be late anyway, so finish many jobs fast), and EDD when they are tight but achievable. Composite rules, such as the modified due date (MDD) rule, which sorts by $\max(d_j, t + p_j)$, interpolate between the two.

Example: four rules on six jobs. The jobs of the previous topic, without weights.

jobs = {1: (4, 9), 2: (2, 6), 3: (6, 14), 4: (3, 8), 5: (5, 24), 6: (1, 3)}   # (p, due date)

def metrics(seq):
    t, C, late = 0, [], []
    for j in seq:
        t += jobs[j][0]
        C.append(t)
        late.append(t - jobs[j][1])            # lateness L_j = C_j - d_j
    tardy = [max(x, 0) for x in late]
    return sum(C) / len(C), max(late), sum(x > 0 for x in late), sum(tardy)

def moore_hodgson(jobs):
    """Minimize the number of tardy jobs: EDD, and when a job is late drop the longest so far."""
    on_time, t = [], 0
    for j in sorted(jobs, key=lambda j: jobs[j][1]):
        on_time.append(j)
        t += jobs[j][0]
        if t > jobs[j][1]:                     # current job is late
            longest = max(on_time, key=lambda k: jobs[k][0])
            on_time.remove(longest)
            t -= jobs[longest][0]
    late = [j for j in jobs if j not in on_time]
    return on_time + late                      # late jobs go at the end, in any order

rules = {"SPT": sorted(jobs, key=lambda j: jobs[j][0]),
         "EDD": sorted(jobs, key=lambda j: jobs[j][1]),
         "Moore-Hodgson": moore_hodgson(jobs),
         "FIFO (job index)": sorted(jobs)}
print(f"{'rule':>16}  sequence             mean flow  Lmax  #tardy  sum T")
for name, seq in rules.items():
    F, Lmax, nT, sT = metrics(seq)
    print(f"{name:>16}  {str(seq):<20} {F:9.2f} {Lmax:5d} {nT:7d} {sT:6d}")
#             rule  sequence             mean flow  Lmax  #tardy  sum T
#              SPT  [6, 2, 4, 1, 5, 3]        9.33     7       2      8
#              EDD  [6, 2, 4, 1, 3, 5]        9.50     2       2      3
#    Moore-Hodgson  [6, 2, 4, 3, 5, 1]       10.00    12       1     12
# FIFO (job index)  [1, 2, 3, 4, 5, 6]       13.00    18       2     25

Each rule wins on its own criterion: SPT has the lowest mean flow time, EDD the lowest maximum lateness, and Moore-Hodgson the fewest tardy jobs. Moore-Hodgson gets there by giving up on job 1, which ends 12 days late. If job 1 belongs to the most important customer, the "optimal" schedule is a bad one. EDD also happens to give the minimum total tardiness here (3, confirmed by brute force), but that is not guaranteed in general. FIFO, the default in many organizations, is worse than every other rule on every measure.

Key takeaways

• SPT minimizes mean flow time, mean lateness and mean work in process.
• EDD minimizes maximum lateness and tardiness, and finds a schedule with no late jobs whenever one exists.
• Moore-Hodgson minimizes the number of tardy jobs, possibly at the cost of a few very late ones.
• Total tardiness is NP-hard. Choose a rule according to the criterion that matters, and the tightness of the due dates.

Practice

1. (Beginner) Five jobs have $(p_j, d_j)$ = (4, 6), (2, 4), (6, 18), (3, 8), (5, 12). Compute the SPT and EDD sequences and their mean flow time and $L_{\max}$.
2. (Intermediate) Generate 1,000 random instances of 8 jobs, with due dates drawn uniformly over a range proportional to $\sum p_j$. Compare SPT and EDD on total tardiness as the due dates go from tight to loose. Where does each rule win?
3. (Interview) "A support team sorts tickets by deadline (EDD). Most tickets are slightly late. The manager wants a rule that keeps most tickets on time. What do you suggest, and what is the risk?" Hint: Moore-Hodgson maximizes on-time tickets but abandons some. Discuss weights, an upper limit on lateness, and how customers perceive abandoned tickets.

Flow Shop Scheduling

Why it matters. Many production lines process every job through the same machines in the same order: cut, then weld, then paint; or mix, then bake, then pack. This is a flow shop. The order in which jobs enter the line determines how much each machine sits idle waiting for the previous one, and hence the makespan. Flow shop scheduling is a classic testbed for heuristics, and NEH, the best-known constructive heuristic, is used directly in industry.

Theory. In $Fm ,||, C_{\max}$, each of $n$ jobs visits machines $1, 2, \dots, m$ in order, with processing time $p_{jk}$ on machine $k$. A machine works on one job at a time. In a permutation schedule all machines process the jobs in the same order. Permutation schedules are optimal for $m \le 3$, and they are the standard assumption ($Fm ,|, prmu ,|, C_{\max}$) because they are simple to run. For a sequence $\pi$, completion times follow the recursion

$$ C_{\pi(i), k} = \max\big(C_{\pi(i-1), k},\ C_{\pi(i), k-1}\big) + p_{\pi(i), k}: $$

job $\pi(i)$ starts on machine $k$ when the machine has finished the previous job and the job has left the previous machine. The makespan is $C_{\pi(n), m}$. Evaluating a sequence takes $O(nm)$ time. The recursion is again a longest path, through a grid:

            machine 1      machine 2      machine 3
 job pi(1)  [p11] ------> [p12] ------> [p13]
              |             |             |
              v             v             v
 job pi(2)  [p21] ------> [p22] ------> [p23]
              |             |             |
              v             v             v
 job pi(3)  [p31] ------> [p32] ------> [p33]      makespan = longest path from top-left to bottom-right

$F2 ,||, C_{\max}$ is solved exactly by Johnson's algorithm (next topic), but $F3 ,||, C_{\max}$ is already strongly NP-hard. With $n!$ sequences, enumeration stops being possible at about 12 jobs, so heuristics are used:

• Palmer's slope index and CDS (Campbell-Dudek-Smith, which applies Johnson's rule to aggregated two-machine problems): fast, moderate quality.
• NEH (Nawaz, Enscore and Ham, 1983): sort jobs by decreasing total processing time, then insert them one at a time into the partial sequence, each at the position that minimizes the partial makespan. With Taillard's acceleration it runs in $O(n^2 m)$ time, and it is usually within a few percent of the optimum.
• Metaheuristics (iterated greedy, tabu search, Chapter 11) start from NEH and improve it.

A lower bound shows how good a sequence is. Machine $k$ must process all its work, and before it starts it must wait at least for the shortest "head" (work of a job on machines $1, \dots, k-1$); after it finishes, the last job still needs at least the shortest "tail":

$$ LB = \max_k \Big( \min_j \sum_{k' < k} p_{jk'} + \sum_j p_{jk} + \min_j \sum_{k' > k} p_{jk'} \Big). $$

(The minimum head and minimum tail can come from different jobs, which is why this is only a bound.) Another bound is the longest job's total processing time.

Example: NEH versus all 720 sequences.

import numpy as np
from itertools import permutations

# p[j][k]: processing time of job j on machine k (3 machines, 6 jobs)
p = np.array([[5, 9, 8], [9, 3, 10], [9, 4, 5], [4, 8, 8], [3, 5, 6], [6, 7, 2]])

def makespan(seq):
    """Permutation flow shop: C[i][k] = max(C[i-1][k], C[i][k-1]) + p."""
    C = np.zeros(p.shape[1])                 # completion of the previous job on each machine
    for j in seq:
        for k in range(p.shape[1]):
            prev = C[k - 1] if k > 0 else 0  # this job finished on machine k-1 (already updated)
            C[k] = max(C[k], prev) + p[j, k]
    return int(C[-1])

def neh(p):
    """NEH: sort by total work, insert each job at the best position of the partial sequence."""
    order = sorted(range(len(p)), key=lambda j: -p[j].sum())
    seq = [order[0]]
    for j in order[1:]:
        candidates = [seq[:i] + [j] + seq[i:] for i in range(len(seq) + 1)]
        seq = min(candidates, key=makespan)
    return seq

seq = neh(p)
print("NEH sequence:", [j + 1 for j in seq], "makespan", makespan(seq))
spans = [makespan(s) for s in permutations(range(len(p)))]
print(f"all 720 sequences: best {min(spans)}, mean {np.mean(spans):.1f}, worst {max(spans)}")
lb = max(max(p[:, k].sum() + min(p[:, :k].sum(axis=1), default=0)
             + min(p[:, k + 1:].sum(axis=1), default=0) for k in range(3)), p.sum(axis=1).max())
print("machine-based lower bound:", lb)
# NEH sequence: [5, 4, 2, 1, 6, 3] makespan 49
# all 720 sequences: best 48, mean 55.6, worst 63
# machine-based lower bound: 47

NEH evaluates only 20 partial sequences, yet its makespan of 49 is one unit (2%) above the best of the 720 sequences. A random sequence averages 55.6 and the worst reaches 63, so the sequence matters. NEH starts with job 5, which is short on the first machine and gets the line going quickly. The lower bound of 47 certifies, without enumerating anything, that no sequence can beat NEH by more than 2 units (about 4%).

Key takeaways

• In a flow shop all jobs visit the machines in the same order. Permutation schedules are standard, and their makespan follows a simple max-plus recursion.
• $F2 ,||, C_{\max}$ is polynomial, but three or more machines make the problem strongly NP-hard.
• NEH (sort by total work, insert each job at its best position) is the standard constructive heuristic and is usually within a few percent of optimal.
• Lower bounds based on machine loads measure the quality of a heuristic solution when the optimum is unknown.

Practice

1. (Beginner) Compute by hand the makespan of the sequence 1-2-3-4-5-6 for the example data, filling in the completion-time table row by row.
2. (Intermediate) Add a simple local search after NEH: try every swap of two jobs and every reinsertion of one job at another position, and accept improvements until none is left. Does it reach 48?
3. (Interview) "The heuristic says this sequence is 2% from optimal. How can you claim that without knowing the optimum?" Hint: compare the makespan with a lower bound. Explain how the machine-based bound is built and why it can be loose.

Johnson's Algorithm

Why it matters. The two-machine flow shop, such as printing then binding, or machining then inspection, is one of the few scheduling problems with a simple, exact, fast solution. Johnson's rule (1954) is the first optimal scheduling algorithm in the OR literature. It is still used directly, it serves as a subroutine in heuristics for longer lines (CDS), and its proof shows the kind of reasoning that leads to exact sorting rules.

Theory. Each job $j$ needs $a_j$ time on machine 1, then $b_j$ on machine 2. Machine 1 never needs to be idle. Machine 2 is idle whenever it waits for machine 1, and the makespan is the total work on machine 2 plus that idle time. So the goal is to keep machine 2 fed. Johnson's rule:

1. Put jobs with $a_j < b_j$ in a first set, sorted by increasing $a_j$.
2. Put jobs with $a_j \ge b_j$ in a second set, sorted by decreasing $b_j$.
3. Sequence: the first set, then the second set.

The intuition: start with jobs that are short on machine 1, so machine 2 gets work early. End with jobs that are short on machine 2, so that little is left on machine 2 after machine 1 finishes. Equivalently, job $i$ goes before job $j$ if $\min(a_i, b_j) \le \min(a_j, b_i)$, and an exchange argument proves that this order is optimal. The rule runs in $O(n \log n)$ time.

A quick lower bound checks the result: machine 1 is busy for $\sum a_j$, and afterwards the last job needs at least $\min_j b_j$ on machine 2; symmetrically, machine 2 cannot start before $\min_j a_j$:

$$ C_{\max} \ge \max\Big(\sum_j a_j + \min_j b_j,\ \ \min_j a_j + \sum_j b_j\Big). $$

Extension to three machines. If machine 2 is dominated, meaning $\min_j p_{j1} \ge \max_j p_{j2}$ or $\min_j p_{j3} \ge \max_j p_{j2}$, then Johnson's rule applied to the virtual times $a_j = p_{j1} + p_{j2}$ and $b_j = p_{j2} + p_{j3}$ is optimal for the three-machine problem. Otherwise it is a good heuristic.

Example: six jobs on two machines.

from itertools import permutations

jobs = {"J1": (3, 6), "J2": (7, 2), "J3": (4, 4), "J4": (5, 3), "J5": (1, 5), "J6": (6, 6)}  # (M1, M2)

def johnson(jobs):
    first = sorted((j for j, (a, b) in jobs.items() if a < b), key=lambda j: jobs[j][0])
    last = sorted((j for j, (a, b) in jobs.items() if a >= b), key=lambda j: -jobs[j][1])
    return first + last                  # short M1 jobs early, short M2 jobs late

def schedule(seq):
    t1 = t2 = 0
    out = []
    for j in seq:
        a, b = jobs[j]
        s1, t1 = t1, t1 + a              # M1 works without idle time
        s2 = max(t1, t2)                 # M2 waits for the job and for itself
        t2 = s2 + b
        out.append((j, s1, s2))
    return out, t2

seq = johnson(jobs)
plan, cmax = schedule(seq)
print("Johnson:", " ".join(seq), "makespan", cmax)
print("brute force best:", min(schedule(s)[1] for s in permutations(jobs)))
print("index order makespan:", schedule(list(jobs))[1])
for row, k in (("M1", 0), ("M2", 1)):    # text Gantt chart, one character per time unit
    line = [" "] * cmax
    for j, s1, s2 in plan:
        s = (s1, s2)[k]
        for t in range(s, s + jobs[j][k]):
            line[t] = j[1]
    print(f"{row} |{''.join(line)}|")
# Johnson: J5 J1 J6 J3 J4 J2 makespan 28
# brute force best: 28
# index order makespan: 33
# M1 |51116666663333444442222222  |
# M2 | 555551111116666663333444 22|

Johnson's sequence reaches the optimum of 28, against 33 for the jobs in index order. It also meets the lower bound: both machines have 26 units of work, and $\max(26 + 2, 1 + 26) = 28$, which proves optimality without the brute-force check. The Gantt chart shows how the rule works. Job 5 (1 unit on M1) starts M2 almost at once, and M2 is idle only twice: for 1 unit at the start, and for 1 unit before the last job, J2, chosen because it is the shortest on M2.

Key takeaways

• Johnson's rule solves $F2 ,||, C_{\max}$ exactly in $O(n \log n)$: jobs with $a_j < b_j$ first by increasing $a_j$, then the rest by decreasing $b_j$.
• The idea is to feed machine 2 early and leave it little work at the end.
• A simple machine-based lower bound often proves the schedule optimal.
• With a dominated middle machine, the rule extends exactly to three machines. Otherwise it gives good heuristics (CDS) for longer lines.

Practice

1. (Beginner) Apply Johnson's rule to five jobs with $(a_j, b_j)$ = (5, 2), (1, 6), (9, 7), (3, 8), (10, 4). Draw the Gantt chart and give the makespan.
2. (Intermediate) Prove the exchange step: for adjacent jobs $i, j$ with $\min(a_i, b_j) \le \min(a_j, b_i)$, show that putting $i$ before $j$ does not increase the completion time of the pair on machine 2, whatever the starting times of the two machines.
3. (Interview) "A print shop has a printer and a binder. Jobs arrive during the day, not all at once. Does Johnson's rule still apply?" Hint: with release dates, $F2 ,|, r_j ,|, C_{\max}$ is strongly NP-hard. Discuss re-applying the rule to the waiting jobs whenever one arrives, and how you would test whether that works.

Job Shop Scheduling

Why it matters. In a job shop, each job has its own route through the machines: a machine shop where one part goes lathe, mill, drill and another goes drill, lathe. Job shops are the most general and hardest classical scheduling problems. The famous 10-job, 10-machine instance of Fisher and Thompson (1963) resisted every attempt for 26 years before Carlier and Pinson solved it in 1989. Job shop ideas reach well beyond factories: hospital patients routed through tests, and tasks routed through specialized processors.

Theory. In $Jm ,||, C_{\max}$, each job is a fixed sequence of operations, each on a given machine with a given time. A machine handles one operation at a time, and the operations of a job are done in order. The problem is strongly NP-hard, and only small special cases are polynomial (for example $J2 ,||, C_{\max}$ with at most two operations per job, solved by an extension of Johnson's rule due to Jackson).

The disjunctive graph (Roy and Sussmann, 1964) represents the problem:

• one node per operation, plus a source and a sink,
• conjunctive arcs for each job's route (fixed),
• disjunctive edges between every pair of operations on the same machine, which must be given a direction (which one goes first).

          J1: [M1 3] --> [M2 2] --> [M3 2]
         /                                \
 source -- J2: [M1 2] --> [M3 1] --> [M2 4] -- sink
         \ J3: [M2 4] --> [M3 3] --> [M1 1] /
          \J4: [M3 2] --> [M1 3] --> [M2 3]/

 disjunctive pairs: the four operations on M1 must be put in an order, and so must those on M2 and M3

Orienting every disjunctive edge corresponds to choosing a sequence on each machine. The choice is feasible exactly when the resulting graph has no cycle, and its makespan is then the longest source-to-sink path. Every exact method and most heuristics work on this graph:

• MIP with big-M disjunctions: $S_{j} \ge S_{i} + p_i - M (1 - y_{ij})$ and $S_{i} \ge S_{j} + p_j - M y_{ij}$ for each pair on the same machine. Easy to write but weak in LP relaxation, so constraint programming solvers usually do better.
• Branch and bound on the orientation of disjunctive edges, with strong lower bounds from one-machine relaxations.
• Shifting bottleneck (Adams, Balas and Zawack, 1988): repeatedly schedule the most critical machine by solving a $1 ,|, r_j ,|, L_{\max}$ subproblem, then re-optimize the machines already scheduled.
• Dispatching rules with the Giffler-Thompson algorithm (1960). It builds an active schedule: find the operation that can finish earliest, look at all operations that compete with it on the same machine and could start before that time (the conflict set), and let a priority rule pick one. Some active schedule is optimal, so the procedure loses nothing in principle; the quality depends on the rule.
• Metaheuristics (tabu search on critical-path moves, Chapter 11).

Example: dispatching rules versus the exact optimum. The four jobs of the diagram, with machines numbered from 0 in the code. The optimum is found by trying all $(4!)^3 = 13{,}824$ combinations of machine sequences, which is possible only because the instance is tiny.

from itertools import permutations, product

# Each job is a route: list of (machine, processing time), visited in this order
jobs = [[(0, 3), (1, 2), (2, 2)],
        [(0, 2), (2, 1), (1, 4)],
        [(1, 4), (2, 3), (0, 1)],
        [(2, 2), (0, 3), (1, 3)]]
M = 3

def giffler_thompson(rule):
    """Active schedule: repeatedly resolve the conflict on the machine that can finish first."""
    nxt = [0] * len(jobs)                        # index of each job's next operation
    job_ready, mach_ready = [0] * len(jobs), [0] * M
    while any(nxt[j] < len(jobs[j]) for j in range(len(jobs))):
        cand = [j for j in range(len(jobs)) if nxt[j] < len(jobs[j])]
        def est(j):                              # earliest start of job j's next operation
            return max(job_ready[j], mach_ready[jobs[j][nxt[j]][0]])
        jstar = min(cand, key=lambda j: est(j) + jobs[j][nxt[j]][1])
        m = jobs[jstar][nxt[jstar]][0]
        cstar = est(jstar) + jobs[jstar][nxt[jstar]][1]   # earliest completion time overall
        conflict = [j for j in cand if jobs[j][nxt[j]][0] == m and est(j) < cstar]
        j = min(conflict, key=lambda j: rule(j, nxt[j]))  # priority rule breaks the conflict
        start = est(j)
        job_ready[j] = mach_ready[m] = start + jobs[j][nxt[j]][1]
        nxt[j] += 1
    return max(job_ready)

def remaining(j, k):
    return sum(t for _, t in jobs[j][k:])
rules = {"SPT":  lambda j, k: jobs[j][k][1],                     # shortest operation
         "LPT":  lambda j, k: -jobs[j][k][1],
         "MWKR": lambda j, k: -remaining(j, k),                  # most work remaining
         "LWKR": lambda j, k: remaining(j, k)}
for name, rule in rules.items():
    print(f"{name:5} makespan {giffler_thompson(rule)}")

def semi_active(orders):
    """Makespan for given machine sequences, or None if they create a cycle."""
    nxt, pos = [0] * len(jobs), [0] * M
    job_ready, mach_ready = [0] * len(jobs), [0] * M
    for _ in range(sum(map(len, jobs))):
        for m in range(M):                       # find an operation first in both queues
            if pos[m] < len(orders[m]):
                j = orders[m][pos[m]]
                if nxt[j] < len(jobs[j]) and jobs[j][nxt[j]][0] == m:
                    t = max(job_ready[j], mach_ready[m]) + jobs[j][nxt[j]][1]
                    job_ready[j] = mach_ready[m] = t
                    nxt[j] += 1; pos[m] += 1
                    break
        else:
            return None                          # deadlock: the sequences are inconsistent
    return max(job_ready)

spans = [semi_active(o) for o in product(permutations(range(len(jobs))), repeat=M)]
feasible = [s for s in spans if s is not None]
print(f"{len(spans)} machine orderings, {len(feasible)} feasible, optimum {min(feasible)}")
loads = [sum(t for route in jobs for m2, t in route if m2 == m) for m in range(M)]
print("machine loads:", loads, "job lengths:", [remaining(j, 0) for j in range(len(jobs))])
# SPT   makespan 20
# LPT   makespan 14
# MWKR  makespan 13
# LWKR  makespan 21
# 13824 machine orderings, 1860 feasible, optimum 13
# machine loads: [9, 13, 8] job lengths: [7, 7, 8, 8]

Only 1,860 of the 13,824 combinations of machine sequences are feasible. The others contain a cycle. For example, "J3 before J1 on M1" together with "J1 before J3 on M2" is impossible: J3 reaches M1 only after its M2 operation, which would have to wait for J1 on M2, which in turn waits for J1 on M1. The optimum is 13, the load of the bottleneck machine M2 (index 1 in the code), so no schedule can do better. MWKR (most work remaining) reaches it by starting the jobs with the most work left first. SPT is 54% worse, because it favors operations that are quick right now, not the jobs that will keep the bottleneck busy later. In job shops, rules based on the remaining work generally beat rules based on the current operation only.

Key takeaways

• In a job shop each job has its own route. $Jm ,||, C_{\max}$ is strongly NP-hard and notoriously hard in practice.
• The disjunctive graph turns a schedule into an orientation of machine-conflict edges. Feasible means acyclic, and the makespan is the longest path.
• Giffler-Thompson with a priority rule builds active schedules quickly. Shifting bottleneck, branch and bound, constraint programming and tabu search do better.
• Look-ahead rules such as MWKR (most work remaining) usually beat myopic ones such as SPT for the makespan.

Practice

1. (Beginner) For the example, write the disjunctive pairs on machine M2 (code index 1). How many orientations do they have, and why are some of them infeasible once combined with the other machines?
2. (Intermediate) Formulate the example as a MIP with big-M disjunctive constraints and solve it with any MIP solver (for example PuLP with CBC). Check that it finds makespan 13, and try $M = 1000$ versus $M = \sum p$. Does the choice of $M$ change the solve time?
3. (Interview) "Our plant scheduler uses FIFO at every machine. How would you show management that better dispatching is worth it, before buying software?" Hint: simulate the shop (Chapter 12) with historical orders under FIFO, SPT, MWKR and a due-date rule. Compare makespan, tardiness and WIP, and check the robustness of the result to machine breakdowns.

Parallel Machine Scheduling

Why it matters. Many systems have several identical resources that can each do any job: check-in desks, CPU cores, delivery vans, operating rooms, cloud servers. The question is which job goes on which machine. Balancing the load (minimizing the makespan) is a classic NP-hard problem. It is also the setting of the first worst-case analysis of a heuristic, by Graham in 1966, which started the field of approximation algorithms.

Theory. In $Pm ,||, C_{\max}$, $n$ jobs are assigned to $m$ identical machines. The order on each machine does not matter for the makespan, so the problem is to partition the jobs into $m$ groups with the smallest maximum sum. Even $P2 ,||, C_{\max}$ is NP-hard, since it contains the partition problem. A simple lower bound holds for any schedule:

$$ C_{\max}^* \ge LB = \max\Big(\max_j p_j,\ \frac{1}{m} \sum_j p_j\Big). $$

Two greedy heuristics:

• List scheduling (LS): take jobs in any order and put each on the machine that becomes free first. Graham proved $C_{\max}^{LS} \le (2 - 1/m), C_{\max}^$. The argument: the job that finishes last started when every machine was busy, so it started no later than $\sum p / m$, and it lasts at most $\max p_j$. Both are at most $C^_{\max}$.
• LPT (longest processing time first): list scheduling with jobs sorted by decreasing $p_j$. Long jobs are placed while there is still room to balance them with short ones. Graham (1969) proved $C_{\max}^{LPT} \le \left(\frac{4}{3} - \frac{1}{3m}\right) C_{\max}^*$, and the bound is tight. For $m = 2$, jobs (3, 3, 2, 2, 2) reach it:

LPT:  M1  [ 3 ][2][2]   = 7          OPT:  M1  [ 3 ][ 3 ]   = 6
      M2  [ 3 ][2]      = 5                M2  [2][2][2]    = 6
                                     ratio 7/6 = 4/3 - 1/6

In practice LPT is usually within a few percent of the optimum. A polynomial-time approximation scheme exists (Hochbaum and Shmoys, 1987), and exact methods use branch and bound or MIP.

Other objectives and machine types:

Problem	Solution
$Pm ,\vert, pmtn ,\vert, C_{\max}$ (preemption allowed)	McNaughton's wrap-around rule: optimal value $LB$, computed in $O(n)$
$Pm ,\vert\vert, \sum C_j$	SPT, dealing jobs to machines in round-robin order: optimal
$Pm ,\vert\vert, \sum w_j C_j$	NP-hard; WSPT list scheduling is a good heuristic
$Qm ,\vert\vert, C_{\max}$ (uniform: machines with different speeds)	LPT variants assign each job to the machine that would finish it first
$Rm ,\vert\vert, \sum C_j$ (unrelated: any time matrix)	assignment problem (Section 3.6), polynomial
$Rm ,\vert\vert, C_{\max}$	NP-hard; LP rounding gives a 2-approximation

McNaughton's rule fills machine 1 up to $LB$, splits the job that crosses $LB$, and continues it at time 0 on the next machine. The two pieces of a split job never overlap, because $LB \ge p_j$.

Example: list scheduling, LPT, optimum and preemption.

import heapq
from itertools import product

def list_schedule(p, m):
    """Graham's list scheduling: next job goes to the machine that frees up first."""
    loads = [(0, i) for i in range(m)]
    for t in p:
        load, i = heapq.heappop(loads)
        heapq.heappush(loads, (load + t, i))
    return max(l for l, _ in loads)

def lpt(p, m):
    return list_schedule(sorted(p, reverse=True), m)   # longest processing time first

def optimum(p, m):                                     # brute force: m^n assignments
    best = sum(p)
    for assign in product(range(m), repeat=len(p)):
        loads = [0] * m
        for t, i in zip(p, assign):
            loads[i] += t
        best = min(best, max(loads))
    return best

for name, p, m in [("tight LPT instance", [3, 3, 2, 2, 2], 2),
                   ("ten jobs", [7, 3, 9, 5, 4, 8, 2, 6, 5, 4], 3)]:
    lb = max(max(p), sum(p) / m)                       # every schedule needs at least this
    print(f"{name}: m={m}, LB={lb:.2f}, list={list_schedule(p, m)}, "
          f"LPT={lpt(p, m)}, OPT={optimum(p, m)}")

# Preemption allowed: McNaughton's wrap-around rule reaches the lower bound exactly
p, m = [7, 3, 9, 5, 4, 8, 2, 6, 5, 4], 3
C = max(max(p), sum(p) / m)
machine, t, pieces = 0, 0.0, []
for j, rem in enumerate(p, start=1):
    while rem > 1e-9:
        piece = min(rem, C - t)                        # fill the machine up to C, wrap the rest
        pieces.append((machine, j, t, t + piece))
        rem -= piece; t += piece
        if C - t < 1e-9:
            machine, t = machine + 1, 0.0
print(f"McNaughton: makespan {C:.2f}")
for i in range(m):
    print(f"  M{i + 1}: " + ", ".join(f"J{j}[{s:.2f}-{e:.2f}]" for k, j, s, e in pieces if k == i))
# tight LPT instance: m=2, LB=6.00, list=7, LPT=7, OPT=6
# ten jobs: m=3, LB=17.67, list=20, LPT=19, OPT=18
# McNaughton: makespan 17.67
#   M1: J1[0.00-7.00], J2[7.00-10.00], J3[10.00-17.67]
#   M2: J3[0.00-1.33], J4[1.33-6.33], J5[6.33-10.33], J6[10.33-17.67]
#   M3: J6[0.00-0.67], J7[0.67-2.67], J8[2.67-8.67], J9[8.67-13.67], J10[13.67-17.67]

On the ten-job instance, list scheduling in the given order ends at 20, LPT at 19 and the optimum at 18. With integer processing times no schedule can end at 17.67, so 18 is the best possible here, and the optimum reaches it. With preemption, McNaughton's rule reaches 17.67 exactly. Job 3 runs on M2 from 0 to 1.33 and on M1 from 10 to 17.67, and job 6 is split in the same way. Preemption is realistic for CPU tasks, and much less so for a surgery.

Key takeaways

• $Pm ,||, C_{\max}$ is a load-balancing problem. It is NP-hard even for two machines, and $\max(\max_j p_j, \sum_j p_j / m)$ is a lower bound.
• List scheduling is within a factor $2 - 1/m$ of optimal. LPT is within $4/3 - 1/(3m)$ and usually much closer.
• With preemption, McNaughton's wrap-around rule is optimal and reaches the lower bound.
• SPT dealt round-robin minimizes $\sum C_j$ on identical machines, and unrelated machines with $\sum C_j$ reduce to an assignment problem.

Practice

1. (Beginner) Schedule jobs with times 8, 7, 6, 5, 4, 3, 3 on three identical machines with LPT. Compute the lower bound and say whether LPT is optimal.
2. (Intermediate) Generate 1,000 random instances of 12 jobs on 3 machines, with times uniform on 1 to 20. Compute the ratio LPT / LB and list scheduling / LB for each. Plot the two distributions and compare them with the worst-case guarantees.
3. (Interview) "A cloud service assigns incoming requests to the least-loaded of 50 servers. Why does this work well, and when does it fail?" Hint: this is online list scheduling with the $2 - 1/m$ guarantee. Think about unknown request lengths, a few very long jobs, and the "power of two choices" (sampling two random servers and picking the less loaded) used at large scale.

<a id="104-forecasting-for-operations"></a>

10.4 Forecasting for Operations

Every model in this chapter needs demand as an input: $D$ in the EOQ, the demand distribution in the newsvendor, $\sigma$ in the safety stock formula. Forecasting provides these inputs. This subchapter covers the methods that run inside most planning systems: moving averages and exponential smoothing for stable series, Holt-Winters for trend and seasonality, and regression when demand depends on known drivers such as promotions. It ends with error measures, which select among methods, detect biased forecasts and size safety stocks. The first two topics use a weekly series whose level jumps from 100 to 130 at week 20. The last three use five years of monthly demand with trend, seasonality and promotions, fitted on the first four years and tested on the fifth.

Moving Averages

Why it matters. The simplest way to forecast a stable demand is to average recent history. Averaging cancels out random noise, and using only recent history lets the forecast follow changes. Moving averages are still widely used for stable items, as a benchmark for other methods, and inside classical seasonal decomposition. They also show in its simplest form the trade-off that every forecasting method faces: stability versus responsiveness.

Theory. Assume demand fluctuates around a level that changes slowly: $y_t = \ell_t + \varepsilon_t$, with noise $\varepsilon_t$ of variance $\sigma^2$. The simple moving average (SMA) of order $N$ forecasts the next period as the mean of the last $N$ observations:

$$ F_{t+1} = \frac{1}{N} \sum_{i=0}^{N-1} y_{t-i}. $$

The choice of $N$ is a trade-off:

• Noise. If the level is constant, the forecast has variance $\sigma^2 / N$, and the forecast error $y_{t+1} - F_{t+1}$ has variance $\sigma^2 (1 + 1/N)$. Larger $N$ gives steadier forecasts.
• Lag. After a step change, the SMA needs $N$ periods to fully adjust. With a linear trend of slope $b$ per period, the average is centered $(N+1)/2$ periods in the past, so the forecast lags behind by $b(N+1)/2$ forever. Larger $N$ gives slower reactions.

A weighted moving average (WMA) gives more weight to recent data, $F_{t+1} = \sum_i w_i y_{t-i} / \sum_i w_i$ with $w_0 \ge w_1 \ge \dots$, and reacts faster for the same $N$.

The analogy is judging a basketball player's form. The last game is noisy. The whole career average misses that they have improved this season. The last ten games are a compromise.

Moving averages also have a second use: as a smoother rather than a forecaster. A centered moving average over one full season (for monthly data, a "2 x 12" average of two consecutive 12-month averages) removes seasonality and noise and leaves the trend. Dividing the data by it gives the raw seasonal indices of classical decomposition (the ratio-to-moving-average method). Never use a centered average as a forecast: it needs future data.

Example: noise versus lag. Weekly demand is 100 plus noise (standard deviation 10) until week 19, then jumps to 130. The code measures each method's error while the level is stable, how much its forecast fluctuates, and how many weeks it needs after the jump to reach 125.

import numpy as np

rng = np.random.default_rng(seed=5)
weeks = 40
level = np.where(np.arange(weeks) < 20, 100.0, 130.0)   # demand jumps from 100 to 130 at week 20
y = level + rng.normal(0, 10, weeks)                     # observed weekly demand

def moving_average(y, N, weights=None):
    """Forecast for week t = (weighted) mean of weeks t-N .. t-1."""
    w = np.ones(N) if weights is None else np.asarray(weights, float)
    f = np.full(len(y), np.nan)
    for t in range(N, len(y)):
        f[t] = np.dot(w, y[t - N:t]) / w.sum()          # weights[-1] applies to the latest week
    return f

def report(name, f):
    stable = np.abs(f[8:20] - y[8:20]).mean()            # error while the level is constant
    reach = next(t for t in range(20, weeks) if f[t] >= 125) - 20   # weeks to catch up after the jump
    print(f"{name:14} MAE stable = {stable:5.2f}  forecast sd = {np.std(f[8:20]):5.2f}"
          f"  weeks to reach 125 = {reach}")

for N in (2, 4, 8):
    report(f"SMA N={N}", moving_average(y, N))
report("WMA 1-2-3-4", moving_average(y, 4, weights=[1, 2, 3, 4]))
# SMA N=2        MAE stable = 11.21  forecast sd =  7.83  weeks to reach 125 = 2
# SMA N=4        MAE stable =  9.10  forecast sd =  3.75  weeks to reach 125 = 4
# SMA N=8        MAE stable =  8.94  forecast sd =  2.56  weeks to reach 125 = 8
# WMA 1-2-3-4    MAE stable =  9.85  forecast sd =  4.79  weeks to reach 125 = 3

Going from $N = 2$ to $N = 8$ cuts the fluctuation of the forecast by two thirds and the stable-period error by 20%, but quadruples the time needed to catch up with the new level. The error cannot fall below the noise itself (about 8 for a standard deviation of 10), so beyond some $N$ there is little left to gain and only lag to lose. The WMA with weights 1-2-3-4 sits between the SMAs with $N = 2$ and $N = 4$: faster than $N = 4$ and steadier than $N = 2$.

Key takeaways

• A moving average forecasts the next period as the mean of the last $N$ observations. Averaging reduces the noise variance by a factor $N$.
• Larger $N$ means steadier forecasts but slower reaction: a step change takes $N$ periods to absorb, and a trend is followed with a lag of $b(N+1)/2$.
• Weighted moving averages react faster for the same $N$.
• Centered moving averages are smoothers for decomposition, not forecasts.

Practice

1. (Beginner) Monthly sales were 120, 135, 128, 140, 150, 145. Compute the 3-month SMA and the 3-month WMA (weights 1, 2, 3) forecasts for month 7.
2. (Intermediate) Generate a series with a linear trend of $b = 2$ units per period plus noise. Verify numerically that the SMA forecast error (actual minus forecast) averages $b(N+1)/2$ for $N = 3, 6, 12$.
3. (Interview) "Our planners use a 12-month moving average for every product. What problems would you look for?" Hint: seasonal products (a 12-month average removes the seasonality the plan needs), growing or declining products (systematic lag), new products with short histories, and promotions inside the window.

Exponential Smoothing

Why it matters. Exponential smoothing, introduced by Brown and Holt in the 1950s for inventory control, is the most widely used forecasting method in operations. It weights all past data, with weights that decline geometrically, it needs to store only one number per item, and it has a single parameter. Planning systems that forecast hundreds of thousands of items every week rely on it. Its extensions to trend and seasonality (next topic) are among the most accurate methods in forecasting competitions.

Theory. Simple exponential smoothing (SES) updates the forecast by a fraction $\alpha \in (0, 1]$ of the last error:

$$ F_{t+1} = F_t + \alpha,(y_t - F_t) = \alpha, y_t + (1 - \alpha), F_t. $$

Unrolling the recursion shows the weights:

$$ F_{t+1} = \alpha, y_t + \alpha (1-\alpha), y_{t-1} + \alpha (1-\alpha)^2, y_{t-2} + \dots $$

Every past observation counts, with geometrically decreasing weight, which is where the name comes from. The analogy is a thermostat that corrects the temperature by a fixed fraction of the gap each time: a large $\alpha$ corrects aggressively and reacts to every fluctuation, a small $\alpha$ corrects gently and stays steady.

Properties:

• Same trade-off as the SMA. The average age of the data is $(1 - \alpha)/\alpha$, while for an SMA it is $(N - 1)/2$. Equating them gives the rule of thumb $\alpha = 2/(N + 1)$: $\alpha = 0.2$ behaves like a 9-period SMA. Typical values are 0.1 to 0.3.
• Choosing $\alpha$. Minimize the sum of squared one-step-ahead errors over the history. The optimum is reliable only if the history is representative: a past level shift drives the fitted $\alpha$ up.
• Initialization. Set $F_1$ to the first observation or to the mean of the first few. Its influence decays like $(1 - \alpha)^t$.
• Flat forecasts. The forecast is the same for every future period. SES is the optimal forecast when the data are a random walk plus noise (an ARIMA(0,1,1) process), and it lags any trend by $b(1 - \alpha)/\alpha$. For trends, use Holt's method (next topic).

Example: SES on the level-shift series. The same weekly data as the Moving Averages topic.

import numpy as np

rng = np.random.default_rng(seed=5)
weeks = 40
level = np.where(np.arange(weeks) < 20, 100.0, 130.0)   # the series of the Moving Averages topic
y = level + rng.normal(0, 10, weeks)

def ses(y, alpha, l0=None):
    """Simple exponential smoothing: f[t] is the forecast for t made at the end of t-1."""
    f = np.empty(len(y))
    l = y[0] if l0 is None else l0                       # initial level
    for t in range(len(y)):
        f[t] = l
        l = l + alpha * (y[t] - l)                       # move a fraction alpha toward the new data
    return f

for alpha in (0.1, 0.2, 0.4, 0.7):
    f = ses(y, alpha)
    stable = np.abs(f[8:20] - y[8:20]).mean()
    reach = next((t for t in range(20, weeks) if f[t] >= 125), None)
    print(f"alpha={alpha}: MAE stable = {stable:5.2f}, forecast sd = {np.std(f[8:20]):5.2f}, "
          f"weeks to reach 125 = {reach - 20 if reach else '>20'}, "
          f"equivalent SMA N = {2 / alpha - 1:4.1f}")

alphas = np.linspace(0.01, 0.99, 99)
sse = [np.sum((y[1:] - ses(y, a)[1:]) ** 2) for a in alphas]   # one-step-ahead errors
print(f"alpha minimizing in-sample SSE: {alphas[np.argmin(sse)]:.2f}")
print("weights on the last 5 observations for alpha=0.3:",
      np.round([0.3 * 0.7 ** k for k in range(5)], 3))
# alpha=0.1: MAE stable =  9.09, forecast sd =  1.15, weeks to reach 125 = >20, equivalent SMA N = 19.0
# alpha=0.2: MAE stable =  8.69, forecast sd =  2.32, weeks to reach 125 = 14, equivalent SMA N =  9.0
# alpha=0.4: MAE stable =  9.47, forecast sd =  4.58, weeks to reach 125 = 5, equivalent SMA N =  4.0
# alpha=0.7: MAE stable = 11.08, forecast sd =  7.61, weeks to reach 125 = 2, equivalent SMA N =  1.9
# alpha minimizing in-sample SSE: 0.53
# weights on the last 5 observations for alpha=0.3: [0.3   0.21  0.147 0.103 0.072]

The results mirror the moving averages: $\alpha = 0.4$ behaves like the 4-week SMA, and $\alpha = 0.7$ like the 2-week one. With $\alpha = 0.1$ the forecast is very steady (standard deviation 1.15) but still has not reached the new level 20 weeks after the jump. The fitted $\alpha = 0.53$ is high because the single level shift dominates the squared errors. On a stable stretch of data the same procedure would choose a much smaller value. The fitted parameter reflects the history it was fitted on.

Key takeaways

• SES updates the forecast by a fraction $\alpha$ of the latest error. Past data get geometrically decreasing weights.
• $\alpha$ plays the role of $1/N$: $\alpha \approx 2/(N+1)$ matches an $N$-period moving average. Typical values are 0.1 to 0.3.
• Fit $\alpha$ by minimizing one-step-ahead squared errors, and check that the history is representative.
• SES gives flat forecasts and lags trends. Use Holt or Holt-Winters for trended or seasonal series.

Practice

1. (Beginner) With $\alpha = 0.3$ and an initial forecast of 100, compute the SES forecasts after observing 110, 95, 120, 105.
2. (Intermediate) Implement Holt's linear method (the level and trend equations of the next topic) and apply it to a series with a trend of 2 units per period. Compare its long-run forecast error with SES for $\alpha = 0.2$, and check SES's lag against $b(1-\alpha)/\alpha$.
3. (Interview) "Why would a planning system with 200,000 items prefer exponential smoothing to a moving average or a machine-learning model?" Hint: storage (one number per item versus $N$), computation, robustness, the ease of explaining and overriding it, and the evidence from forecasting competitions on simple methods.

Holt-Winters Method

Why it matters. Most operational demand has a trend (growth or decline) and seasonality (weekly, monthly or yearly cycles: ice cream in summer, toys in December, call volume on Monday mornings). SES cannot follow either one. Holt-Winters extends exponential smoothing with a trend and a seasonal component. Each component is updated by its own smoothing equation. It remains one of the most accurate and widely used methods for seasonal operational series.

Theory. The method keeps three components, each updated by exponential smoothing:

• the level $\ell_t$ (deseasonalized demand now),
• the trend $b_t$ (change of the level per period),
• the seasonal indices $s_t$ (one per season position, for a season of length $m$, e.g. $m = 12$ for monthly data).

Holt's linear method (trend, no seasonality) comes first:

$$ \ell_t = \alpha, y_t + (1 - \alpha)(\ell_{t-1} + b_{t-1}), \qquad b_t = \beta,(\ell_t - \ell_{t-1}) + (1 - \beta), b_{t-1}, \qquad \hat y_{t+h} = \ell_t + h, b_t. $$

Each equation is a weighted average of new evidence and the previous estimate. The new level estimate is the observation; the previous one is last period's level advanced by the trend. The new trend estimate is the latest change in level.

Holt-Winters adds seasonality. In the multiplicative version, seasonal indices are ratios (December = 1.3 times an average month), suitable when the seasonal swings grow with the level:

$$ \begin{aligned} \ell_t &= \alpha, \frac{y_t}{s_{t-m}} + (1 - \alpha)(\ell_{t-1} + b_{t-1}) && \text{(deseasonalize, then smooth)} \ b_t &= \beta,(\ell_t - \ell_{t-1}) + (1 - \beta), b_{t-1} \ s_t &= \gamma, \frac{y_t}{\ell_t} + (1 - \gamma), s_{t-m} && \text{(update this season's index)} \ \hat y_{t+h} &= (\ell_t + h, b_t); s_{t+h-m} && (h \le m) \end{aligned} $$

The additive version replaces the divisions by subtractions and the product by a sum, and suits seasonal swings of constant size.

Choice	Use when
Additive seasonality	the seasonal swing has a constant size in units
Multiplicative seasonality	the swing grows in proportion to the level (common in sales)
Damped trend ($\hat y_{t+h} = \ell_t + (\phi + \dots + \phi^h) b_t$, $\phi < 1$)	long horizons: trends rarely continue forever, and damping is often the most accurate choice

Initialization uses the first seasons: the level is the mean of the first season, the trend is the change between the means of the first two seasons divided by $m$, and the seasonal indices are the first season's ratios to its mean. Parameters $\alpha, \beta, \gamma \in (0, 1)$ are fitted by minimizing one-step-ahead squared errors. Modern software treats these methods as ETS state-space models (error, trend, seasonal), which add likelihood-based fitting, model selection by AIC and prediction intervals.

Example: forecasting the fifth year. The monthly series has a trend, a December peak and occasional promotions (20% of months, +25% demand). Holt-Winters is fitted on four years and forecasts the fifth, compared with the seasonal naive forecast (same month last year).

import numpy as np
from scipy.optimize import minimize

# Five years of monthly demand: trend, multiplicative seasonality (December peak), promotions, noise
rng = np.random.default_rng(seed=21)
n = 60
t = np.arange(n)
season = np.array([0.80, 0.78, 0.92, 0.98, 1.02, 1.05, 0.95, 0.90, 1.00, 1.08, 1.20, 1.32])
promo = (rng.random(n) < 0.2).astype(float)            # promotion months, planned in advance
y = (400 + 5 * t) * season[t % 12] * (1 + 0.25 * promo) * rng.normal(1, 0.04, n)
train, test = slice(0, 48), slice(48, 60)              # fit on 4 years, evaluate on the 5th

def holt_winters(y, alpha, beta, gamma, m=12, horizon=0):
    """Multiplicative Holt-Winters. Returns one-step-ahead fitted values and h-step forecasts."""
    level = y[:m].mean()                               # initial level: mean of the first season
    trend = (y[m:2 * m].mean() - y[:m].mean()) / m     # initial trend: year-over-year change per month
    s = list(y[:m] / level)                            # initial seasonal indices
    fitted = np.full(len(y), np.nan)
    for t in range(m, len(y)):
        fitted[t] = (level + trend) * s[t - m]         # forecast for t made at t-1
        prev = level
        level = alpha * y[t] / s[t - m] + (1 - alpha) * (level + trend)   # deseasonalized update
        trend = beta * (level - prev) + (1 - beta) * trend
        s.append(gamma * y[t] / level + (1 - gamma) * s[t - m])
    n = len(y)
    forecast = np.array([(level + h * trend) * s[n - m + (h - 1) % m] for h in range(1, horizon + 1)])
    return fitted, forecast

def sse(params):
    fitted, _ = holt_winters(y[train], *params)
    return np.nansum((y[train] - fitted) ** 2)

res = minimize(sse, x0=[0.3, 0.1, 0.1], bounds=[(0.01, 0.99)] * 3)
alpha, beta, gamma = res.x
_, fc = holt_winters(y[train], alpha, beta, gamma, horizon=12)
snaive = y[36:48]                                      # seasonal naive: same month last year
print(f"alpha = {alpha:.2f}, beta = {beta:.2f}, gamma = {gamma:.2f}")
print("month  actual  Holt-Winters  seasonal naive")
for k in range(12):
    print(f"{k + 49:5d} {y[48 + k]:7.0f} {fc[k]:13.0f} {snaive[k]:15.0f}")
for name, f in (("Holt-Winters", fc), ("seasonal naive", snaive)):
    print(f"{name:15} test MAPE = {np.mean(np.abs(y[test] - f) / y[test]):.1%}")
# alpha = 0.01, beta = 0.01, gamma = 0.56
# month  actual  Holt-Winters  seasonal naive
#    49     529           504             438
#    50     672           575             568     <- promotion
#    51     608           602             537
#    52     621           681             606
#    53     697           744             707
#    54     872           864             782     <- promotion
#    55     631           629             576
#    56     608           692             585
#    57     715           670             593
#    58     891           691             596     <- promotion
#    59     993           845             714     <- promotion
#    60    1144           969             847     <- promotion
# Holt-Winters    test MAPE = 9.2%
# seasonal naive  test MAPE = 14.6%

Holt-Winters cuts the error of the seasonal naive forecast by more than a third, mainly because it projects the trend that the naive method ignores. The fitted parameters are informative. With $\alpha = \beta = 0.01$, the level and trend barely move after initialization, because the underlying growth is steady. The high $\gamma = 0.56$ makes the seasonal indices chase last year's values, including last year's promotions. May and June of year 4 had promotions, so the forecast for May of year 5 is too high (744 against 697). The last three months have promotions the method cannot know about, and it under-forecasts them by 15 to 22%. Promotions are planned in advance, so this information is available. The next topic uses it.

Key takeaways

• Holt's method adds a smoothed trend to SES, and Holt-Winters adds smoothed seasonal indices, each with its own parameter.
• Use multiplicative seasonality when swings grow with the level, additive when they are constant. Damped trends are safer for long horizons.
• Parameters are fitted by minimizing one-step-ahead errors. Their values say how stable each component is.
• Holt-Winters extrapolates the series' own past. It cannot anticipate known future events such as promotions, and it may absorb past events into the seasonal indices.

Practice

1. (Beginner) A multiplicative Holt-Winters model currently has level 500, trend 4 per month and a December index of 1.30. Today is the end of September. What is the forecast for December?
2. (Intermediate) Implement the additive version and fit it to the same training data. Compare the test MAPE with the multiplicative version and explain the difference using the way the data were generated.
3. (Interview) "A Holt-Winters model forecast strong growth for a product that then flattened, and the warehouse is now full. What went wrong, and how would you guard against it?" Hint: undamped trends extrapolate indefinitely. Consider damped trends, limits on the trend parameter, monitoring with a tracking signal (last topic), and judgmental review of long-horizon forecasts.

Regression-Based Forecasting

Why it matters. Time-series methods extrapolate a series from its own past. But demand often depends on things that are known or planned in advance: price, promotions, holidays, advertising, the weather forecast, the number of working days in the month. Regression links demand to these drivers. It forecasts the effect of a planned promotion, and it answers "what if" questions ("what happens if we cut the price 10%?") that pure time-series methods cannot.

Theory. A linear regression model writes demand as a combination of explanatory variables plus noise:

$$ y_t = \beta_0 + \beta_1 x_{1t} + \dots + \beta_k x_{kt} + \varepsilon_t, $$

and estimates the $\beta$ coefficients by least squares, minimizing $\sum_t \big(y_t - \hat y_t\big)^2$. This is an unconstrained convex quadratic problem (Section 5.1) with the closed-form solution $\hat\beta = (X^\top X)^{-1} X^\top y$. Numerically, a QR or SVD-based solver such as numpy.linalg.lstsq is used instead of the inverse. Three types of variables appear in operational models:

• Trend: $x_t = t$, possibly with a slope change at known dates.
• Seasonality: $m - 1$ dummy variables (for monthly data: 11 indicators, one per month except a reference month). Using all 12 together with an intercept makes $X$ singular, the "dummy variable trap".
• Causal drivers: price, promotion flags, holidays, temperature.

Taking logs, $\log y_t = \dots$, turns the effects into percentages: a coefficient $\delta$ on a promotion flag means demand is multiplied by $e^{\delta}$. This matches multiplicative seasonality.

An approximate 95% prediction interval is $\hat y \pm 1.96, s$, where $s$ is the residual standard error. This ignores the uncertainty in the estimated coefficients, which is small when there are many more observations than coefficients.

Pitfalls:

Pitfall	Symptom	Remedy
Future drivers unknown	the model needs next month's temperature or competitor price	forecast the driver, use scenarios, or use only planned drivers
Autocorrelated residuals	residuals of consecutive periods are correlated; intervals too narrow	add lagged terms, model the errors (ARIMA errors), check the trend form
Extrapolation	drivers outside the historical range (a 50% discount never tried before)	treat as a guess; run experiments
Overfitting	many variables, short history; excellent fit, poor forecasts	hold-out testing, fewer variables, regularization
Correlation is not causation	price was cut because demand was low	design experiments or use causal-inference methods

Example: trend, seasonal dummies and a promotion flag. The model is fitted to the log of the first four years of the Holt-Winters series and forecasts the fifth year using the planned promotion calendar.

import numpy as np
# Reuses y, promo, n, train and test from the Holt-Winters topic.

t = np.arange(n)
X = np.column_stack([np.ones(n), t] +                       # intercept and linear trend
                    [(t % 12 == k).astype(float) for k in range(1, 12)] +   # 11 month dummies
                    [promo])                                # causal driver: promotion flag
logy = np.log(y)                                            # log: seasonality and promo act as % effects

coef, *_ = np.linalg.lstsq(X[train], logy[train], rcond=None)
resid = logy[train] - X[train] @ coef
s = np.sqrt(resid @ resid / (48 - X.shape[1]))              # residual standard error (log scale)
print(f"trend = {np.expm1(coef[1]):+.2%} per month, promotion uplift = {np.expm1(coef[-1]):+.1%}")
print(f"December vs January: {np.exp(coef[12]):.2f}x,  residual sd = {s:.3f}")
print(f"lag-1 autocorrelation of residuals = {np.corrcoef(resid[:-1], resid[1:])[0, 1]:+.2f}")

pred = np.exp(X[test] @ coef)
lo, hi = pred * np.exp(-1.96 * s), pred * np.exp(1.96 * s)  # approximate 95% interval
print("month  promo  actual  forecast   95% interval")
for k in range(12):
    print(f"{49 + k:5d} {int(promo[48 + k]):6d} {y[48 + k]:7.0f} {pred[k]:9.0f}   [{lo[k]:4.0f}, {hi[k]:4.0f}]")
print(f"test MAPE = {np.mean(np.abs(y[test] - pred) / y[test]):.1%}, "
      f"coverage = {np.mean((y[test] >= lo) & (y[test] <= hi)):.0%}")

coef0, *_ = np.linalg.lstsq(X[train][:, :-1], logy[train], rcond=None)   # same model without promo
pred0 = np.exp(X[test][:, :-1] @ coef0)
print(f"without the promotion variable: test MAPE = {np.mean(np.abs(y[test] - pred0) / y[test]):.1%}")
# trend = +0.94% per month, promotion uplift = +22.2%
# December vs January: 1.69x,  residual sd = 0.040
# lag-1 autocorrelation of residuals = +0.34
# month  promo  actual  forecast   95% interval
#    49      0     529       509   [ 471,  551]
#    50      1     672       620   [ 573,  671]
#    51      0     608       600   [ 555,  650]
#    52      0     621       663   [ 612,  717]
#    53      0     697       684   [ 632,  740]
#    54      1     872       896   [ 828,  969]
#    55      0     631       634   [ 586,  685]
#    56      0     608       639   [ 591,  691]
#    57      0     715       686   [ 634,  743]
#    58      1     891       878   [ 812,  950]
#    59      1     993      1047   [ 967, 1132]
#    60      1    1144      1165   [1077, 1261]
# test MAPE = 3.5%, coverage = 92%
# without the promotion variable: test MAPE = 8.7%

The regression recovers the structure of the data: a promotion uplift of 22% (the true value is 25%) and a December-to-January ratio of 1.69 (true: 1.32 / 0.80 = 1.65). With the promotion calendar, the test error falls to 3.5%, compared with 9.2% for Holt-Winters. Without the promotion variable, the same regression is no better than Holt-Winters (8.7%). The gain comes from the information, not from the method. The 95% intervals cover 11 of the 12 months (month 50 is just outside). The residual autocorrelation of +0.34 is a warning sign: the data grow linearly, $400 + 5t$, but a linear trend in logs is exponential growth. The misspecified trend leaves a pattern in the residuals, which a plot of residuals over time would show.

Key takeaways

• Regression forecasts demand from drivers: trend, seasonal dummies ($m - 1$ of them) and causal variables such as price and promotions.
• Logs turn coefficients into percentage effects and match multiplicative seasonality.
• Regression wins when drivers are known in advance and matter. The future values of the drivers must be known or forecast.
• Check residuals for autocorrelation and test on held-out data. A good in-sample fit proves little.

Practice

1. (Beginner) A fitted model is $\log y = 6.1 + 0.01,t + 0.18,\text{promo} + s_{\text{month}}$, with the December coefficient $s_{\text{Dec}} = 0.35$. By what percentage does a promotion raise demand? How much higher is December than the reference month?
2. (Intermediate) The residual autocorrelation suggests that the trend has the wrong form. Fit two variants: the model in logs with an extra quadratic trend term $t^2$, and the model in levels ($y$ instead of $\log y$, so all effects are additive). Which variant removes the residual autocorrelation, and which has the lower test MAPE?
3. (Interview) "Marketing claims their promotions lift sales by 40%, based on comparing promotion weeks with other weeks. Why might this be wrong?" Hint: promotions are scheduled in high-season weeks (confounding with seasonality), sales are pulled forward from the following weeks (post-promotion dip), and other products are cannibalized. Think about what a regression with seasonal terms and lagged promotion effects would show.

Forecast Error Measures

Why it matters. A forecast is useful only if its error is known. Error measures choose between methods, tune their parameters, detect forecasts that have drifted into bias, and size safety stock: the standard deviation of forecast errors is the $\sigma$ of the reorder point formula (Section 10.1). Choosing the wrong measure can make a bad method look good, and some popular measures break down on common data, such as intermittent demand.

Theory. With error $e_t = y_t - F_t$ (actual minus forecast, positive when the forecast is too low):

Measure	Formula	Measures	Notes
ME (mean error, bias)	$\frac{1}{n}\sum e_t$	systematic over- or under-forecasting	should be near 0; errors of opposite sign cancel
MAE / MAD	$\frac{1}{n}\sum \lvert e_t\rvert$	typical error size, in units	robust; optimized by the median
RMSE	$\sqrt{\frac{1}{n}\sum e_t^2}$	error size, penalizing large errors	optimized by the mean; feeds safety stock
MAPE	$\frac{100}{n}\sum \lvert e_t\rvert / y_t$	error in percent	scale-free; undefined when $y_t = 0$; puts a heavier penalty on over-forecasts
MASE	$\text{MAE} / \text{MAE}_{\text{naive, in-sample}}$	error relative to a naive benchmark	scale-free, works with zeros; below 1 beats the benchmark

Some rules for using them:

• Use out-of-sample errors. Fit on one period, measure on a later one (a hold-out set or rolling origin). In-sample errors reward overfitting.
• Match the measure to the decision. RMSE is minimized by the mean forecast and MAE by the median. If the costs of being wrong are asymmetric, as in the newsvendor, use a quantile (pinball) loss.
• Compare across items with scale-free measures (MAPE when demand is never close to zero, otherwise MASE).
• For normal errors, $\sigma \approx 1.25 \times \text{MAD}$ ($\sqrt{\pi/2} \approx 1.2533$). Older systems estimated $\sigma$ for safety stock this way because the MAD is easy to update with exponential smoothing.

Monitoring bias. The tracking signal is the running sum of errors divided by the running MAD:

$$ TS_t = \frac{\sum_{i \le t} e_i}{\text{MAD}_t}. $$

An unbiased forecast has errors of both signs, and the signal wanders around zero. A biased one has errors that pile up on one side, and the signal drifts steadily. A common rule flags a forecast when $|TS|$ exceeds about 4 (limits between 3 and 6 are used). The first few values are mechanical: if every error so far has the same sign, $|TS_t| = t$, so the signal needs a few periods before it means anything.

Example: comparing five methods on the test year. The forecasts of the two previous topics, plus two naive benchmarks: the last observed value (December of year 4) for every month, and the same month last year.

import numpy as np
# Reuses y, train, test and the Holt-Winters forecasts fc from the Holt-Winters topic,
# and pred (with promotions) and pred0 (without) from the Regression-Based Forecasting topic.

actual = y[test]
scale = np.mean(np.abs(y[12:48] - y[0:36]))   # in-sample MAE of the seasonal naive method

def measures(f):
    e = actual - f                            # error = actual - forecast (positive: under-forecast)
    return {"ME": e.mean(), "MAE": np.abs(e).mean(), "RMSE": np.sqrt((e ** 2).mean()),
            "MAPE": np.mean(np.abs(e) / actual) * 100, "MASE": np.abs(e).mean() / scale}

methods = {"naive (last value)": np.full(12, y[47]), "seasonal naive": y[36:48],
           "Holt-Winters": fc, "regression": pred0, "regression + promo": pred}
print(f"{'method':>19} {'ME':>7} {'MAE':>6} {'RMSE':>6} {'MAPE%':>6} {'MASE':>5}")
for name, f in methods.items():
    m = measures(f)
    print(f"{name:>19} " + " ".join(f"{m[k]:{w}.{d}f}" for k, w, d in
          (("ME", 7, 1), ("MAE", 6, 1), ("RMSE", 6, 1), ("MAPE", 6, 1), ("MASE", 5, 2))))

# Tracking signal = running sum of errors / running MAD; |TS| > 4 flags a biased forecast
for name in ("seasonal naive", "regression + promo"):
    e = actual - methods[name]
    ts = np.cumsum(e) / (np.cumsum(np.abs(e)) / np.arange(1, 13))
    print(f"tracking signal, {name:>18}: " + " ".join(f"{v:+.1f}" for v in ts))

# MAPE breaks down on intermittent demand
spare = np.array([0, 2, 0, 1, 0, 0, 3, 0, 1, 0])
flat = np.full(10, 0.7)
with np.errstate(divide="ignore"):
    mape = np.mean(np.abs(spare - flat) / spare)       # division by zero in the zero-demand weeks
print(f"spare parts: MAE = {np.abs(spare - flat).mean():.2f}, MAPE = {mape}")
#              method      ME    MAE   RMSE  MAPE%  MASE
#  naive (last value)   -98.6  183.9  203.2   27.0  2.50
#      seasonal naive   119.3  121.0  159.2   14.6  1.64
#        Holt-Winters    42.8   74.6   98.8    9.2  1.01
#          regression    35.2   68.1   88.9    8.7  0.92
#  regression + promo    -3.3   25.6   30.1    3.5  0.35
# tracking signal,     seasonal naive: +1.0 +2.0 +3.0 +4.0 +4.7 +5.7 +6.7 +7.7 +8.7 +9.8 +10.8 +11.8
# tracking signal, regression + promo: +1.0 +2.0 +3.0 +1.3 +1.9 +1.1 +1.1 -0.2 +0.9 +1.6 -0.7 -1.5
# spare parts: MAE = 0.84, MAPE = inf

All measures rank the methods the same way here, but the bias column explains why they differ. Repeating December's peak all year over-forecasts by 99 units a month on average (ME = -98.6). The seasonal naive method ignores growth and under-forecasts by 119. Holt-Winters and the regression without promotions still under-forecast by 35 to 43 units, mostly in the promotion months. Only the regression with promotions is essentially unbiased. The tracking signal of the seasonal naive forecast reaches 4 in the fourth month and keeps climbing, a clear alarm. The regression's signal reaches 3 by mechanical start-up, then falls back and stays within ±2. A MASE of 0.35 means the regression's error is about a third of the benchmark's. On the spare-parts series the MAPE is infinite because of the zero-demand weeks, while the MAE remains meaningful.

Key takeaways

• Measure errors out of sample. Use ME for bias, MAE or RMSE for size in units, and MAPE or MASE to compare across items.
• RMSE favors the mean and MAE the median. Pick the measure that matches the cost of errors.
• MAPE fails on zero or near-zero demand. MASE is the robust scale-free alternative.
• The tracking signal detects bias early. The standard deviation of forecast errors, not of demand, should drive safety stock.

Practice

1. (Beginner) Actual demand was 100, 120, 90, 110 and the forecasts were 105, 110, 100, 100. Compute ME, MAE, RMSE, MAPE and the tracking signal after each period.
2. (Intermediate) Using the RMSE of the regression with promotions, compute the safety stock for a 95% cycle service level with a one-month lead time (Section 10.1). Repeat with the RMSE of the seasonal naive forecast, and convert the difference into annual holding cost at 25% of a unit value of 20.
3. (Interview) "Our forecast accuracy KPI is 1 - MAPE, and the team's bonus depends on it. What behavior might this encourage?" Hint: MAPE penalizes over-forecasts more than under-forecasts of the same size, explodes on low-volume items, and can be improved by excluding hard items. Discuss a bias measure, MASE and weighting by volume or value.

<a id="11-heuristics-and-metaheuristics"></a>