Building a cart-picking optimizer with MILP, aisle-level dynamic programming and an async FastAPI service

In a cart picking system, empty boxes ride on a cart, orders are mapped to boxes, and a worker pushes the cart through the warehouse picking SKUs from many locations. Picking time per SKU is largely fixed, so productivity depends on how far workers walk and how long they wait for each other at busy locations. This is not a routing problem alone: which boxes share a cart changes the routes, and the routes change which grouping is good. Because one order can need several locations, it does not reduce to a travelling-salesman or vehicle-routing formulation.

In this post we walk through the optimization service we built. It ranks carts by a lexicographic objective — zones, then aisles, then travel time. In real-time mode it chooses the next cart with an exhaustive minimum-aisle-set search and lets a MILP re-partition (SciPy with HiGHS) improve the next three carts, keeping the result only if it is never worse. In whole-wave mode it packs every box into trips and orders their release so concurrent carts avoid each other's aisles. Both modes share a two-level dynamic program that routes a cart exactly on aisle geometry inferred from the travel-time matrix alone. The service runs asynchronously behind the site's execution system and answers by callback.

Solution overview

The optimizer is one Python service. The site's execution system sends chunked requests — boxes, box-type capacity, a location-to-location travel-time matrix, cart status and parameters — and receives the plan by HTTP callback. Inside, preparation fixes the mode and capacity rule and rebuilds the aisle layout, composition decides which boxes share a cart, a shared routing core scores every candidate cart, and in wave mode a release step orders the trips.

The numbered steps in the diagram:

1. A FastAPI service validates chunked requests with pydantic models, acknowledges them and runs the optimization in the background.
2. Each request is set to single-cart or whole-wave mode, a capacity rule is chosen — quantity per box type, count plus volume, or count — and aisle depths and positions are reconstructed from the travel-time matrix with NumPy.
3. An exhaustive minimum-aisle-set search and a greedy fill build three carts on a rolling horizon; a SciPy milp re-partition solved by HiGHS replaces them only if it wins on the same key, and only the first cart is released.
4. Orders are grouped by the zones they need, packed with few aisles, merged first-fit-decreasing, and improved by 1:1 order swaps between carts compiled with numba.
5. An in-aisle sweep and a between-aisle subset dynamic program give the exact route under aisle and zone rules; nearest neighbour with 2-opt, or-opt and relocate moves is the fallback.
6. A greedy release order, ranked with np.lexsort, keeps carts that run at the same time from sharing locations and aisles.
7. Header and detail rows with lead times are posted back with httpx and tenacity retries; the service ships as a Docker image through CI to a cloud container service.

Technology stack

Step 1: Rank carts with a lexicographic key

Every decision in the system — which aisle set to take, whether to keep a swap, whether to accept the MILP's answer — compares candidate carts, and every comparison uses one key with three levels: zones, then aisles, then walking time.

Zone-switch and aisle-switch penalties are inputs from the site. The zone level is counted only when the zone-switch penalty is positive, so a site without zone rules falls back to aisles, then travel.

def cart_key(boxes, route_time, zone_of, zones_matter):
    """Lexicographic cost of one cart. Python compares tuples element
    by element, so cart_key(a) < cart_key(b) means cart a is better."""
    aisles = {loc.aisle for box in boxes for loc in box.locations}
    zones = {zone_of[a] for a in aisles} if zones_matter else set()
    return (len(zones), len(aisles), route_time(boxes))

# Every search step compares carts with this key:
#   improvement moves are kept if  cart_key(new) <= cart_key(old)
#   the MILP answer is kept only if  cart_key(milp) < cart_key(constructive)

Simplified. route_time is the routing core from Step 5.

Why an ordering instead of one weighted sum? The floor's priorities are strict: no saving in metres justifies an extra zone. An ordering encodes that directly, and it turns every acceptance test in the system into a plain comparison.

Step 2: Accept work asynchronously and rebuild the layout from travel times

