Placing stock across a distribution network with a mixed-integer quadratic model in PySCIPOpt

A food distributor's inventory decisions ripple into centre capacity, transport cost, order efficiency, stock-outs and overstock at once. Conservative policies keep shelves full but fill storage, multiply safety stock and move goods between centres that did not need to hold them. Simple rules each see one part of the problem: a reorder point per centre ignores that demand pooled at one hub needs less buffer than the same demand split across several, and a cheapest-lane assignment ignores pallets and capacity.

In this post we walk through the model built for the study. Recency-weighted demand statistics per SKU and centre feed a per-SKU feasibility pre-check, and then one mixed-integer quadratically constrained program (MIQCP) decides storage, single sourcing and whole-pallet stock together, with a pooled safety-stock constraint that makes the trade between consolidation and transport explicit. The model is written in PySCIPOpt and solved with SCIP, which handles the nonconvex quadratic constraint directly; a second-order-cone reformulation is the next step for scale. This was a feasibility study run in Jupyter notebooks, not a deployed service.

Solution overview

The study is a notebook pipeline rather than a service. Pre-built demand history per SKU and centre, lane cost and permission matrices and the item master are loaded with pandas; demand statistics are computed; SKUs that no allowed set of centres can serve are set aside; and the remaining SKUs go into one model, whose solution is unpacked into tables for the business.

The numbered steps in the diagram:

1. Daily demand per SKU × centre over a few months, lane costs, allowed and blocked lanes, lead times, minimum order quantities, case packs, storage conditions and pallet capacities are loaded with pandas.
2. NumPy computes a recency-weighted mean and standard deviation of daily demand per SKU and centre, counted from the first day with nonzero demand.
3. For each SKU, itertools combinations of one, two and three centres are checked for lane coverage and capacity; SKUs that no combination can serve are set aside.
4. Binary storage and sourcing variables, integer pallets and continuous safety-stock and stock-level variables are built in PySCIPOpt, with a quadratic pooled-safety-stock constraint.
5. SCIP solves the nonconvex MIQCP under a time limit and an optimality-gap target.
6. The solution is unpacked into pandas tables — placement, sourcing, pallets and capacity use per centre — and summarized in the notebook.

Technology stack

Step 1: Turn demand history into recency-weighted statistics

The model needs, for each SKU at each centre, a daily demand mean μ and standard deviation σ: lead-time demand comes from μ, safety stock from σ. Both are computed from a few months of daily history, with recent months weighted more heavily so the statistics follow current behaviour rather than averaging it away.

History is counted from the first day with nonzero demand. An item that started selling at a centre partway through the window would otherwise carry weeks of structural zeros, which would understate its mean and overstate its variability.

d_t     daily demand of one SKU at one centre, t from its first nonzero day to today
ω_t     recency weight of day t (by month, recent months heavier),   Σ_t ω_t = 1

μ  = Σ_t ω_t · d_t                   σ² = Σ_t ω_t · (d_t − μ)²

import numpy as np
import pandas as pd

def weighted_stats(daily, month_weights, days_per_month=30):
    """Recency-weighted mean and sd of daily demand (oldest day first),
    counted from the first day with nonzero demand.
    month_weights[0] applies to the most recent month."""
    nz = np.flatnonzero(daily)
    if nz.size == 0:
        return 0.0, 0.0
    d = np.asarray(daily[nz[0]:], dtype=float)
    month = np.arange(len(d))[::-1] // days_per_month     # 0 = most recent
    w = np.asarray(month_weights)[np.minimum(month, len(month_weights) - 1)]
    mu = np.average(d, weights=w)
    sd = np.sqrt(np.average((d - mu) ** 2, weights=w))
    return mu, sd

def demand_stats(history: pd.DataFrame, month_weights):
    """history: one row per (sku, centre, day) with a qty column, days complete."""
    rows = []
    for (sku, centre), g in history.sort_values("day").groupby(["sku", "centre"]):
        mu, sd = weighted_stats(g["qty"].to_numpy(), month_weights)
        rows.append((sku, centre, mu, sd))
    return pd.DataFrame(rows, columns=["sku", "centre", "mu", "sd"])

