kishorkukreja/awesome-supply-chain/route-optimization

from pulp import *
import numpy as np
def vrp_mip_model(distances, demands, vehicle_capacity, num_vehicles):
    """
    Basic VRP using Mixed-Integer Programming
    Parameters:
    - distances: NxN distance matrix
    - demands: list of demands at each location
    - vehicle_capacity: max capacity per vehicle
    - num_vehicles: number of available vehicles
    """
    n = len(distances)  # includes depot at index 0
    customers = range(1, n)
    all_nodes = range(n)
    vehicles = range(num_vehicles)
    # Create problem
    prob = LpProblem("VRP", LpMinimize)
    # Decision variables
    # x[i,j,k] = 1 if vehicle k travels from i to j
    x = LpVariable.dicts("route",
                         [(i,j,k) for i in all_nodes
                          for j in all_nodes
                          for k in vehicles if i != j],
                         cat='Binary')
    # u[i,k] = position of node i in route of vehicle k (for subtour elimination)
    u = LpVariable.dicts("order",
                         [(i,k) for i in customers for k in vehicles],
                         lowBound=1,
                         upBound=n-1,
                         cat='Integer')
    # Objective: minimize total distance
    prob += lpSum([distances[i][j] * x[i,j,k]
                   for i in all_nodes
                   for j in all_nodes
                   for k in vehicles
                   if i != j])
    # Constraints
    # 1. Each customer visited exactly once
    for j in customers:
        prob += lpSum([x[i,j,k]
                      for i in all_nodes
                      for k in vehicles
                      if i != j]) == 1
    # 2. Flow conservation: enter and leave each node
    for k in vehicles:
        for j in all_nodes:
            prob += (lpSum([x[i,j,k] for i in all_nodes if i != j]) ==
                    lpSum([x[j,i,k] for i in all_nodes if i != j]))
    # 3. Each vehicle starts from depot
    for k in vehicles:
        prob += lpSum([x[0,j,k] for j in customers]) <= 1
    # 4. Each vehicle returns to depot
    for k in vehicles:
        prob += lpSum([x[i,0,k] for i in customers]) <= 1
    # 5. Vehicle capacity constraint
    for k in vehicles:
        prob += lpSum([demands[j] * x[i,j,k]
                      for i in all_nodes
                      for j in customers
                      if i != j]) <= vehicle_capacity
    # 6. Subtour elimination (Miller-Tucker-Zemlin)
    for i in customers:
        for j in customers:
            for k in vehicles:
                if i != j:
                    prob += u[i,k] - u[j,k] + n * x[i,j,k] <= n - 1
    # Solve
    prob.solve(PULP_CBC_CMD(msg=1))
    # Extract routes
    routes = extract_routes(x, vehicles, n)
    return {
        'status': LpStatus[prob.status],
        'total_distance': value(prob.objective),
        'routes': routes
    }
def extract_routes(x, vehicles, n):
    """Extract route sequences from solution"""
    routes = []
    for k in vehicles:
        route = [0]  # start at depot
        current = 0
        while True:
            for j in range(n):
                if j != current and (current,j,k) in x:
                    if x[current,j,k].varValue > 0.5:
                        if j == 0:  # back to depot
                            route.append(0)
                            if len(route) > 2:  # non-empty route
                                routes.append(route)
                            current = -1
                            break
                        else:
                            route.append(j)
                            current = j
                            break
            if current == -1 or current == 0:
                break
    return routes

