Building SKU allocation and injection-order engines for a QPS with knapsack DP, simulation and Tabu Search

In a QPS, picking stations sit around a ring conveyor. Boxes are injected at fixed intervals, travel the ring and stop at every station that holds one of their SKUs, where a worker picks into them. Two decisions dominate throughput. Where each SKU sits decides how many stations a box must visit and how evenly the work spreads. The order in which boxes enter decides whether stations queue up — and when a station's queue is full, the box is bypassed and has to go round again. Operations had relied on experience and manual judgment, and the two decisions run on different clocks: allocation once per batch, injection order again and again while the batch runs.

In this post we walk through the two engines we built. The allocation engine measures SKU co-ordering with a sparse matrix product, seeds each station with an anchor SKU and fills stations in rounds with an exact two-dimensional 0/1 knapsack DP, for several station counts. The reorder engine scores an injection order with a discrete-event simulation of the ring and improves it with random multi-start and Tabu Search, re-planning on a rolling horizon with released boxes fixed. Both run behind a FastAPI service that reads its inputs from a Redis cache and answers by callback. Deployed in regular operations at the Dongtan logistics centre, the system delivered a 20% productivity improvement.

Solution overview

Both engines live in one Python service with one endpoint per engine. A call is acknowledged at once with HTTP 202; the engine runs as a background task on a snapshot read from a Redis cache and posts a compressed JSON result back to the operating system. The allocation engine runs once per batch, the reorder engine repeatedly during operation, and both share one sequencing core: a ring-conveyor simulation with a search on top.

The numbered steps in the diagram:

1. A FastAPI endpoint per engine answers 202 and runs the work as a background task; masters, inventory, order lines and task status are read from a Redis cache.
2. pandas builds order × SKU and order × zone workloads, and a SciPy sparse product XᵀX gives SKU co-order counts.
3. A deterministic multi-step heuristic places every SKU on one station, with an exact two-dimensional 0/1 knapsack DP compiled with numba at its core.
4. A discrete-event simulation of the ring scores an injection order by makespan and bypasses.
5. Random starts seed the search; Tabu Search then evaluates sampled position swaps in parallel with numba prange.
6. During operation released boxes stay fixed and the rest are re-sequenced batch by batch on a rolling horizon, then grouped into box-type runs.
7. Results are serialized to JSON, gzip-compressed, base64-encoded and POSTed back; the service ships as a Docker image.

Technology stack

Step 1: Serve two engines behind an asynchronous API

Allocation runs when stations are set up for a batch; reordering runs again and again while the batch is on the floor. Both are endpoints of one FastAPI service. A call is acknowledged immediately with HTTP 202 and the engine runs as a background task. Its inputs — station and zone masters, inventory, order lines with SKU, quantity and box type, task status, batch definitions and station priorities — are read from a Redis cache that the service treats as read-only, so each call works on the operating system's current snapshot.

import base64, gzip, json
import httpx, redis
from fastapi import BackgroundTasks, FastAPI

app = FastAPI()
cache = redis.Redis(decode_responses=True)       # read-only snapshot source

def pack(result: dict) -> str:
    return base64.b64encode(gzip.compress(json.dumps(result).encode())).decode()

def run(engine, req):
    snap = read_snapshot(cache, req.scope)       # masters, stock, orders, task status
    result = engine(snap, seed=req.seed)         # fixed seed: same input, same answer
    httpx.post(req.callback_url, json={"job": req.job_id, "result": pack(result)},
               timeout=60)

@app.post("/jobs/allocation", status_code=202)
async def post_allocation(req: JobRequest, tasks: BackgroundTasks):
    tasks.add_task(run, allocate, req)           # once per batch
    return {"job": req.job_id}

@app.post("/jobs/reorder", status_code=202)
async def post_reorder(req: JobRequest, tasks: BackgroundTasks):
    tasks.add_task(run, resequence, req)         # repeatedly during operation
    return {"job": req.job_id}

Simplified. Paths, field names and the snapshot reader are generic; authentication and error handling are omitted.

Results are serialized to JSON, gzip-compressed and base64-encoded before the callback, which keeps long sequences small on the wire. Seeds are fixed, so the same snapshot produces the same answer. The service is containerized with Docker and goes through CI with static analysis.

Step 2: Measure which SKUs are ordered together