Simplified. The month weights are a parameter of the study and are not shown.

Step 2: Check each SKU's feasibility before the solve

A single SKU that no allowed centre can serve makes the whole model infeasible — and an infeasible model over thousands of SKUs and tens of centres does not say which SKU is at fault. So feasibility is checked per SKU before the model is built. The check first tries one storage centre: can it reach every demand point over allowed lanes and hold the required stock within capacity? If not, it tries pairs, then triples of centres.

from itertools import combinations

def can_serve(sku, combo, points, allowed, need, cap):
    """Every demand point reachable from some centre in combo over an
    allowed lane, and the stock each centre would hold fits its capacity."""
    covered = all(any(allowed[sku, c, j] for c in combo) for j in points)
    return covered and all(need(sku, c, combo) <= cap[c] for c in combo)

def precheck(skus, centres, points, allowed, need, cap, max_k=3):
    """Try one storage centre, then pairs, then triples."""
    feasible, excluded = [], []
    for sku in skus:
        ok = any(can_serve(sku, combo, points[sku], allowed, need, cap)
                 for k in range(1, max_k + 1)
                 for combo in combinations(centres, k))
        (feasible if ok else excluded).append(sku)
    return feasible, excluded

Simplified. need() estimates the pallets a centre would hold if the combination served the SKU.

SKUs that fail are set aside and listed separately; everything else goes into the model.

Why not let the solver find out? An infeasibility certificate for a model this size points at constraints, not at a SKU. A per-SKU check is cheap, explains itself, and keeps the full model solvable.

Step 3: Formulate storage, sourcing and pallets as one MIQCP

For every SKU the model decides three things together: storage — which centres hold it, at most K of them; sourcing — which stocking centre serves each demand point, exactly one, over allowed lanes only; and size — how many whole pallets each stocking centre keeps.

i ∈ SKUs,   c ∈ centres that can store,   j ∈ centres with demand (demand points)

y_icj ∈ {0,1}   c supplies j with item i           z_ic ∈ {0,1}   item i stocked at c
p_ic  ∈ ℤ₊      pallets of i at c                  S_ic, R_ic ≥ 0  safety stock, stock level

min   Σ cost_cj · w_i · μ_ij · H · y_icj   +   w_s · Σ p_ic

s.t.  Σ_c y_icj = 1                          one source per demand point (allowed lanes only)
      y_icj ≤ z_ic,     z_ic ≤ y_icc         only stocking centres supply; they serve themselves
      Σ_c z_ic ≤ K                           at most K stocking centres per SKU
      S_ic² = LT_i · Σ_j σ_ij² · y_icj       pooled safety stock
      R_ic ≥ LT_i · Σ_j μ_ij · y_icj + 3 · S_ic
      R_ic ≥ MOQ_i · b_i · z_ic              never below the minimum order
      b_i · p_ic ≥ R_ic,   b_i · p_ic ≤ R_ic + b_i − 1        p_ic = ⌈R_ic / b_i⌉
      p_ic ≤ M · z_ic,     Σ_i p_ic ≤ cap_c  pallets only where stocked; centre capacity

The objective is transport cost over a planning horizon H plus a holding term per pallet. The storage weight w_s is the lever between the two: the higher it is, the harder the model consolidates.

Why pooling drives the trade

Safety stock at a hub covers the combined uncertainty of the demand points it serves. For independent demand, variances add, so the buffer grows with the square root of the summed variances rather than with their sum. Two points with equal σ need 3·σ·√LT each if stocked separately — 6·σ·√LT in total — but 3·σ·√(2·LT) ≈ 4.24·σ·√LT if one hub serves both, about 29% less. Consolidating a SKU into fewer hubs therefore saves pallets even when it adds transport, and the model weighs exactly that trade, SKU by SKU, under shared capacity.