import numpy as np
def nearest_neighbor_vrp(distance_matrix, demands, capacity, depot=0):
    """
    Nearest neighbor heuristic for VRP
    Simple construction heuristic - builds routes by always
    selecting the nearest unvisited customer
    """
    n = len(distance_matrix)
    unvisited = set(range(1, n))  # exclude depot
    routes = []
    while unvisited:
        route = [depot]
        current_load = 0
        current_location = depot
        while unvisited:
            # Find nearest feasible customer
            nearest = None
            nearest_dist = float('inf')
            for customer in unvisited:
                if current_load + demands[customer] <= capacity:
                    dist = distance_matrix[current_location][customer]
                    if dist < nearest_dist:
                        nearest = customer
                        nearest_dist = dist
            if nearest is None:
                break  # No feasible customer, return to depot
            route.append(nearest)
            current_load += demands[nearest]
            current_location = nearest
            unvisited.remove(nearest)
        route.append(depot)  # return to depot
        routes.append(route)
    total_distance = sum(calculate_route_distance(route, distance_matrix)
                        for route in routes)
    return {
        'routes': routes,
        'total_distance': total_distance,
        'num_vehicles': len(routes)
    }
def calculate_route_distance(route, distance_matrix):
    """Calculate total distance of a route"""
    distance = 0
    for i in range(len(route) - 1):
        distance += distance_matrix[route[i]][route[i+1]]
    return distance

def clarke_wright_savings(distance_matrix, demands, capacity):
    """
    Clarke-Wright Savings Algorithm for VRP
    Starts with individual routes to each customer,
    then merges routes based on savings
    """
    n = len(distance_matrix)
    depot = 0
    # Calculate savings for all customer pairs
    savings = []
    for i in range(1, n):
        for j in range(i+1, n):
            saving = (distance_matrix[depot][i] +
                     distance_matrix[depot][j] -
                     distance_matrix[i][j])
            savings.append((saving, i, j))
    # Sort savings in descending order
    savings.sort(reverse=True, key=lambda x: x[0])
    # Initialize: each customer has its own route
    routes = {i: [depot, i, depot] for i in range(1, n)}
    route_loads = {i: demands[i] for i in range(1, n)}
    # Try to merge routes based on savings
    for saving_value, i, j in savings:
        # Find routes containing i and j
        route_i = None
        route_j = None
        for route_id, route in routes.items():
            if i in route:
                route_i = route_id
            if j in route:
                route_j = route_id
        # Can only merge if i and j are in different routes
        if route_i != route_j:
            route_i_obj = routes[route_i]
            route_j_obj = routes[route_j]
            # Check if merging is feasible
            combined_load = route_loads[route_i] + route_loads[route_j]
            if combined_load <= capacity:
                # Check if i and j are at ends of their routes
                i_at_end = (route_i_obj[1] == i or route_i_obj[-2] == i)
                j_at_end = (route_j_obj[1] == j or route_j_obj[-2] == j)
                if i_at_end and j_at_end:
                    # Merge routes
                    new_route = merge_routes(route_i_obj, route_j_obj, i, j)
                    # Update routes
                    routes[route_i] = new_route
                    route_loads[route_i] = combined_load
                    del routes[route_j]
                    del route_loads[route_j]
    return {
        'routes': list(routes.values()),
        'total_distance': sum(calculate_route_distance(route, distance_matrix)
                            for route in routes.values()),
        'num_vehicles': len(routes)
    }
def merge_routes(route1, route2, i, j):
    """Merge two routes through nodes i and j"""
    # Remove depot endpoints
    r1 = route1[1:-1]
    r2 = route2[1:-1]
    # Orient routes so i and j connect properly
    if r1[-1] == i and r2[0] == j:
        merged = [0] + r1 + r2 + [0]
    elif r1[-1] == i and r2[-1] == j:
        merged = [0] + r1 + r2[::-1] + [0]
    elif r1[0] == i and r2[0] == j:
        merged = [0] + r1[::-1] + r2 + [0]
    elif r1[0] == i and r2[-1] == j:
        merged = [0] + r2 + r1 + [0]
    else:
        merged = [0] + r1 + r2 + [0]  # fallback
    return merged

import math
def sweep_algorithm(coordinates, demands, capacity, depot_idx=0):
    """
    Sweep Algorithm for VRP
    Sorts customers by polar angle from depot,
    then builds routes by sweeping clockwise
    """
    depot = coordinates[depot_idx]
    n = len(coordinates)