Alternative Approaches to Assignment Problem

Question

I have a real world use case where I aim to optimize the assignment of trucks to orders of different sizes over a time window by minimizing some cost function (e.g. total distance traveled).

The below model (coded in Python with Google OR-Tools) illustrates my problem, and solves quickly on my toy dataset. My issue is that the real problem will contain approximately 1000 locations and 1000+ trucks, which causes the number of path variables (denoted r[i,t,k,k'] in the sample code) to explode to 1000 trucks * 15 time periods * choose(1000,2) distance matrix entries. Setting up this many variables and associated constraints in the MIP solver will be time prohibitive, and I have my doubts as to whether the solver (or a commercial solver) will be able to give a solution in a reasonable amount of time.

import pandas as pd
import numpy as np
from ortools.linear_solver import pywraplp

#CREATE DATA
distance_matrix = pd.DataFrame([(str(x),str(y)) for x in range(1,6) for y in range(1,6)], columns=['ORIGIN','DESTINATION'])
distance_matrix['DISTANCE'] = abs(distance_matrix.ORIGIN.astype(int)-distance_matrix.DESTINATION.astype(int))

orders = pd.DataFrame(
         {'ORDER': ['A', 'B', 'C'],
          'PICKUP_LOCATION': ['1', '4', '5'],
          'DELIVER_LOCATION': ['5', '1', '2'],
          'SIZE': [1, 2, 3],
          'PICKUP_EARLIEST': [0, 0, 0],
          'PICKUP_LATEST': [3, 3, 3],
          'DELIVER_EARLIEST': [4, 4, 4],
          'DELIVER_LATEST': [9, 9, 9]}
          )

trucks = pd.DataFrame(
         {'TRUCK': ['A', 'B', 'C'],
         'CAPACITY': [1, 2, 3],
         'CURRENT_LOCATION': ['4', '2', '3']}
         )


#DEFINE MODEL PARAMETERS
max_distance_per_time = 1 #MAX NUMBER OF UNITS A TRUCK CAN MOVE IN ONE TIME PERIOD
max_time = 9 #NUMBER OF TIME PERIODS TO OPTIMIZE OVER


#PRE-FILTER DISTANCE MATRIX TO ONLY ALLOWABLE DISTANCES OVER A TIME PERIOD
distance_matrix = distance_matrix.loc[distance_matrix.DISTANCE<=max_distance_per_time].reset_index(inplace=False, drop=True)

#DEFINE VECTOR OF POSSIBLE LOCATIONS FOR LOOPING THROUGH
all_locations = np.unique(np.append(distance_matrix.ORIGIN.unique(), distance_matrix.DESTINATION.unique()))


###BUILD MODEL
solver = pywraplp.Solver.CreateSolver('SCIP')


##CREATE DECISION VARIABLES
x={} #x[i,j,t]=1 if truck i picks up order j at time t, 0 otherwise
for truck_index, truck in trucks.iterrows():
    for order_index, order in orders.iterrows():
        for t in range(max_time+1):
            if t>=order.PICKUP_EARLIEST and t<=order.PICKUP_LATEST:
                x[truck_index, order_index, t] = solver.BoolVar(f'x[{truck_index},{order_index},{t}]')

y={} #y[i,j,t]=1 if truck i delivers order j at time t, 0 otherwise
for truck_index, truck in trucks.iterrows():
    for order_index, order in orders.iterrows():
        for t in range(max_time+1):
            if t>=order.DELIVER_EARLIEST and t<=order.DELIVER_LATEST:
                y[truck_index, order_index, t] = solver.BoolVar(f'y[{truck_index},{order_index},{t}]')

z={} #z[i,t,k]=1 if truck i is positioned at location k at time t, 0 otherwise
for truck_index, truck in trucks.iterrows():
    for t in range(max_time+1):
        for location_index, location in enumerate(all_locations):
            z[truck_index, t, location] = solver.BoolVar(f'z[{truck_index},{t},{location}]')

r = {} #r[i,t,k,k']=1 if z[i,t,k]*z[i,t-1,k']=1, 0 otherwise. REPRESENTS POINTS TRAVERSED BETWEEN t-1 AND t.
for truck_index, truck in trucks.iterrows():
    for t in range(1, max_time+1):
        for location_index, location in enumerate(all_locations):
            for od_index, od in distance_matrix[distance_matrix.DESTINATION==location].iterrows():
                r[truck_index, t, location, od.ORIGIN] = solver.BoolVar(f'r[{truck_index},{t},{location},{od.ORIGIN}]')
                solver.Add(r[truck_index, t, location, od.ORIGIN]-z[truck_index, t, location]<=0)
                solver.Add(r[truck_index, t, location, od.ORIGIN]-z[truck_index, t-1, od.ORIGIN]<=0)
                solver.Add(r[truck_index, t, location, od.ORIGIN] - z[truck_index, t, location] - z[truck_index, t-1, od.ORIGIN] >= -1)
                
                #CANNOT TRAVEL MORE THAN max_distance_per_time BETWEEN TIME STEPS (ADDING CONSTRAINT HERE TO AVOID LOOPING THROUGH AGAIN)
                solver.Add(r[truck_index, t, location, od.ORIGIN]*od.DISTANCE<=max_distance_per_time)
        
        #TRUCK MUST MAKE 1 MOVE EACH TIME PERIOD ALONG VALID PATH (INCLUDING MOVES TO SAME LOCATION)
        solver.Add(solver.Sum(r[truck_index, t, od.DESTINATION, od.ORIGIN] for od_index, od in distance_matrix.iterrows())==1)



##CREATE CONSTRAINTS
for truck_index, truck in trucks.iterrows():
    #TRUCK AT CURRENT_LOCATION AT TIME 0
    solver.Add(z[truck_index, 0, truck.CURRENT_LOCATION]==1)
    
    for t in range(max_time+1):
        #TOTAL SIZE LOADED ON TRUCK i AT TIME t <= CAPACITY[i]
        solver.Add(solver.Sum(order.SIZE*x[truck_index, order_index, time_index] for order_index, order in orders.iterrows() for time_index in range(order.PICKUP_EARLIEST, min(order.PICKUP_LATEST, t)+1) if t>=order.PICKUP_EARLIEST and t<=order.DELIVER_LATEST) -
                   solver.Sum(order.SIZE*y[truck_index, order_index, time_index] for order_index, order in orders.iterrows() for time_index in range(order.DELIVER_EARLIEST, min(order.DELIVER_LATEST, t)+1) if t>=order.DELIVER_EARLIEST and t<=order.DELIVER_LATEST)
                   <=truck.CAPACITY)
        
        #TRUCK MUST BE AT 1 AND ONLY 1 LOCATION AT EACH TIME
        solver.Add(solver.Sum(z[truck_index, t, location] for location in all_locations)==1)


for order_index, order in orders.iterrows():
    #ALL ORDERS MUST BE PICKED UP (EXACTLY ONCE)
    solver.Add(solver.Sum(x[truck_index, order_index, t] for truck_index in range(trucks.shape[0]) for t in range(order.PICKUP_EARLIEST, order.PICKUP_LATEST+1))==1)
    
    for truck_index, truck in trucks.iterrows():
        #ALL ORDERS PICKED UP ARE DELIVERED
        solver.Add(solver.Sum(x[truck_index, order_index, t] for t in range(order.PICKUP_EARLIEST, order.PICKUP_LATEST+1)) -
                   solver.Sum(y[truck_index, order_index, t] for t in range(order.DELIVER_EARLIEST, order.DELIVER_LATEST+1))
                   ==0)
        
        #TRUCK MUST BE AT ORIGIN AT TIME OF PICKUP
        for t in range(order.PICKUP_EARLIEST, order.PICKUP_LATEST+1):
            solver.Add(x[truck_index, order_index, t]-z[truck_index, t, order.PICKUP_LOCATION]<=0)
        
        #TRUCK MUST BE AT DESTINATION AT TIME OF DELIVERY
        for t in range(order.DELIVER_EARLIEST, order.DELIVER_LATEST+1):
            solver.Add(y[truck_index, order_index, t]-z[truck_index, t, order.DELIVER_LOCATION]<=0)




##DEFINE OBJECTIVE
objective_terms = []
for truck_index, truck in trucks.iterrows():
    for t in range(1,max_time+1):
        for od_index, od in distance_matrix.iterrows():
            objective_terms.append(r[truck_index, t, od.DESTINATION, od.ORIGIN]*od.DISTANCE)



#SOLVE THE MODEL
solver.Minimize(sum(objective_terms))
status = solver.Solve()


#PRINT THE SOLUTION
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    for truck_index, truck in trucks.iterrows():
        for t in range(1, max_time+1):
            for location_index, location in enumerate(all_locations):
                print_ind=0
                for prev_location in distance_matrix[distance_matrix.DESTINATION==location].ORIGIN:
                    if r[truck_index, t, location, prev_location].solution_value()>0.5:
                        print(f'Truck {truck.TRUCK} moved from {prev_location} to {location} at time {t}')
                
                for order_index, order in orders.iterrows():
                    if t>=order.PICKUP_EARLIEST and t<=order.PICKUP_LATEST and x[truck_index, order_index, t].solution_value() > 0.5 and order.PICKUP_LOCATION==location:
                        print(f'Truck {truck.TRUCK} picked up order {order.ORDER} at time {t} from location {order.PICKUP_LOCATION}')
                        print_ind+=1
                    if t>=order.DELIVER_EARLIEST and t<=order.DELIVER_LATEST and y[truck_index, order_index, t].solution_value() > 0.5 and order.DELIVER_LOCATION==location:
                         print(f'Truck {truck.TRUCK} delivered order {order.ORDER} at time {t} at location {order.DELIVER_LOCATION}')
                         print_ind+=1

Answer 1

If you are open to heuristic (approximate) solutions, one possibility is to generate a set of routes and then rewrite the problem to assign trucks and customers/orders to routes. In problems of this sort, a common tactic is to start with an initial set of routes created manually with an eye toward ensuring you have enough routes to make the problem feasible. The master problem is then relaxed from an IP to an LP and solved, and the dual variables are used to create objective functions in a subproblem that generates additional routes that could be used to improve the LP solution. At some point, you stop generating new routes (either because you cannot find any that look attractive, or because you hit a time limit, or because the master problem is getting too big), restore the integrality restrictions in the master problem and solve it.

Answer 2

It seems that you are trying to solve a pickup and delivery problem. MIP solvers are often not the best tool for that, and getting good results with them requires special models and tricks.

There are many heuristics that do the job. Since you are using OR-tools already, you may want to look at their pickup and delivery example. Other solvers provide direct support for this kind of problems too (LocalSolver and OptaPlanner do at least).

You can roll your own heuristic too, or look for better MIP models (column generation...), but that is no easy task.