The site's execution system sends work as chunked requests and expects the plan back by callback. The service validates each request with pydantic, acknowledges it and runs the optimization as a background task; httpx posts the result and tenacity retries on transient failures. Preparation fixes the mode — one cart per request with the rest carried over, or a whole wave — and the capacity rule: quantity per box type, count plus volume (derived from the cart's cells and tiers), or count.

import httpx
from fastapi import BackgroundTasks, FastAPI
from pydantic import BaseModel
from tenacity import retry, stop_after_attempt, wait_exponential

app = FastAPI()


class PlanRequest(BaseModel):
    request_id: str
    mode: str                              # "next_cart" | "wave"
    boxes: list[Box]
    capacity: Capacity
    travel_time: list[list[float]]         # location x location, seconds
    callback_url: str


@retry(stop=stop_after_attempt(5), wait=wait_exponential(max=30))
def deliver(url: str, plan: dict) -> None:
    httpx.post(url, json=plan, timeout=30).raise_for_status()


def solve_and_reply(req: PlanRequest) -> None:
    deliver(req.callback_url, solve(req))  # steps 2-6, then the callback


@app.post("/plans")
async def submit(req: PlanRequest, tasks: BackgroundTasks):
    tasks.add_task(solve_and_reply, req)
    return {"accepted": req.request_id}

Simplified. Field names are generic; chunk handling, authentication and logging are omitted.

The request carries no floor map, only travel times. The routing core needs each location's aisle, its depth, and where each aisle meets the cross aisle, so the service infers that geometry from the matrix itself.

corridor model    d(p, q) = h(p) + | x(a_p) − x(a_q) | + h(q)        p, q in different aisles
                  h(p) = depth of location p in its aisle
                  x(a) = position of aisle a along the cross aisle

depths            fit the star metric  d(p, q) = h(p) + h(q)
positions         (p | q)_o = ½ · [ d(o, p) + d(o, q) − d(p, q) ] = min( x(a_p), x(a_q) )
                  the Gromov product against a point o at the end of the cross aisle
in-aisle times    rebuilt from the fitted depths, removing the matrix's rounding noise

Because the layout comes from data the site already maintains, a centre with a different floor plan needs no extra input — a large part of what made the inputs standard across centres.

Step 3: Choose the next cart from the smallest aisle set

In real-time mode the question “what should the next cart carry?” is asked repeatedly as orders arrive and answered within seconds. The model builds the next three carts on a rolling horizon and releases only the first; the other two are a lookahead.

Construction starts with an exhaustive search for the minimum aisle set: try every combination of one aisle, then of two, and so on, and stop at the first size whose boxes can fill a cart. Because sizes are tried in increasing order, that first size is the proven minimum — as long as the search finished inside its work budget. Among sets of that size, the one that opens fewer zones and has the shorter trial route wins.

from itertools import combinations

def min_aisle_set(boxes, capacity, max_k, budget, zone_count, trial_route):
    """Try aisle sets of size 1, 2, ... The first size whose boxes can fill
    a cart is the minimum, provided the work budget was not exhausted."""
    aisles = sorted({a for box in boxes for a in box.aisles})
    work = 0
    for k in range(1, max_k + 1):
        best_key, best_set = None, None
        for subset in combinations(aisles, k):
            work += 1
            if work > budget:
                return best_set, False             # best found, not proven
            chosen = set(subset)
            pool = [box for box in boxes if box.aisles <= chosen]
            if not capacity.can_fill(pool):
                continue
            key = (zone_count(chosen), trial_route(pool))
            if best_key is None or key < best_key:
                best_key, best_set = key, chosen
        if best_set is not None:
            return best_set, True                  # proven minimum size
    return None, True

Simplified. capacity.can_fill applies the active capacity rule; trial_route is a quick route estimate.

A greedy fill completes the cart from boxes inside the chosen aisles. Improvement passes — add, move, swap between carts, swap with the waiting pool — then run for a bounded number of rounds, and each move is kept only if it does not worsen the key.

