kishorkukreja/awesome-supply-chain/capacitated-vrp

Minimize: Σ_{k∈K} Σ_{i∈V} Σ_{j∈V, i≠j} c_{ij} * x_{ij}^k

1. Each customer visited exactly once:
   Σ_{k∈K} Σ_{i∈V, i≠j} x_{ij}^k = 1,  ∀j ∈ V\{0}
 
2. Flow conservation:
   Σ_{i∈V, i≠h} x_{ih}^k = Σ_{j∈V, j≠h} x_{hj}^k,  ∀h ∈ V, ∀k ∈ K
 
3. Each vehicle leaves depot at most once:
   Σ_{j∈V\{0}} x_{0j}^k ≤ 1,  ∀k ∈ K
 
4. Each vehicle returns to depot at most once:
   Σ_{i∈V\{0}} x_{i0}^k ≤ 1,  ∀k ∈ K
 
5. Capacity constraint (load-based):
   u_j^k ≥ u_i^k + q_j - Q*(1 - x_{ij}^k),  ∀i,j ∈ V\{0}, ∀k ∈ K
   q_i ≤ u_i^k ≤ Q,  ∀i ∈ V\{0}, ∀k ∈ K
 
6. Binary variables:
   x_{ij}^k ∈ {0,1},  ∀i,j ∈ V, ∀k ∈ K

Minimize: Σ_{i∈V} Σ_{j∈V, i≠j} c_{ij} * x_{ij}

1. Degree constraints:
   Σ_{j∈V, j≠i} x_{ij} = 1,  ∀i ∈ V\{0}
   Σ_{i∈V, i≠j} x_{ij} = 1,  ∀j ∈ V\{0}
 
2. Capacity constraints:
   u_j ≥ u_i + q_j - Q*(1 - x_{ij}),  ∀i,j ∈ V\{0}, i≠j
   q_i ≤ u_i ≤ Q,  ∀i ∈ V\{0}
 
3. Valid cuts can be added for strengthening

from pulp import *
import numpy as np
def cvrp_branch_and_cut(dist_matrix, demands, vehicle_capacity,
                        num_vehicles=None, depot=0):
    """
    CVRP using branch-and-cut with capacity inequalities
    Args:
        dist_matrix: n x n distance matrix
        demands: list of customer demands (depot = 0)
        vehicle_capacity: vehicle capacity
        num_vehicles: max number of vehicles (None = unlimited)
        depot: depot index
    Returns:
        solution dictionary
    """
    n = len(dist_matrix)
    customers = [i for i in range(n) if i != depot]
    # Calculate minimum vehicles needed
    total_demand = sum(demands[i] for i in customers)
    min_vehicles = int(np.ceil(total_demand / vehicle_capacity))
    if num_vehicles is None:
        num_vehicles = min_vehicles
    elif num_vehicles < min_vehicles:
        return {
            'status': 'Infeasible',
            'message': f'Need at least {min_vehicles} vehicles for total demand {total_demand}'
        }
    # Create problem
    prob = LpProblem("CVRP", LpMinimize)
    # Decision variables: x[i,j,k] for each vehicle
    x = {}
    for i in range(n):
        for j in range(n):
            if i != j:
                for k in range(num_vehicles):
                    x[i,j,k] = LpVariable(f"x_{i}_{j}_{k}", cat='Binary')
    # Load variables for capacity tracking
    u = {}
    for i in customers:
        for k in range(num_vehicles):
            u[i,k] = LpVariable(f"u_{i}_{k}",
                               lowBound=demands[i],
                               upBound=vehicle_capacity,
                               cat='Continuous')
    # Objective: Minimize total distance
    prob += lpSum([dist_matrix[i][j] * x[i,j,k]
                   for i in range(n) for j in range(n) if i != j
                   for k in range(num_vehicles)]), "Total_Distance"
    # Constraints
    # 1. Each customer visited exactly once
    for j in customers:
        prob += lpSum([x[i,j,k] for i in range(n) if i != j
                      for k in range(num_vehicles)]) == 1, f"Visit_{j}"
    # 2. Flow conservation
    for h in range(n):
        for k in range(num_vehicles):
            prob += (lpSum([x[i,h,k] for i in range(n) if i != h]) ==
                    lpSum([x[h,j,k] for j in range(n) if j != h])), \
                    f"Flow_{h}_{k}"
    # 3. Vehicle capacity (load-based subtour elimination)
    for k in range(num_vehicles):
        for i in customers:
            for j in customers:
                if i != j:
                    prob += (u[i,k] + demands[j] <=
                            u[j,k] + vehicle_capacity * (1 - x[i,j,k])), \
                            f"Load_{i}_{j}_{k}"
    # 4. Each vehicle used at most once
    for k in range(num_vehicles):
        prob += lpSum([x[depot,j,k] for j in customers]) <= 1, \
                f"Leave_Depot_{k}"
    # Solve
    import time
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=1, timeLimit=300))
    solve_time = time.time() - start_time
    # Extract solution
    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        routes = []
        route_loads = []
        for k in range(num_vehicles):
            route = [depot]
            current = depot
            for _ in range(n):
                next_node = None
                for j in range(n):
                    if j != current and (current,j,k) in x:
                        if x[current,j,k].varValue > 0.5:
                            next_node = j
                            break
                if next_node is None or next_node == depot:
                    route.append(depot)
                    break
                route.append(next_node)
                current = next_node
            if len(route) > 2:  # Route has customers
                routes.append(route)
                load = sum(demands[i] for i in route[1:-1])
                route_loads.append(load)
        return {
            'status': LpStatus[prob.status],
            'total_distance': value(prob.objective),
            'routes': routes,
            'route_loads': route_loads,
            'num_vehicles_used': len(routes),
            'solve_time': solve_time
        }
    else:
        return {
            'status': LpStatus[prob.status],
            'solve_time': solve_time
        }
