1
$\begingroup$

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
$\endgroup$

2 Answers 2

1
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Please correct me if I'm wrong, but this appears to be effectively what is happening behind the scenes in the OR-Tools VRP example that @Alnurn linked. I'm accepting this answer since it's more general and gives me a path if the OR-Tools solution doesn't fit my needs out of the box. $\endgroup$
    – Jordan C
    Mar 9, 2023 at 18:07
  • $\begingroup$ Not at all. The routing library in or-tools is not based on column generation or MIP solvers. $\endgroup$ Mar 9, 2023 at 18:32
  • $\begingroup$ @LaurentPerron in that case, do you have a recommendation on which path I should go down given the size of my problem? One detail I omitted in the original question is that eventually I want to create dummy orders and investigate the shadow prices associated with them. Assuming the or-tools routing library will allow that, are there other factors I should consider when choosing between the two approaches? $\endgroup$
    – Jordan C
    Mar 9, 2023 at 19:27
  • $\begingroup$ The routing library will not give you shadow price. Not, if you think you can can write the column generation approach, and interpret correctly the dual information to get the shadow price, go for it. Personally, I would not go that way. $\endgroup$ Mar 11, 2023 at 10:13
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.