Step 4: Re-partition the lookahead with a MILP that can only help

Greedy construction is fast but myopic, so the three carts are handed to an integer program that may move boxes between them. Boxes are first aggregated into signatures — the same aisle set, zone set and box type — so the model counts boxes per signature per cart instead of deciding every box separately.

signature   s = (aisle set A_s, zone set G_s, box type),   cnt_s boxes waiting
variables   n[s,c] ∈ ℤ₊      boxes of signature s in cart c,   c ∈ {1, 2, 3}
            y[a,c] ∈ {0,1}   aisle a open in cart c
            z[g,c] ∈ {0,1}   zone g open in cart c

min   −M · Σ n[s,c]  +  Σ y[a,c]  +  w_z · Σ z[g,c]  +  ε · (aisle and zone terms of cart 1)

s.t.  Σ_c n[s,c] ≤ cnt_s                               each box at most once
      n[s,c] ≤ ub_s · y[a,c]        ∀ a ∈ A_s          a box opens its aisles
      n[s,c] ≤ ub_s · z[g,c]        ∀ g ∈ G_s          … and its zones
      load_c(n) ≤ cap_c                                count, per type or volume
      load_1(n) ≥ cap_1                                cart 1 leaves full
      symmetry-breaking constraints between the interchangeable carts

The large weight M makes filling carts the first priority; open aisles come next, then zones weighted by w_z, and a small ε on cart 1's terms favours the cart that is actually released. Cart 1 must leave full. The model is built as a sparse constraint matrix and solved with scipy.optimize.milp, which calls HiGHS.

import numpy as np
from scipy.optimize import Bounds, LinearConstraint, milp

def repartition(sig_aisles, cnt, n_aisles, cap, C=3, big=1e3, node_limit=1000):
    """Variables: n[s, c] (boxes of signature s in cart c), then y[a, c]
    (aisle a open in cart c). Count capacity only."""
    S = len(sig_aisles)
    N = lambda s, c: s * C + c
    Y = lambda a, c: S * C + a * C + c
    nv = S * C + n_aisles * C
    rows, lo, hi = [], [], []
    def add(terms, l, h):
        r = np.zeros(nv)
        for i, v in terms:
            r[i] = v
        rows.append(r); lo.append(l); hi.append(h)
    for s, aisles in enumerate(sig_aisles):
        add([(N(s, c), 1) for c in range(C)], 0, cnt[s])          # supply
        for c in range(C):
            for a in aisles:                                       # opens its aisles
                add([(N(s, c), 1), (Y(a, c), -cnt[s])], -np.inf, 0)
    for c in range(C):                                             # cart 1 leaves full
        add([(N(s, c), 1) for s in range(S)], cap if c == 0 else 0, cap)
    cost = np.r_[np.full(S * C, -big), np.ones(n_aisles * C)]
    ub = np.r_[np.repeat(cnt, C), np.ones(n_aisles * C)]
    return milp(cost, integrality=np.ones(nv), bounds=Bounds(0, ub),
                constraints=LinearConstraint(np.array(rows), lo, hi),
                options={"node_limit": node_limit})

Simplified. Zones, per-type and volume capacity, the ε term and symmetry breaking are omitted.

The MILP's carts replace the constructive ones only if they win on the same lexicographic key from Step 1; otherwise the constructive answer stands, so the MILP can never make the plan worse. It runs under a node limit with a pinned solver version, and it is skipped when the number of signatures exceeds a cap, so large requests fall back to the constructive result.

Why aggregate into signatures? Boxes with the same aisle set, zone set and box type are interchangeable for this objective. One integer variable per signature and cart shrinks the model and removes a large source of symmetry.

Step 5: Route each cart exactly with a two-level aisle DP

The routing core is called for every candidate cart, so it has to be both exact and fast. It is split into two dynamic programs that match how a cart moves: along an aisle, and between aisles.