# Example usage
if __name__ == "__main__":
    np.random.seed(42)
    # Generate random problem
    n = 16  # 1 depot + 15 customers
    coords = np.random.rand(n, 2) * 100
    # Distance matrix
    dist_matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            dist_matrix[i][j] = np.linalg.norm(coords[i] - coords[j])
    # Demands (depot has 0)
    demands = [0] + [np.random.randint(5, 25) for _ in range(n-1)]
    vehicle_capacity = 50
    print(f"Total demand: {sum(demands)}")
    print(f"Min vehicles needed: {int(np.ceil(sum(demands) / vehicle_capacity))}")
    result = cvrp_branch_and_cut(dist_matrix, demands, vehicle_capacity)
    if result['status'] in ['Optimal', 'Feasible']:
        print(f"\nStatus: {result['status']}")
        print(f"Total Distance: {result['total_distance']:.2f}")
        print(f"Vehicles Used: {result['num_vehicles_used']}")
        for i, (route, load) in enumerate(zip(result['routes'],
                                              result['route_loads'])):
            print(f"\nVehicle {i+1}: {route}")
            print(f"  Load: {load}/{vehicle_capacity} ({load/vehicle_capacity*100:.1f}%)")
    else:
        print(f"Status: {result['status']}")

def clarke_wright_parallel_cvrp(dist_matrix, demands, vehicle_capacity, depot=0):
    """
    Clarke-Wright Savings Algorithm (Parallel version) for CVRP
    Most famous CVRP heuristic - merges routes based on savings
    Time complexity: O(n^2 log n)
    Solution quality: Typically within 5-10% of optimal
    Args:
        dist_matrix: distance matrix
        demands: customer demands
        vehicle_capacity: vehicle capacity
        depot: depot index
    Returns:
        solution dictionary with routes and costs
    """
    n = len(dist_matrix)
    customers = [i for i in range(n) if i != depot]
    # Calculate savings s_{ij} = d_{0i} + d_{0j} - d_{ij}
    savings = []
    for i in customers:
        for j in customers:
            if i < j:
                saving = (dist_matrix[depot][i] +
                         dist_matrix[depot][j] -
                         dist_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 in separate route
    routes = [[depot, c, depot] for c in customers]
    route_loads = [demands[c] for c in customers]
    # Merge routes based on savings
    for saving_value, i, j in savings:
        # Find routes containing i and j
        route_i_idx = None
        route_j_idx = None
        for idx, route in enumerate(routes):
            if i in route:
                route_i_idx = idx
            if j in route:
                route_j_idx = idx
        if route_i_idx is None or route_j_idx is None:
            continue
        if route_i_idx == route_j_idx:
            continue  # Already in same route
        # Check if i and j are at route ends (interior = 1 and -2 positions)
        route_i = routes[route_i_idx]
        route_j = routes[route_j_idx]
        i_at_end = (route_i[1] == i or route_i[-2] == i)
        j_at_end = (route_j[1] == j or route_j[-2] == j)
        if not (i_at_end and j_at_end):
            continue
        # Check capacity constraint
        combined_load = route_loads[route_i_idx] + route_loads[route_j_idx]
        if combined_load > vehicle_capacity:
            continue
        # Merge routes
        route_i_interior = route_i[1:-1]
        route_j_interior = route_j[1:-1]
        # Determine correct merge order
        if route_i_interior[-1] == i and route_j_interior[0] == j:
            new_route = [depot] + route_i_interior + route_j_interior + [depot]
        elif route_i_interior[-1] == i and route_j_interior[-1] == j:
            new_route = [depot] + route_i_interior + route_j_interior[::-1] + [depot]
        elif route_i_interior[0] == i and route_j_interior[0] == j:
            new_route = [depot] + route_i_interior[::-1] + route_j_interior + [depot]
        elif route_i_interior[0] == i and route_j_interior[-1] == j:
            new_route = [depot] + route_j_interior + route_i_interior + [depot]
        else:
            continue