A box visits every station that holds one of its SKUs, so an allocation is good when SKUs ordered together sit together and when picking load spreads evenly. Both quantities come from the order lines. A binary order × SKU matrix X gives each SKU's load — the number of distinct orders containing it — and, in one sparse product, every pair's co-order count.

X ∈ {0,1}^(orders × SKUs)      X_os = 1 if order o contains SKU s
load_s  = Σ_o X_os               distinct orders containing s
C       = XᵀX                    C_ij = orders containing both i and j      keep C_ij ≥ 2
J_ij    = C_ij / (load_i + load_j − C_ij)
aff_ij  = J_ij · ln(1 + C_ij)

target load per station    T = Σ_s load_s / Z            Z = stations opened

Jaccard overlap alone favours rare pairs that happen to coincide; raw counts favour SKUs that are simply popular. Their product rewards pairs that are both consistent and frequent, and pairs seen only once are dropped as noise.

import numpy as np
import pandas as pd
from scipy import sparse

def affinity(lines: pd.DataFrame, min_count=2):
    """lines: one row per (order, sku). Returns the load of each SKU and a
    sparse affinity matrix aff_ij = Jaccard_ij * ln(1 + C_ij)."""
    o = lines["order"].astype("category").cat.codes.to_numpy()
    s = lines["sku"].astype("category").cat.codes.to_numpy()
    X = sparse.csr_matrix((np.ones(len(o)), (o, s)))
    X.data[:] = 1                                   # an order counts once per SKU
    C = (X.T @ X).tocoo()                           # co-order counts
    load = np.asarray(X.sum(axis=0)).ravel()        # distinct orders per SKU
    keep = (C.row != C.col) & (C.data >= min_count)
    i, j, c = C.row[keep], C.col[keep], C.data[keep]
    jaccard = c / (load[i] + load[j] - c)
    aff = sparse.csr_matrix((jaccard * np.log1p(c), (i, j)), shape=C.shape)
    return load, aff

Simplified. Column names are generic; zone workloads are built the same way.

Step 3: Allocate SKUs to stations with rounds of an exact knapsack

Each SKU goes to exactly one station. Instead of one monolithic model, the allocation is built in deterministic stages, each enforcing cell capacity and the pallet-zone rule:

• Fixed SKUs that operations has already placed are pre-assigned.
• Top-decile SKUs by load go to pallet cells.
• Anchors: one seed SKU per station, strongly connected to other SKUs but weakly connected to the other anchors, so stations start apart.
• Knapsack rounds: stations are filled to 20%, 40%, 60%, 80% and 100% in turn, each step an exact knapsack.
• Polish: greedy placement of what is left, affinity-driven moves, and 1:1 swaps between heavy and light stations to balance load.

for r in (0.2, 0.4, 0.6, 0.8, 1.0), for each station z in turn:

  max   Σ_s (load_s + α_r · aff(s, z)) · x_s           x_s ∈ {0,1}, unassigned SKUs only
  s.t.  Σ_s load_s · x_s ≤ r · T − load_z               load budget
        Σ_s x_s        ≤ r · cells_z − used_z            cell budget

  aff(s, z) = affinity of s to the SKUs station z already holds;  α_r grows as stations fill

For one station in one round, the knapsack picks the SKUs that maximize load plus α times their affinity to what the station already holds, within two budgets: the station's share of the target load and its share of cells. Loads are integer order counts, so a two-dimensional dynamic program solves each knapsack exactly, and numba compiles it.

import numpy as np
from numba import njit

@njit(cache=True)
def knapsack_2d(value, weight, L, K):
    """Exact 0/1 knapsack with two budgets: sum of weight <= L, count <= K.
    value = load + alpha * affinity to the station; weight = integer load."""
    n = len(value)
    best = np.zeros((L + 1, K + 1))
    take = np.zeros((n, L + 1, K + 1), dtype=np.bool_)
    for i in range(n):
        w = weight[i]
        for l in range(L, w - 1, -1):              # backwards: each SKU once
            for k in range(K, 0, -1):
                cand = best[l - w, k - 1] + value[i]
                if cand > best[l, k]:
                    best[l, k] = cand
                    take[i, l, k] = True
    pick = np.zeros(n, dtype=np.bool_)
    l, k = L, K
    for i in range(n - 1, -1, -1):                 # walk the choices back
        if take[i, l, k]:
            pick[i] = True
            l -= weight[i]
            k -= 1
    return pick

Simplified. Budgets are the round's load and cell allowances for one station.

