@ecotner Solution #2

SOLUTION

ROUTING SUCCESS
30.0 sec elapsed
Depot: [0]
Customers: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100]
Force Return: [101, 102]
Dropped Nodes:  [] 
Dropped distance: []
Dropped TWs: []
Routes:
        0: [0, 93, 20, 43, 82, 10, 41, 81, 69, 24, 6, 18, 39, 11, 96, 17, 61, 60, 90, 64, 54, 74, 48, 75, 97, 0]
        1: [0, 100, 32, 36, 89, 15, 67, 80, 79, 37, 55, 71, 65, 19, 8, 98, 76, 12, 28, 46, 30, 57, 40, 4, 44, 26, 66, 91, 13, 22, 63, 85, 5, 0]
        2: [0, 86, 52, 92, 95, 23, 45, 72, 29, 70, 35, 83, 62, 56, 88, 33, 87, 1, 27, 59, 3, 78, 53, 0]
        3: [0, 101, 0]
        4: [0, 49, 14, 68, 42, 50, 7, 31, 94, 73, 9, 47, 58, 51, 84, 16, 99, 77, 34, 25, 38, 21, 2, 102, 0]

solver.py

import time, re

from colorama import Fore, Back, Style
import numpy as np
from scipy.spatial.distance import pdist, squareform
from ortools.constraint_solver import pywrapcp, routing_enums_pb2

RNG_SEED = 42
PRECISION = 1_000

NUM_STOPS = 100
MEAN_DIST = 35 # "Kilometers"
NUM_VEHICLES = 5
NUM_MUST_RETURN = 2
VEHICLE_SPEED = 65 # "KPH"

MAX_TIME = 7 # "Hours" that we serve customers in.
EXTRA_TIME = 1 # Time at EOD that we wont schedule customers in.
MIN_WINDOW = 1 # Minimum size of customer windows in hours.

DROP_PENALTY = 10000 * PRECISION

SOLVER_INFINTIY = int(2 ** 63 - 1)

NODES_REAL = NUM_STOPS + 1
NODES_TOTAL = NUM_STOPS + NUM_MUST_RETURN + 1

def color_substr(string, substr, color):
    start = getattr(Fore, color)
    end = Fore.RESET
    return re.sub(substr, start + substr + end, string)


