Solving Order Quantity and Safety Stock Together
Every warehouse manager eventually runs into the same tension: order in large batches and you pay too much to hold stock; order in tiny batches and you pay too much in ordering fees; keep too little buffer stock and you risk running out during the lead time. This article works through a concrete numerical example that solves the order quantity and safety stock simultaneously, rather than treating them as two separate problems, and visualizes the full cost landscape in 3D.
The Business Scenario
A distribution center sells a mid-volume SKU with the following characteristics:
- Annual demand: 12,000 units/year
- Ordering cost: $45 per purchase order
- Holding cost: $4 per unit per year
- Shortage (backorder) cost: $15 per unit short
- Average daily demand: ~32.9 units, with a daily demand standard deviation of 5 units
- Supplier lead time: 14 days
The question: how many units should we order each time, and how much safety stock should we carry, to minimize total annual cost?
Mathematical Formulation
Classical Economic Order Quantity (EOQ)
Ignoring uncertainty for a moment, the trade-off between ordering cost and holding cost gives the well-known EOQ formula:
$$
Q^{*} = \sqrt{\frac{2DS}{H}}
$$
where $D$ is annual demand, $S$ is the fixed cost per order, and $H$ is the holding cost per unit per year.
Adding Demand Uncertainty: Safety Stock
Because daily demand is random, the demand realized during the lead time $L$ is also random. If lead-time demand has standard deviation $\sigma_L$, we define a safety factor $z$ (a standard normal quantile) so that:
$$
SS = z,\sigma_L, \qquad R = \bar{d}L + z,\sigma_L
$$
where $SS$ is the safety stock, $R$ is the reorder point, and $\bar{d}$ is average daily demand. A larger $z$ means a higher target service level, but also more holding cost.
Total Annual Cost
Combining ordering cost, cycle-stock holding cost, safety-stock holding cost, and the expected cost of running out during the lead time (the Hadley–Whitin formulation), the total annual cost as a function of both decision variables $Q$ and $z$ is:
where $C_s$ is the shortage cost per unit, and $L(z)$ is the standard normal loss function:
$$
L(z) = \phi(z) - z\big(1-\Phi(z)\big)
$$
with $\phi$ the standard normal PDF and $\Phi$ the standard normal CDF. $L(z)$ represents the expected number of units short (in standard-deviation units) per replenishment cycle.
The goal is to jointly minimize $TC(Q,z)$ over both $Q$ and $z$ — this is what makes the problem genuinely two-dimensional rather than two separate one-dimensional problems.
Source Code
1 | # ============================================================ |
Console Output
================================================== INVENTORY OPTIMIZATION RESULTS ================================================== Classical EOQ (no uncertainty) : 519.62 units Optimal order quantity Q* : 526.11 units Optimal safety factor z* : 2.2671 Implied service level : 98.83 % Optimal safety stock SS* : 42.41 units Optimal reorder point R* : 502.69 units Minimum total annual cost : $ 2274.07 ==================================================
Code Walkthrough
Section 1 – Parameters. All business inputs are declared as plain floats at the top: demand, ordering cost, holding cost, shortage cost, and the demand-variability figures (daily demand std dev and lead time). Keeping these separate from the model logic makes the script easy to re-run with a different SKU’s numbers.
Section 2 – Cost model. loss_function(z) implements the standard normal unit loss function $L(z)$, which converts a safety factor into an expected number of units short per cycle. total_cost(x) is the scalar objective function that scipy.optimize.minimize calls; it unpacks x = [Q, z] and sums the four cost terms described in the math section above. A guard clause returns a very large penalty if Q is non-positive, which keeps the optimizer away from degenerate values. total_cost_vec(Q, z) is a NumPy-vectorized twin of the same function — it accepts full arrays for Q and z and returns an array of costs with no Python-level loop, which is what makes the 3D surface below build instantly instead of looping over thousands of grid points one at a time.
Section 3 – Joint optimization. Rather than solving order quantity and safety stock as two independent problems, minimize searches over both Q and z at once, using the classical EOQ value and a 95%-service-level guess (z = 1.65) as the starting point. L-BFGS-B is used because it supports box constraints (bounds), which keeps Q and z in economically sensible ranges during the search. The results are then translated into safety stock and reorder point using the formulas from the math section.
Section 4 – Cost surface. A 120×120 grid of (Q, z) combinations is built with np.meshgrid, and the entire cost surface is evaluated in one vectorized call to total_cost_vec. This avoids a nested double loop (120×120 = 14,400 evaluations) and keeps the whole computation well under a second.
Section 5 – 3D plot. plot_surface draws the cost landscape, and the optimum found in Section 3 is marked as a single cyan point. The surface makes the trade-off visually obvious: moving along the $Q$-axis away from the optimum increases cost because of the ordering/holding trade-off, while moving along the $z$-axis away from the optimum increases cost because of the safety-stock/shortage trade-off. The dark theme (background, panes, tick colors) matches a typical technical-blog dark layout.
Section 6 – 2D breakdown. Fixing $z$ at its optimal value, this chart decomposes total cost into its four components as a function of $Q$ alone — this is the classic “U-shaped EOQ curve” but now shown alongside the safety-stock and shortage-cost lines so the reader can see how much of the total cost each component actually contributes at the optimum.
The 3D surface is bowl-shaped, with a single global minimum in the interior of the plotted region — that’s the optimum reported by the optimizer. Moving toward small $Q$ makes the surface rise steeply because ordering cost explodes as $D/Q$; moving toward small $z$ makes it rise because expected shortage cost grows. The optimum sits at a moderate order quantity with a fairly high safety factor, since in this example the shortage cost is significant relative to holding cost.
In the 2D breakdown, the white “Total Cost” curve is clearly the sum of the other four curves, and its minimum lines up exactly with the cyan marker. Ordering cost falls steadily as $Q$ grows, cycle holding cost rises linearly, and the safety-stock holding cost stays flat (since $z$ is fixed in this chart) while shortage cost falls as larger, more frequent… actually less frequent orders reduce the number of cycles per year in which a shortage can occur.
Interpretation of the Numbers
For this example, the model finds an optimal order quantity of roughly 526 units per order — close to, but not identical to, the classical EOQ of about 520 units, because uncertainty slightly shifts the ideal batch size. The optimal safety factor comes out to about $z \approx 2.27$, which corresponds to an implied service level of roughly 98.8%. That translates into a safety stock of about 42 units and a reorder point of about 503 units, at a minimum total annual cost of roughly $2,274.
The relatively high service level here is a direct consequence of the numbers chosen: a shortage cost of $15 per unit is considerably higher than the $4 holding cost per unit, so the optimizer prefers to carry a bit more safety stock rather than risk running short. Changing that ratio — for example, lowering the shortage cost or raising the holding cost — would pull the optimal $z$ down and reduce the safety stock accordingly. That sensitivity is exactly what the 3D surface makes visible at a glance: the cost landscape’s shape along the $z$-axis directly reflects how expensive shortages are relative to holding inventory.
