Ship from store and save

As the world undergoes the digital transformations, accelerated by the global pandemics, omnichannel shopping experience becomes a new norm. Ship from Store turns out to be an important step on the path of establishing effective and convenient retail practice.

But the transformation brings its' own challenges and demands for new ways of business practices and business optimizations. Large retailers that operate hundreds of stores now face a question on how to choose the best way of fulfillment that is both cost effective and convenient. Traditionally, Order Management Systems (OMS) choose the optimal route for each order independently which often leads to sub-optimal results. In this article we want to present an alternative approach, leveraged by one of the largest US department stores, where we looked for an optimal fulfillment considering all in-flight orders over a time window.

By doing so, they were able to achieve 50% less order splits, double-digit EBIT improvement for optimized orders and multi-million annual savings.

Order cool down period

Important fact that makes this approach feasible is that typically eCommerce sites have certain cool down period while customers are allowed to make changes to their orders before those are pushed down to fulfillment. This allows to accumulate a larger number of orders for simultaneous optimization, without a need for immediate processing.

Example

Let's take a look at a small example. Assume that we have three incoming orders, the same item in each, but with different shipping destination. And those items are available in few stores, with different costs associated with each store and destination combination.

If we would follow "per order" optimization (in red), then the first customer order would be shipped from store A, 2nd - B and third - C; with a total cost of 2 + 19 + 18 = 39.

But if we try to optimize all three orders together, we can see a better option 1 - C, 2 - A, 3 - B (marked in green), that has cost of 7 + 15 + 6 = 28.

In this example we can see that leveraging global optimization allowed to introduce some savings, but when we multiply it over thousands of orders this might make a noticeable increase in business efficiency. Let's talk about how we can implement that.

Cost

While mathematical approach requires us to present a cost of shipping as a single number, it does not have to be an exact monetary cost of shipment, but can be a derived score that combines many business incentives associated with the shipment:

• The actual shipping cost C 𝑖,𝑗, which ultimately can be a non-linear function depending on total weight of items in shipment and distance,
• Projected time on the shelf in the store, 𝑃𝑖,𝑗
• Promised delivery time 𝐸𝑖,𝑗
• Count weeks before the item expires or becomes obsolete
• Number of total shipments in order

And a certain weight can be associated with importance of each goal, so for example it could look like:

𝐶𝑜𝑠𝑡 𝑖,𝑗 = C 𝑖,𝑗∗𝑊𝑐+𝑃𝑖,𝑗∗𝑊𝑝+𝑆(𝐸𝑖,𝑗)∗𝑊𝑒, where S(x) is a sigmoid function

Brute Force Combinability

Naive approach could be to just try every possible combination and choose one with the least cost. Unfortunately but predictably that becomes quite expensive as soon as number of orders grow. With the example below (3 orders and 3 stores), we only had to compare 27 combinations, but let's take a look how this number grows.

Even if we could check 1 million combinations in a second, it would take almost 3 hours to check all solutions for 10 stores and 10 customers. For 10 stores and 20 customers, it could need more than 3 million years with this approach.

If we were lucky as species we could probably calculate 14 orders fulfilling in 50 stores before the Sun burns out.

MILP optimization

Luckily there are some optimizations we can leverage. Let's try Mixed-Integer Linear Programming (MILP); this approach appeared in the middle of 20th century, and found it's place in theoretical and practical computer science.

This solution applies to problems when we need to find optimal distribution and minimize some loss, but the distribution parameters can take only allowed values, like we can't ship half of item from store A and another half from store B. Mathematically it is known to be a NP-complete problem; and the specific algorithm we will utilize is called branch-and-cut, you can find some more details below; essentially the idea behind it is that if we know that particular branch of combinations is known to be worse than current best we can skip it entirely.

If we take a look back at our example:

And assume we have already found across solution 1-C, 2-A, 3-B with cost of 29, we can skip checking any combinations with 2-B, 3-C and 2-C, 3-B, as we know that they can't be better than one we already have.

Using Solvers

There's a number of available MILP solvers both free and commercial. Let's take a look how we can use free CBC solver that is available in the pywraplp module of Google OR-TOOLS package.

First let's load the library and initialize the solver. We also limit max calculation time to 5 minutes. For such problems, it is possible that checking all combinations might actually take years even with the solver, so we need to put some reasonable limit to get 'good enough' estimate.