Inside an aisle, picks are ordered along the aisle axis and swept with forward and backward chains, allowing at most one reversal. The result is a small table: for every entry pick and exit pick, the cheapest time to visit all of the aisle's picks. Across aisles, a subset DP adds one aisle at a time, keeping the set of aisles visited and the pick the cart last left from as its state. A bitmask check forbids re-entering a zone the cart has already left, so each zone is entered once.

inner_j(e, x′)    cheapest sweep of aisle j: enter at pick e, visit every pick,
                  leave at pick x′, with at most one reversal

g({j}, x′)      = min over e        D(start, e) + inner_j(e, x′)
g(S ∪ {j}, x′)  = min over (x, e)   g(S, x) + D(x, e) + inner_j(e, x′)
                  allowed only if zone(j) is the current zone or one not yet entered
route time      = min over x        g(all aisles, x) + D(x, end)

INF = float("inf")

def relax(table, key, value):
    if value < table.get(key, INF):
        table[key] = value

def reenters_zone(mask, i, j, zone):
    """Stay in the current zone, or move to a zone not visited yet."""
    if zone[j] == zone[i]:
        return False
    return any(mask >> k & 1 and zone[k] == zone[j] for k in range(len(zone)))

def route_aisles(inner, D, start, end, zone):
    """inner[j][(e, x)] = time to sweep aisle j entering at pick e and leaving
    at pick x. State: (aisles visited as a bitmask, last aisle, exit pick)."""
    m = len(inner)
    g = [dict() for _ in range(1 << m)]
    for j in range(m):
        for (e, x), t in inner[j].items():
            relax(g[1 << j], (j, x), D[start][e] + t)
    for mask in range(1, 1 << m):                    # supersets come later
        for (i, x), base in g[mask].items():
            for j in range(m):
                if mask >> j & 1 or reenters_zone(mask, i, j, zone):
                    continue
                for (e, x2), t in inner[j].items():
                    relax(g[mask | 1 << j], (j, x2), base + D[x][e] + t)
    return min(t + D[x][end] for (j, x), t in g[(1 << m) - 1].items())

Simplified. Outer DP only; inner comes from the in-aisle sweep and D is the travel-time matrix.

The state space grows exponentially with the number of aisles, so it is capped. A zone-sweep variant of the DP avoids the blow-up when a cart spans many aisles, and if a cap is still exceeded the route falls back to nearest neighbour improved by 2-opt, or-opt and relocate moves, compiled with numba.

Step 6: Plan a whole wave and stagger its release

In wave mode the service receives every box at once. Orders are grouped by zone signature — the set of zones they need — and each group is packed into carts starting from the smallest aisle set that still fills a cart. Leftovers from different groups are merged first-fit-decreasing, and 1:1 order swaps between carts run while the key improves.

Bottleneck waiting is handled at release time. Trips are released in a greedy order that, at each step, takes the trip sharing the fewest locations, and then the fewest aisles, with the trips still running on other workers. np.lexsort ranks the candidates in one call.

import numpy as np

def release_order(trips, workers):
    """Greedy release: the next trip shares the fewest locations, then the
    fewest aisles, with the trips other workers are still running."""
    left, order = list(range(len(trips))), []
    while left:
        running = order[-(workers - 1):] if workers > 1 else []
        locs = set().union(*(trips[t].locations for t in running))
        aisles = set().union(*(trips[t].aisles for t in running))
        shared_locs = np.array([len(trips[t].locations & locs) for t in left])
        shared_aisles = np.array([len(trips[t].aisles & aisles) for t in left])
        k = np.lexsort((shared_aisles, shared_locs))[0]   # last key is primary
        order.append(left.pop(k))
    return order

Simplified. The trips still running are approximated by the last workers − 1 released.

Step 7: Keep it fast and reproducible

Two choices made the search affordable. First, the swap loop scores candidates with an approximate route cost and runs exact routing only once, on the final carts — tens of times faster than routing every candidate. Second, the inner loops — swap evaluation, routing and the fallback heuristics — are compiled with numba and warmed up at start-up, so the first real request does not pay the compile time.