Whole pallets matter as much as the square root. Lead-time demand plus three standard deviations gives a stock level R; for integer quantities, the pair b·p ≥ R and b·p ≤ R + b − 1 forces p to be exactly ⌈R/b⌉. The minimum-order floor keeps small SKUs realistic, p ≤ M·z ties pallets to the storage decision, and a stocking centre always serves its own demand.

Step 4: Build and solve the model in PySCIPOpt

The model is written directly in PySCIPOpt. Sourcing variables are created only for allowed lanes, so blocked lanes never enter the model, and the quadratic safety-stock equality is added like any other constraint.

from pyscipopt import Model, quicksum
def build(I, C, J, d, K, w_s, M):
    m = Model("placement")
    IC = [(i, c) for i in I for c in C]
    y = {(i, c, j): m.addVar(vtype="B") for i, c in IC for j in J[i] if d.allowed[i, c, j]}
    z = {k: m.addVar(vtype="B") for k in IC}
    p = {k: m.addVar(vtype="I", lb=0) for k in IC}
    S = {k: m.addVar(lb=0) for k in IC}
    R = {k: m.addVar(lb=0) for k in IC}
    for i in I:
        lt, b, mu, var = d.lead[i], d.per_pallet[i], d.mu[i], d.var[i]
        for j in J[i]:
            m.addCons(quicksum(y[i, c, j] for c in C if (i, c, j) in y) == 1)
        m.addCons(quicksum(z[i, c] for c in C) <= K)
        for c in C:
            Y = [(j, y[i, c, j]) for j in J[i] if (i, c, j) in y]
            for _, v in Y:
                m.addCons(v <= z[i, c])
            m.addCons(S[i, c] * S[i, c] == lt * quicksum(var[j] * v for j, v in Y))
            m.addCons(R[i, c] >= lt * quicksum(mu[j] * v for j, v in Y) + 3 * S[i, c])
            m.addCons(R[i, c] >= d.moq[i] * b * z[i, c])
            m.addCons(b * p[i, c] >= R[i, c])
            m.addCons(b * p[i, c] <= R[i, c] + b - 1)
            m.addCons(p[i, c] <= M * z[i, c])
    for c in C:
        m.addCons(quicksum(p[i, c] for i in I) <= d.cap[c])
    transport = quicksum(d.cost[c, j] * d.w[i] * d.mu[i][j] * d.H * v
                         for (i, c, j), v in y.items())
    m.setObjective(transport + w_s * quicksum(p.values()), "minimize")
    return m, y, z, p

Simplified. The self-supply constraint and storage-condition filters are omitted; d holds the prepared data.

SCIP solves the model under a time limit and a relative-gap target. Because the safety-stock constraint is a nonconvex quadratic equality, SCIP branches on it as well as on the integer variables; the gap target lets a large instance stop at a provably good solution instead of running to proven optimality.

import pandas as pd

m, y, z, p = build(I, C, J, data, K, w_s, M)
m.setParam("limits/time", time_limit)          # seconds
m.setParam("limits/gap", gap_target)           # relative optimality gap
m.optimize()

if m.getNSols() == 0:
    raise RuntimeError(f"no solution: {m.getStatus()}")
sol = m.getBestSol()
val = lambda v: m.getSolVal(sol, v)
placement = pd.DataFrame(
    [(i, c, round(val(p[i, c]))) for (i, c) in p if val(z[i, c]) > 0.5],
    columns=["sku", "centre", "pallets"])
sourcing = pd.DataFrame(
    [(i, c, j) for (i, c, j), v in y.items() if val(v) > 0.5],
    columns=["sku", "centre", "demand_point"])
usage = placement.groupby("centre")["pallets"].sum()
print(m.getStatus(), m.getObjVal(), m.getGap())

Simplified. The limits are study settings and are not shown.

Step 5: Convexify pooled safety stock as a second-order cone