problem_type = pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING
solver = pywraplp.Solver('Ship From Store: optimal fulfillment', problem_type)
solver.SetTimeLimit(5 * 60 * 1000)

Now, let's define contraints and target for optimization:

X = {}
for i in range(order_count):
    for j in range(store_count):
        for k in range(catalog_count):
            X[i, j, k] = s.IntVar(0, int(orders[i, k]), 'x_{}_{}_{}'.format(i, j, k))

# Create u_ij variables
U = {}
for i in range(order_count):
    for j in range(store_count):
        U[i, j] = s.IntVar(0, 1, 'u_{}_{}'.format(i, j))

# Add constraint: order must be satisfied
for i in range(order_count):
    for k in range(catalog_count):
        s.Add(s.Sum(X[i, j, k] for j in range(store_count)) == orders[i, k])

# Add constraint: do not ship more SKUs from store than available inventory in this store
for j in range(store_count):
    for k in range(catalog_count):
        s.Add(s.Sum(X[i, j, k] for i in range(order_count)) <= availability[j, k])

# Big const
T = catalog_count * np.max(orders)

# Add constraint: set u_ij to one if anything is shipped from store j to client i
for i in range(order_count):
    for j in range(store_count):
        s.Add(s.Sum(X[i, j, k] for k in range(catalog_count)) <= T * U[i, j])

# Set target function
s.Minimize(s.Sum(cost[i, j] * U[i, j] for i in range(order_count) for j in range(store_count)))

Finally, let's run the solver and print the results

res = s.Solve()

sol = {}
# Extract the solution from solver
for k, x in X.items():
    if SolVal(x) > 0:
        sol[k] = SolVal(x)
        
s.Iterations(), s.Objective().Value(), res, sol

    
def solve_problem_cbc(order_count, store_count, catalog_count, orders, availability, cost):
    problem_desc = 'SFS simple model'
    
    # pywraplp supports several optimizers, proprietary and open-source
    # We should specify which solver are we going to use
    problem_type = pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING

    # Initializing solver
    s = pywraplp.Solver(problem_desc, problem_type)

    # Create x_ijk variables from model description above
    X = {}
    for i in range(order_count):
        for j in range(store_count):
            for k in range(catalog_count):
                X[i, j, k] = s.IntVar(0, int(orders[i, k]), 'x_{}_{}_{}'.format(i, j, k))

    # Create u_ij variables
    U = {}
    for i in range(order_count):
        for j in range(store_count):
            U[i, j] = s.IntVar(0, 1, 'u_{}_{}'.format(i, j))
    
    # Add constraint: order must be satisfied
    for i in range(order_count):
        for k in range(catalog_count):
            s.Add(s.Sum(X[i, j, k] for j in range(store_count)) == orders[i, k])

    # Add constraint: do not ship more SKUs from store than available inventory in this store
    for j in range(store_count):
        for k in range(catalog_count):
            s.Add(s.Sum(X[i, j, k] for i in range(order_count)) <= availability[j, k])
    
    # Big const
    T = catalog_count * np.max(orders)
    
    # Add constraint: set u_ij to one if anything is shipped from store j to client i
    for i in range(order_count):
        for j in range(store_count):
            s.Add(s.Sum(X[i, j, k] for k in range(catalog_count)) <= T * U[i, j])
    
    # Set target function
    s.Minimize(s.Sum(cost[i, j] * U[i, j] for i in range(order_count) for j in range(store_count)))

    # For such kind of problems it's a common situation that solution may take years, even with solvers
    # So it's a good idea to set a time limit for solution: when the time runs out 
    # the best found solution returned, without any guarantees that it's a best one existing.
    # If time runs out, it's also possible to obtain best found estimate for the optimal solution
    # Here we set time limit to 300 seconds
    s.SetTimeLimit(300000)
    
    # Solve the problem. "res" variable contains return code for the solution. Description of codes
    # can be found here 
    # https://google.github.io/or-tools/java/enumcom_1_1google_1_1ortools_1_1linearsolver_1_1MPSolver_1_1ResultStatus.html
    res = s.Solve()
    
    sol = {}
    # Extract the solution from solver
    for k, x in X.items():
        if SolVal(x) > 0:
            sol[k] = SolVal(x)
            
    return s.Iterations(), s.Objective().Value(), res, sol    
    