Why rounds? Filling one station to 100% first would let it take the best SKUs and leave the last station with leftovers. Raising every station's allowance together lets all of them compete for good SKUs.

The whole build is repeated for several station counts, so the floor sees alternatives — open one more station or not — each with its predicted duration and man-hours from the simulation in the next step.

Step 4: Score an injection order with a discrete-event simulation

The injection order is a permutation of boxes on fixed release slots. Its quality depends on interactions no simple proxy captures — queues, a single worker per station, bypasses when a queue is full — so the objective is a simulation. A box arriving at a station with nothing to pick moves on after a hop. If it has picks and the queue has room, it waits for the worker, whose service time is proportional to the number of distinct SKUs to pick. If the queue is full, it is bypassed and must come round again.

decision    π = order of boxes on fixed injection slots (slot k at time k·Δ)
service     t_pick(b, s) = τ · (distinct SKUs of box b at station s)
movement    t_hop per station;  one worker and a FIFO queue of length q per station
bypass      queue full on arrival  →  the box recirculates and tries again next lap

objective   min_π   makespan(π) + bypasses(π)

import heapq

def simulate(order, stops, n_st, slot, hop, per_sku, qlen):
    """Ring conveyor, event by event. stops[b] = {station: distinct SKUs}.
    One worker and a finite FIFO queue per station. Returns (makespan, bypasses)."""
    free_at = [0.0] * n_st                      # when each worker is next free
    queued = [[] for _ in range(n_st)]          # finish times of boxes at a station
    done = {b: set() for b in order}
    events = [(k * slot, b, 0) for k, b in enumerate(order)]   # (time, box, station)
    heapq.heapify(events)
    makespan, bypasses = 0.0, 0
    while events:
        t, b, s = heapq.heappop(events)
        if s in stops[b] and s not in done[b]:
            queued[s] = [f for f in queued[s] if f > t]
            if len(queued[s]) < qlen:           # join the queue, get picked
                free_at[s] = max(t, free_at[s]) + per_sku * stops[b][s]
                queued[s].append(free_at[s])
                done[b].add(s)
                t = free_at[s]
                makespan = max(makespan, t)
            else:
                bypasses += 1                   # queue full: come round again
        if len(done[b]) < len(stops[b]):
            heapq.heappush(events, (t + hop, b, (s + 1) % n_st))
    return makespan, bypasses

Simplified. In production the simulation is compiled with numba.

In the live model the parameters are ten seconds per SKU, three seconds per hop and five queue slots per station.

Step 5: Search the order with multi-start and Tabu Search

The search is deliberately plain, because every candidate costs a full simulation. Many random permutations are simulated first and the best is kept, stopping early once new starts stop helping. Tabu Search then samples position-pair swaps around the current order, simulates each one and moves to the best even if it is worse; a short tabu tenure forbids undoing recent swaps, which lets the search walk out of local minima.

import numpy as np

def search_order(boxes, score, evaluate, rng, starts, patience, iters, neigh, tenure):
    """score(order) and evaluate(order, moves) run the simulation and
    return makespan + bypasses."""
    best, best_val, stale = None, np.inf, 0
    for _ in range(starts):                            # random multi-start
        cand = rng.permutation(boxes)
        val = score(cand)
        if val < best_val:
            best, best_val, stale = cand, val, 0
        elif (stale := stale + 1) >= patience:         # early stop
            break
    cur, tabu = best.copy(), {}
    for it in range(iters):                            # Tabu Search
        pairs = {tuple(sorted(rng.choice(len(cur), 2, replace=False)))
                 for _ in range(neigh)}
        moves = np.array([m for m in pairs if tabu.get(m, -1) < it])
        if len(moves) == 0:
            continue
        vals = evaluate(cur, moves)                    # parallel kernel below
        k = int(np.argmin(vals))
        i, j = moves[k]
        cur[[i, j]] = cur[[j, i]]                      # take it, even if worse
        tabu[(i, j)] = it + tenure
        if vals[k] < best_val:
            best, best_val = cur.copy(), vals[k]
    return best, best_val

Simplified. Aspiration rules and budgets are omitted.

Swap evaluation is the hot loop, so it is a numba kernel that simulates the sampled swaps in parallel with prange.

import numpy as np
from numba import njit, prange