import numpy as np
from numba import njit

@njit(cache=True)
def n_open(row):
    return np.count_nonzero(row)

@njit(cache=True)
def best_swap(cnt, need, size, cart_of, load, cap, max_pairs):
    """Best 1:1 order swap between two carts by the change in open aisles.
    cnt[c, a]: picks of cart c in aisle a; need[o, a]: picks of order o."""
    best, bo, bp, tried = 0, -1, -1, 0
    for o in range(need.shape[0]):
        for p in range(o + 1, need.shape[0]):
            c1, c2 = cart_of[o], cart_of[p]
            if c1 == c2:
                continue
            if (load[c1] - size[o] + size[p] > cap or
                    load[c2] - size[p] + size[o] > cap):
                continue
            tried += 1
            if tried > max_pairs:                  # a work limit, not a time limit
                return bo, bp, best
            before = n_open(cnt[c1]) + n_open(cnt[c2])
            after = (n_open(cnt[c1] - need[o] + need[p]) +
                     n_open(cnt[c2] - need[p] + need[o]))
            if after - before < best:
                best, bo, bp = after - before, o, p
    return bo, bp, best

Simplified. Only the aisle level of the key is shown; zones and an approximate travel term are compared the same way.

Every limit counts work, not time: a combination budget for the enumeration, a pair budget for swaps, a node limit for the MILP and a state cap for the DP, with the solver version pinned. Wall-clock limits would make the answer depend on machine load; work limits mean that replaying a request returns the same plan.

An earlier reinforcement-learning policy (PPO) was tried in the pipeline and later removed.

Results

In simulation against existing assignment and routing practice, the approach showed an average productivity improvement of around 20% within limited computation time. The gain comes from planning order allocation, movement and bottleneck waiting together rather than optimizing each cart's route in isolation.

The more durable contribution was the standardized structure. Order information, SKU locations, cart capacity, travel times and operating constraints are defined once as common inputs, and the layout itself is inferred from the travel times, so a centre with a different floor plugs into the same framework instead of getting a new model.

Lessons learned

• Encode the floor's priorities as an order. A lexicographic key made every comparison in the system unambiguous and kept travel savings from buying zone changes.
• Let exact methods prove what they can, and bound the rest. Enumeration proves the minimum aisle set and the DP gives exact routes; both carry work caps and fallbacks for the cases they cannot finish.
• Make the expensive optimizer earn its place. Accepting the MILP only when it wins on the same key let us add it without any risk of a worse plan.
• Approximate inside loops, be exact once. An approximate route cost in the swap loop and exact routing at the end made the search tens of times faster.
• Budget work, not time. Iteration, node and state limits plus a pinned solver made every plan reproducible.

Conclusion

Cart picking looks like routing but is an assignment problem with routing inside it. Ranking carts by zones, aisles and then travel, choosing them from the smallest aisle set, improving them with a MILP that can only help, and routing them exactly with a two-level DP turned that into a service that answers within seconds and gives the same answer every time.

The same structure — a strict objective, exact components where they are cheap, guarded improvement steps and inputs any centre can supply — is what lets the framework move to other picking layouts.

Limitations

• The productivity figure is a simulation estimate against existing assignment and routing practice, not a field A/B test.
• Picking time is treated as fixed, and congestion is handled at planning time by separating concurrent trips, not by re-routing carts on the floor.
• Routes are exact only within the model — the corridor-and-aisle geometry inferred from travel times and the aisle and zone rules — and beyond the DP's state caps the route falls back to heuristics.
• The live model is a JavaScript re-implementation on a generated two-zone warehouse with reduced search budgets; the embedded view runs without the MILP re-partition.

About the demo and confidentiality

The warehouse, aisles, orders, SKUs and locations in the embedded model are generated from a seed. No centre layout, shipper, site system, interface field, penalty value, operational count or replay result appears in this post; code is simplified and written for illustration.