def main():
    np.random.seed(RNG_SEED)
    locations = np.random.randn(NODES_TOTAL, 2)
    locations *= (MEAN_DIST / np.sqrt(2))
    for i in [0] + list(range(-NUM_MUST_RETURN, 0)):
        locations[i] = (0, 0)

    # Set up routing manager/problem
    manager = pywrapcp.RoutingIndexManager(len(locations), NUM_VEHICLES, 0)
    routing = pywrapcp.RoutingModel(manager)

    # Matrixify and convert to integers
    D = squareform(pdist(locations))
    D = (PRECISION * D).round().astype(int)
    DEPT_IDX = [0]
    CUST_IDX = slice(1, NUM_STOPS + 1)
    RETURN_IDX = slice(NUM_STOPS + 1, len(locations))

    # Depot -> RETURN is free
    D[DEPT_IDX, RETURN_IDX] = 0
    # RETURN -> Depot is free
    D[RETURN_IDX, DEPT_IDX] = 0
    # Customer -> Return == Depot Cost
    D[CUST_IDX, RETURN_IDX] = D[CUST_IDX, DEPT_IDX]
    # Return -> Customer is impossible
    D[RETURN_IDX, CUST_IDX] = SOLVER_INFINTIY
    # RETURN -> RETURN is impossible
    D[RETURN_IDX, RETURN_IDX] = SOLVER_INFINTIY
    # Zero the self-routes
    D[np.eye(len(D)).astype(bool)] = 0

    def dist_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return int(D[from_node, to_node])


    def time_callback(from_index, to_index):
        return dist_callback(from_index, to_index) // VEHICLE_SPEED

    DIST_CALLBACK_IDX = routing.RegisterTransitCallback(dist_callback)
    TIME_CALLBACK_IDX = routing.RegisterTransitCallback(time_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(DIST_CALLBACK_IDX)

    MAX_TW = PRECISION * MAX_TIME
    EXTRA_TW = PRECISION * EXTRA_TIME

    VEH_CAP = np.full(NUM_VEHICLES, MAX_TW)
    VEH_CAP += EXTRA_TW

    routing.AddDimensionWithVehicleCapacity(
        TIME_CALLBACK_IDX,
        int(MAX_TW),  # time slack
        VEH_CAP,  # Max travel distance per vehicle
        False,  # don't force start cumul to zero
        "Time",  # Name of the dimension
    )

    # Create customer time windows
    TW = np.random.uniform(size=(NODES_REAL, 2))
    TW = (MAX_TW * TW).round().astype(int)
    TW.sort(axis=1)

    # Force customer windows to be at least MIN_WINDOW big
    TW[:, 0] -= int(MIN_WINDOW * PRECISION) // 2
    TW[:, 1] += int(MIN_WINDOW * PRECISION) // 2
    TW = np.clip(TW, 0, MAX_TW)

    # Apply windows and disjunctions
    # Customers
    time_dimension = routing.GetDimensionOrDie("Time")
    for loc_idx, tw in enumerate(TW):
        if loc_idx == 0:
            continue
        idx = manager.NodeToIndex(loc_idx)
        routing.AddDisjunction([idx], DROP_PENALTY)
        time_dimension.CumulVar(idx).SetRange(int(tw[0]), int(tw[1]))
    
    # Return nodes
    for ret_idx in range(NODES_REAL, NODES_TOTAL):
        idx = manager.NodeToIndex(ret_idx)
        routing.AddDisjunction([idx], DROP_PENALTY)
        time_dimension.CumulVar(idx).SetRange(0, (MAX_TW + EXTRA_TW) // 2)

    # Minimize Global Time
    # time_dimension.SetGlobalSpanCostCoefficient(1)

    #Solve
    search_params = pywrapcp.DefaultRoutingSearchParameters()
    search_params.local_search_metaheuristic = (routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    search_params.log_search = True
    search_params.time_limit.seconds = 30
    tic = time.time()
    solution = routing.SolveWithParameters(search_params)
    toc = time.time()

    # Print Results
    print(
        [
            "ROUTING NOT SOLVED",
            "ROUTING SUCCESS",
            "ROUTING FAIL",
            "ROUTING FAIL TIMEOUT",
            "ROUTING INVALID",
        ][routing.status()]
    )
    print(f"{toc-tic:.1f} sec elapsed")
    stops = []
    for route_nr in range(NUM_VEHICLES):
        index = routing.Start(route_nr)
        node_idx = manager.IndexToNode(index)
        stops.append([node_idx])
        while not routing.IsEnd(index):
            index = solution.Value(routing.NextVar(index))
            node_idx = manager.IndexToNode(index)
            stops[-1].append(node_idx)
    
    dept_idx = [0]
    cust_idx = list(range(1, NODES_REAL))
    ret_idx = list(range(NODES_REAL, NODES_TOTAL))
    drop_idx = []

    for node in range(routing.Size()):
        if routing.IsStart(node) or routing.IsEnd(node):
            continue
        if solution.Value(routing.NextVar(node)) == node:
            drop_idx += [manager.IndexToNode(node)]


    print("Depot:", dept_idx)
    print("Customers:", cust_idx)
    print("Force Return:", Fore.YELLOW, ret_idx, Fore.RESET)
    print("Dropped Nodes:", Fore.RED, (drop_idx), Fore.RESET)
    print("Dropped distance:", D[dept_idx, drop_idx])
    print("Dropped TWs:", TW[drop_idx])

    print("Routes:")
    for i, route in enumerate(stops):
        for idx in ret_idx:
            route = color_substr(str(route), str(idx), "YELLOW")
        print(f"\t{i}: {route}")

if __name__ == "__main__":
    main()