The nonconvexity comes from one place: the equality S² = LT·Σσ²y. Two observations remove it. Because y is binary, y = y², so the right-hand side is the squared norm of the vector √LT·σ∘y. And outside its own definition, S appears only as a lower bound on the stock level, R ≥ … + 3·S, so relaxing the equality to an inequality does not change the optimal cost: any solution of the relaxed model stays feasible when S is lowered to its exact value.

as solved (nonconvex)    S_ic² = LT_i · Σ_j σ_ij² · y_icj
y binary  ⇒  y = y²      LT_i · Σ_j σ_ij² · y_icj = ‖ √LT_i · σ_i ∘ y_ic ‖₂²
next step (convex)       S_ic ≥ ‖ √LT_i · σ_i ∘ y_ic ‖₂            second-order cone

The result is a mixed-integer second-order-cone program whose continuous relaxation is convex, which is what branch-and-bound needs to scale. The reformulation was identified in the study as the next step; the study itself solved the direct form.

# As solved in the study: a nonconvex quadratic equality
m.addCons(S[i, c] * S[i, c] == lt * quicksum(var[j] * v for j, v in Y))

# Next step: a second-order cone. With y binary, var[j] * v equals
# (sqrt(var[j]) * v) ** 2, and S only needs to be large enough.
m.addCons(quicksum(lt * var[j] * v * v for j, v in Y) <= S[i, c] * S[i, c])

Simplified. With S ≥ 0 the second constraint describes a convex second-order cone that a conic-aware solver can exploit.

Step 6: Wrap the model in a decision environment

The study did not assume the optimization model was the answer. Reinforcement learning, MILP and search-based optimization were reviewed against each other; data marts were built to join inventory, demand and capacity information; and simulation structures were designed to evaluate how a change in inventory policy moves capacity, cost and service level.

Together they form the environment in which inventory policy alternatives can be tested before they are tried on the floor.

Results

As a feasibility study, the outcome was a case rather than a deployment. The optimization model confirmed the potential to secure logistics-centre capacity and reduce transportation cost, and the project proposed a shift from conservative inventory management to demand-based inventory operation. No operational savings are claimed here.

Just as important, it left a decision environment — data marts, the model and simulation structures — in which inventory policy alternatives can be tested before they are tried.

Lessons learned

• Check feasibility per item before solving the whole. A cheap per-SKU pre-check turned an uninformative “infeasible” into a list of SKUs to look at.
• Model the square root. Pooled safety stock is what makes consolidation pay; a model with additive safety stock would never see the trade.
• Keep the business units in the model. Whole pallets and minimum orders change which plan wins for small SKUs, so they belong in the constraints rather than in a rounding step afterwards.
• Solve the direct form first, then convexify. The nonconvex model was quick to write and check; the cone form is the route to scale.

Conclusion

Where to keep stock across a distribution network is a joint decision: storage, sourcing and size interact through pooled safety stock, whole pallets and shared capacity. Writing them as one MIQCP in PySCIPOpt, guarded by a per-SKU feasibility check and fed by recency-weighted statistics, gave the business a quantified case for moving from conservative to demand-based inventory.

Reformulated as a second-order-cone program, the same model is the natural path from a feasibility study to a planning tool for the full network.

Limitations

• This was a feasibility study; the policy was proposed, not rolled out, and no operational savings are claimed.
• Single sourcing per demand point and a fixed service factor of three standard deviations are modelling choices, not business constraints.
• Demand statistics are recency-weighted history; there is no forecast model for promotions or new items.
• The live model solves 24 SKUs × 8 centres by plan enumeration and a linear integer program instead of the nonconvex model — exact only because the network is small — and its holding-cost settings are synthetic.

About the demo and confidentiality

The centres, SKUs, demand histories, lanes and costs in the embedded model are generated from a seed. No network structure, item master, cost table, conversion constant, storage weight or result from the real study appears in this post — the storage weight is written as w_s throughout; code is simplified and written for illustration.