Execution time

Now let's see how much time does it take. First let's start with 20 orders, 5 stores and 10 items, and do a multiple random runs. All of samples are calculated under 3 seconds, with mean execution time about 0.5 seconds.

But if we double the size and try allocating 40 orders, 10 stores and 20 items, solver can't complete full search in 5 minutes anymore.

Commercial Solvers

So the next step is to see what commercial solvers can achieve. For MILPs, the best two options are Gurobi and Ibm's CPLEX. Here's some summary of Gurobi execution.

Now the execution time distribution for second problem looks much better, with mean execution time under 50 seconds and 90% completing under 90 seconds:

Also the worst case is still solved under 155 seconds, less than three minutes.

import gurobipy as gp
from gurobipy import GRB

def solve_problem_grb_core(model, order_count, store_count, sku_count, orders, availability, cost):
    # Create x_ijk variables
    x_vars = model.addVars(range(order_count), range(store_count), range(sku_count),
                       vtype=GRB.BINARY,
                       name="x_assign")
    
    # Create u_ij variables
    u_vars = model.addVars(range(order_count), range(store_count),
                       vtype=GRB.BINARY,
                       name="u_assign")
    
    # Tell the optimizer that we're looking for minimum-cost solution
    model.modelSense = GRB.MINIMIZE
    
    # Add constraint: order must be satisfied
    for i in range(order_count):
        for k in range(sku_count):
            model.addConstr(sum(x_vars[i, j, k] for j in range(store_count)) == orders[i][k])

    # Add constraint: do not ship more SKUs from store than available inventory in this store
    for j in range(store_count):
        for k in range(sku_count):
            model.addConstr(sum(x_vars[i, j, k] for i in range(order_count)) <= availability[j][k])

    # Big const
    T = sku_count * np.max(orders)
    
    # Add constraint: set u_ij to one if anything is shipped from store j to client i
    for i in range(order_count):
        for j in range(store_count):
            model.addConstr(sum(x_vars[i, j, k] for k in range(sku_count)) <= T * u_vars[i, j])
            
    # Set formula for the target function
    model.setObjective(sum(cost[i][j] * u_vars[i, j] for i in range(order_count) for j in range(store_count)))
    
    # Run optimization
    model.optimize()
    
    # Return optimization status code
    # Reference can be found here
    # https://www.gurobi.com/documentation/9.1/refman/optimization_status_codes.html
    return model.Status, model.Runtime

def solve_problem_grb(param_I, param_J, param_K, mat_D, mat_S, mat_C):
    model = gp.Model("sfs")
    return solve_problem_grb_core(model, param_I, param_J, param_K, mat_D, mat_S, mat_C)
    

How can we improve it further?

The performance of the solver heavily depends on how accurate lower and upper bounds are estimated, and also on the order of branch evaluation: the sooner solver comes across a solution close to optimal, the more effective cutting it is able to leverage, the less combinations it has to evaluate, hence getting optimal solution faster.

To direct that process solvers use heuristics, where they spend some time evaluating data in order to find most promising order of evaluation, some of it before all calculations, some along the way. By default Gurobi allocates 5% of time to heuristics, but if we override and increase it to 10% we can see much better performance.

m = gp.Model("sfs")

# Default value of m.Params.Heuristics is 0.05
m.Params.Heuristics = 0.1

Now average execution time increased by 1.5 seconds, but 90% decreased by 5 seconds, and the worst case decreased by almost 10 seconds.

Early Stopping

As we can see majority of cases are solved way under the limit, let's see how much less optimal solution for the worst cases we get if we stop early, maybe we don't even need to solve all cases optimally.

So let's try to set max limit to 1 minute, instead of 5. Now we can see that only 20% of cases were terminated due to time out, and the solutions found within the time limit were just 2% to 3% higher than the optimal.

Next Optimization Steps

We need to mention that further accuracy improvement can be achieved with custom lower boundary calculation procedure, but will keep details out of scope of this article. The same applies to heuristic determination: with custom logic based on the domain specifics it's possible to speed up the solver few times.

Conclusion

We have demonstrated how using MILP techniques we were able to introduce cross order optimization into ecommerce ship-from-store order fulfillment. Free and commercial software is available with ready to use solutions, and tuning of those allows to achieve significant performance improvements.