@njit(parallel=True, cache=True)
def evaluate_swaps(order, moves, stops, params):
    """One full ring simulation per sampled swap, spread over threads."""
    out = np.empty(len(moves))
    for k in prange(len(moves)):
        cand = order.copy()
        i, j = moves[k, 0], moves[k, 1]
        cand[i], cand[j] = order[j], order[i]
        makespan, bypasses = simulate_nb(cand, stops, params)
        out[k] = makespan + bypasses
    return out

Simplified. simulate_nb is the Step 4 simulation over NumPy arrays; evaluate in the search is this kernel with the batch data bound.

The order handed to the floor is not the permutation itself but the order in which boxes first start picking when the best permutation is simulated. Two engineering details mattered in production: numba kernels are warmed up on the first call so compile time is not paid inside a request, and a thread-safe parallel layer (TBB) was needed once several requests could run concurrently.

Step 6: Re-plan during operation on a rolling horizon

During operation the floor changes, so the reorder engine re-sequences what has not been released yet. Released boxes form a fixed prefix; the remaining boxes are re-sequenced batch by batch, each batch searched with everything before it fixed. Calling the engine again never moves a box that is already on its way.

def reorder(released, batches, search_batch):
    """Re-sequence what has not been released. Released boxes keep their
    slots; each batch is searched with everything before it fixed."""
    fixed = list(released)
    for batch in batches:                               # batch definitions, in order
        tail = search_batch(prefix=fixed, boxes=batch)  # simulation + Tabu
        fixed += runs_by_box_type(tail)
    return fixed[len(released):]


def runs_by_box_type(order):
    """Contiguous runs per box type, in order of first appearance; the
    searched order is kept within a type. Not re-simulated."""
    first = {}
    for k, box in enumerate(order):
        first.setdefault(box.box_type, k)
    return sorted(order, key=lambda box: first[box.box_type])

Simplified. search_batch runs Steps 4–5 for one batch behind the fixed prefix.

A final post-process groups boxes of the same box type into contiguous runs, because frequent alternation between box types slows injection on the floor. It is not re-simulated, so it can give back a little of the makespan the search found — a deliberate trade for a sequence the floor can work with.

Results

The system demonstrated a 20% productivity improvement at the Dongtan logistics centre and was deployed into regular field operations. The model design was adjusted around field requirements — response speed, stability, operational constraints and integration feasibility — and outputs were validated against actual productivity, not only against the simulator.

Alongside the engines, an LLM-agent-based AI Floor Supervisor interprets real-time operational conditions, monitors workload and bottleneck signals, and recommends task allocation and priority adjustments to field managers; it is outside the scope of this post. The result supported a planned rollout to a second centre and the wider logistics AI roadmap, including cart picking optimization.

Lessons learned

• Split decisions by their clock. Allocation once per batch and reordering during operation became two engines on one simulator; treating them as one problem would have made both slow.
• Solve clean subproblems exactly. Each knapsack is small and well-defined, so it is solved exactly, while the overall allocation stays a fast, deterministic heuristic.
• Make the simulator the objective. Queues and bypasses are what the floor feels, and they only show up when the ring is simulated.
• Fix what is already committed. A rolling horizon with a fixed prefix made it safe to call the reorder engine repeatedly while a batch runs.
• Plan for concurrency early. JIT warm-up and a thread-safe parallel layer were needed as soon as requests overlapped.

Conclusion

A ring-conveyor picking system has two levers on two clocks. Measuring co-ordering with a sparse product, filling stations with an exact knapsack in rounds, scoring injection orders with a discrete-event simulation and improving them with Tabu Search gave two engines the operation calls in its own rhythm — and a 20% productivity improvement at the Dongtan centre.

The pattern — exact optimization where the structure is clean, simulation plus search where interactions dominate, and an asynchronous service that fits the operating system's cycle — is what moved the work from offline analysis into regular operations.

Limitations

• The reorder simulation approximates boxes already on the ring by re-injecting them from the start of the horizon; it does not track their true position.
• The objective has no due dates or order priorities — only completion time and bypasses.
• Box-type grouping is applied after the search and is not re-evaluated against the simulation.
• The live model is a Python port on a generated batch of 40 boxes, with a reduced search budget and a simplified station-selection rule; the AI Floor Supervisor is not reproduced.

About the demo and confidentiality

The orders, SKUs, stations and box types in the embedded model are generated from a seed. No equipment ID, zone code, cache key, endpoint, port, batch volume or production setting appears in this post, and the AI Floor Supervisor's internals are not described; code is simplified and written for illustration.