Vehicle Routing Problem (VRP) integrated with Google Maps Platform.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "198392de",
   "metadata": {},
   "source": [
    "# Vehicle Routing Problem (VRP) integrated with Google Maps Platform"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29823b8a",
   "metadata": {},
   "source": [
    "### Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ac829f39",
   "metadata": {},
   "outputs": [],
   "source": [
    "import googlemaps\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "import folium\n",
    "import random\n",
    "from googlemaps.convert import decode_polyline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecf04ab2",
   "metadata": {},
   "source": [
    "### Create Matrix of Distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "945b00da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_distance_matrix(addresses, api_key):\n",
    "    \"\"\"\n",
    "    Uses the Google Distance Matrix API to build a 2D list of driving\n",
    "    distances (in meters) between all pairs of addresses.\n",
    "    \"\"\"\n",
    "    gmaps = googlemaps.Client(key=api_key)\n",
    "\n",
    "    dm = gmaps.distance_matrix(\n",
    "        origins=addresses,\n",
    "        destinations=addresses,\n",
    "        mode='driving',\n",
    "        units='metric'\n",
    "    )\n",
    "\n",
    "    n = len(addresses)\n",
    "    dist_matrix = [[0.0]*n for _ in range(n)]\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            elem = dm['rows'][i]['elements'][j]\n",
    "            if elem['status'] == 'OK':\n",
    "                dist_matrix[i][j] = float(elem['distance']['value'])\n",
    "            else:\n",
    "                # If no route is found\n",
    "                dist_matrix[i][j] = float('inf')\n",
    "    return dist_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9aa0929",
   "metadata": {},
   "source": [
    "### VRP Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a393bc98",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_vrp(dist_matrix, demands, vehicle_capacity, num_vehicles):\n",
    "    \"\"\"\n",
    "    Solve a basic Capacitated VRP:\n",
    "    - Single depot (node 0)\n",
    "    - N-1 customers\n",
    "    - 'num_vehicles' identical vehicles\n",
    "    - Each vehicle has capacity = vehicle_capacity\n",
    "    - Each customer must be visited exactly once\n",
    "\n",
    "    Returns (routes, obj_val):\n",
    "      routes: list of routes, routes[k] = the sequence of visited nodes for vehicle k\n",
    "      obj_val: total distance of the solution\n",
    "    \"\"\"\n",
    "    n = len(dist_matrix)\n",
    "    depot = 0\n",
    "    customers = range(1, n)\n",
    "    K = num_vehicles\n",
    "\n",
    "    m = gp.Model(\"CVRP\")\n",
    "\n",
    "    # x[k,i,j] = 1 if vehicle k travels from i to j\n",
    "    x = m.addVars(K, n, n, vtype=GRB.BINARY, name=\"x\")\n",
    "    # u[k,i] = load on vehicle k upon arriving at node i\n",
    "    u = m.addVars(K, n, lb=0, ub=vehicle_capacity, vtype=GRB.CONTINUOUS, name=\"u\")\n",
    "\n",
    "    # Objective: minimize total distance\n",
    "    m.setObjective(\n",
    "        gp.quicksum(dist_matrix[i][j] * x[k,i,j] for k in range(K) for i in range(n) for j in range(n)),\n",
    "        GRB.MINIMIZE\n",
    "    )\n",
    "\n",
    "    # 1) Each customer is visited exactly once (sum over all vehicles, and all outgoing arcs):\n",
    "    for i in customers:\n",
    "        # Ensure that each customer i is visited exactly once by summing over all vehicles and possible departures.\n",
    "        m.addConstr(\n",
    "            gp.quicksum(x[k,i,j] for k in range(K) for j in range(n)) == 1,\n",
    "            name=f\"visit_{i}\"\n",
    "        )\n",
    "\n",
    "    for k in range(K):\n",
    "        # Ensure that each vehicle k leaves the depot exactly once (excluding self-loop from depot).\n",
    "        m.addConstr(\n",
    "            gp.quicksum(x[k,depot,j] for j in range(n) if j != 0) == 1,\n",
    "            name=f\"visit_{i}\"\n",
    "        )\n",
    "\n",
    "    for k in range(K):\n",
    "        # Ensure that each vehicle k returns to the depot exactly once (excluding self-loop to depot).\n",
    "        m.addConstr(\n",
    "            gp.quicksum(x[k,i,depot] for i in range(n) if i != 0) == 1,\n",
    "            name=f\"visit_{i}\"\n",
    "        )\n",
    "\n",
    "    # 2) Flow conservation for each vehicle and each node (in = out):\n",
    "    for k in range(K):\n",
    "        for i in range(n):\n",
    "            # Ensure the number of arcs entering node i equals the number leaving node i for vehicle k.\n",
    "            m.addConstr(\n",
    "                gp.quicksum(x[k,j,i] for j in range(n)) == gp.quicksum(x[k,i,j] for j in range(n)),\n",
    "                name=f\"flow_{k}_{i}\"\n",
    "            )\n",
    "\n",
    "    # 3) Depot in/out must match for each vehicle if it’s used\n",
    "    #    (the number of times a vehicle leaves the depot equals the number of times it returns).\n",
    "    for k in range(K):\n",
    "        # Ensure that for vehicle k, departures from the depot (to customers) equal returns to the depot.\n",
    "        m.addConstr(\n",
    "            gp.quicksum(x[k,depot,j] for j in customers) == gp.quicksum(x[k,i,depot] for i in customers),\n",
    "            name=f\"depot_inout_{k}\"\n",
    "        )\n",
    "\n",
    "    # 4) Capacity constraints (MTZ):\n",
    "    for k in range(K):\n",
    "        # Start each vehicle k with load 0 at the depot.\n",
    "        m.addConstr(u[k,depot] == 0, name=f\"load_depot_{k}\")\n",
    "\n",
    "        for i in customers:\n",
    "            # Check if the customer demand exceeds the vehicle capacity (feasibility check).\n",
    "            if demands[i] > vehicle_capacity:\n",
    "                print(f\"Customer {i} demand = {demands[i]} exceeds vehicle capacity {vehicle_capacity}.\")\n",
    "                # This will likely be infeasible\n",
    "            for j in customers:\n",
    "                if i != j:\n",
    "                    # MTZ constraint: ensure that the load carried is adjusted correctly when traveling from i to j.\n",
    "                    m.addConstr(\n",
    "                        u[k,i] - u[k,j] + vehicle_capacity * x[k,i,j] <= vehicle_capacity - demands[j],\n",
    "                        name=f\"cap_{k}_{i}_{j}\"\n",
    "                    )\n",
    "\n",
    "    # No self-loops:\n",
    "    for k in range(K):\n",
    "        for i in range(n):\n",
    "            # Prevent vehicle k from traveling from node i directly back to itself.\n",
    "            m.addConstr(x[k,i,i] == 0, name=f\"no_loop_{k}_{i}\")\n",
    "\n",
    "    for i in customers:\n",
    "        for k in range(K):\n",
    "            # Ensure that when a vehicle k visits customer i, the load is at least the demand of i.\n",
    "            m.addConstr(u[k,i] >= demands[i])\n",
    "                        \n",
    "    m.optimize()\n",
    "\n",
    "    if m.status == GRB.OPTIMAL:\n",
    "        obj_val = m.objVal\n",
    "\n",
    "        print(\"\\nArc decisions x[k,i,j] in the optimal solution:\")\n",
    "        for k in range(K):\n",
    "            for i in range(n):\n",
    "                for j in range(n):\n",
    "                    if x[k,i,j].X > 0.5:\n",
    "                        print(f\"x[{k}][{i}][{j}] = 1\")\n",
    "\n",
    "        # Reconstruct each vehicle’s route\n",
    "        routes = []\n",
    "        for k in range(K):\n",
    "            route_k = [depot]\n",
    "            current = depot\n",
    "            while True:\n",
    "                next_j = None\n",
    "                # find j where x[k, current, j] == 1\n",
    "                for j in range(n):\n",
    "                    if x[k,current,j].X > 0.5:\n",
    "                        next_j = j\n",
    "                        break\n",
    "                if next_j is None:\n",
    "                    # no next move => route ended (vehicle not used or route finished)\n",
    "                    break\n",
    "                if next_j == depot:\n",
    "                    # returned to depot => route closed\n",
    "                    route_k.append(depot)\n",
    "                    break\n",
    "                else:\n",
    "                    route_k.append(next_j)\n",
    "                    current = next_j\n",
    "\n",
    "            routes.append(route_k)\n",
    "\n",
    "        return routes, obj_val\n",
    "    else:\n",
    "        print(\"No optimal solution found or model infeasible.\")\n",
    "        return None, None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5005c33",
   "metadata": {},
   "source": [
    "### Latitude/longitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "82f448ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "def geocode_addresses(addresses, api_key):\n",
    "    \"\"\"\n",
    "    Get latitude/longitude for each address using Google Geocoding API.\n",
    "    Returns list of (lat, lng) in the same order as 'addresses'.\n",
    "    \"\"\"\n",
    "    gmaps = googlemaps.Client(key=api_key)\n",
    "    coords = []\n",
    "    for addr in addresses:\n",
    "        result = gmaps.geocode(addr)\n",
    "        if result:\n",
    "            loc = result[0]['geometry']['location']\n",
    "            coords.append((loc['lat'], loc['lng']))\n",
    "        else:\n",
    "            coords.append((0.0, 0.0))\n",
    "    return coords"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4375f8fc",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "065c5f51",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_vrp_routes_on_map(addresses, coords, routes, api_key, outfile=\"vrp_route.html\"):\n",
    "    \"\"\"\n",
    "    Plots each vehicle's route with a different color on\n",
    "    a Folium map using the road paths.\n",
    "    \"\"\"\n",
    "    gmaps = googlemaps.Client(key=api_key)\n",
    "\n",
    "    # Initialize map around the depot (node 0)\n",
    "    start_lat, start_lng = coords[0]\n",
    "    fmap = folium.Map(location=[start_lat, start_lng], zoom_start=13)\n",
    "\n",
    "    # Some distinct route colors. Expand if you have many vehicles.\n",
    "    route_colors = [\n",
    "        \"red\", \"blue\", \"green\", \"purple\", \"orange\", \"darkred\", \n",
    "        \"cadetblue\", \"darkblue\", \"darkgreen\", \"pink\"\n",
    "    ]\n",
    "    # Optional\n",
    "    dash_styles = [None, \"5,5\", \"1,5\", \"5,10\"]\n",
    "\n",
    "    for k, route in enumerate(routes):\n",
    "        color = route_colors[k % len(route_colors)]\n",
    "        dash = dash_styles[k % len(dash_styles)] if dash_styles else None\n",
    "\n",
    "        # Place markers for each stop\n",
    "        for order, node_idx in enumerate(route):\n",
    "            lat, lng = coords[node_idx]\n",
    "            popup_text = f\"Vehicle {k+1} - Stop {order+1}: {addresses[node_idx]}\"\n",
    "            # depot marker is different\n",
    "            if node_idx == 0:\n",
    "                folium.Marker([lat, lng], popup=popup_text,\n",
    "                              icon=folium.Icon(color='blue', icon='info-sign')\n",
    "                             ).add_to(fmap)\n",
    "            else:\n",
    "                folium.Marker([lat, lng], popup=popup_text,\n",
    "                              icon=folium.Icon(color='green', icon='info-sign')\n",
    "                             ).add_to(fmap)\n",
    "\n",
    "        # Draw polylines for each pair in the route\n",
    "        for i in range(len(route) - 1):\n",
    "            origin_idx = route[i]\n",
    "            dest_idx = route[i+1]\n",
    "            origin_lat, origin_lng = coords[origin_idx]\n",
    "            dest_lat, dest_lng = coords[dest_idx]\n",
    "\n",
    "            directions_result = gmaps.directions(\n",
    "                (origin_lat, origin_lng),\n",
    "                (dest_lat, dest_lng),\n",
    "                mode=\"driving\"\n",
    "            )\n",
    "            if directions_result:\n",
    "                polyline_str = directions_result[0]['overview_polyline']['points']\n",
    "                polyline_points = decode_polyline(polyline_str)\n",
    "                folium.PolyLine(\n",
    "                    locations=[(pt['lat'], pt['lng']) for pt in polyline_points],\n",
    "                    color=color,\n",
    "                    weight=4,\n",
    "                    opacity=1,\n",
    "                    dash_array=dash  \n",
    "                ).add_to(fmap)\n",
    "\n",
    "    fmap.save(outfile)\n",
    "    print(f\"Map saved to {outfile}. Open 'vrp_route.html'.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d15f16c",
   "metadata": {},
   "source": [
    "### Choose Locations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f88e00f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 8 locations around Midtown Atlanta, with starting point at Georgia Tech\n",
    "addresses = [\n",
    "    \"Georgia Tech, 225 North Ave NW, Atlanta, GA 30332\",\n",
    "    \"Fox Theatre, 660 Peachtree St NE, Atlanta, GA 30308\",\n",
    "    \"Atlanta Botanical Garden, 1345 Piedmont Ave NE, Atlanta, GA 30309\",\n",
    "    \"Ponce City Market, 675 Ponce De Leon Ave NE, Atlanta, GA 30308\",\n",
    "    \"Coca-Cola Headquarters, 1 Coca Cola Plz NW, Atlanta, GA 30313\",\n",
    "    \"Mercedes-Benz Stadium, 1 AMB Dr NW, Atlanta, GA 30313\",\n",
    "    \"Georgia Aquarium, 225 Baker St NW, Atlanta, GA 30313\",\n",
    "    \"World of Coca-Cola, 121 Baker St NW, Atlanta, GA 30313\"\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af98dc6c",
   "metadata": {},
   "source": [
    "### API Key"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6a27fe0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Replace with your valid Google Cloud API key\n",
    "api_key = \"YOUR_API\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0de5e5ae",
   "metadata": {},
   "source": [
    "### Solve the Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0b0ad6a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Demands: [0, 3, 1, 1, 3, 2, 1, 1]\n",
      "Set parameter Username\n",
      "Academic license - for non-commercial use only - expires 2026-03-13\n",
      "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (win64 - Windows 11+.0 (26100.2))\n",
      "\n",
      "CPU model: AMD Ryzen 7 4800H with Radeon Graphics, instruction set [SSE2|AVX|AVX2]\n",
      "Thread count: 8 physical cores, 16 logical processors, using up to 16 threads\n",
      "\n",
      "Optimize a model with 214 rows, 216 columns and 1014 nonzeros\n",
      "Model fingerprint: 0xe357e0b5\n",
      "Variable types: 24 continuous, 192 integer (192 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 4e+00]\n",
      "  Objective range  [2e+02, 6e+03]\n",
      "  Bounds range     [1e+00, 4e+00]\n",
      "  RHS range        [1e+00, 3e+00]\n",
      "Presolve removed 72 rows and 45 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 142 rows, 171 columns, 1221 nonzeros\n",
      "Variable types: 21 continuous, 150 integer (150 binary)\n",
      "\n",
      "Root relaxation: objective 1.476900e+04, 72 iterations, 0.00 seconds (0.00 work units)\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 14769.0000    0    6          - 14769.0000      -     -    0s\n",
      "H    0     0                    20669.000000 14769.0000  28.5%     -    0s\n",
      "     0     0 14769.0000    0    6 20669.0000 14769.0000  28.5%     -    0s\n",
      "     0     0 14826.0000    0   15 20669.0000 14826.0000  28.3%     -    0s\n",
      "H    0     0                    17687.000000 14854.5000  16.0%     -    0s\n",
      "     0     0 14854.5000    0   15 17687.0000 14854.5000  16.0%     -    0s\n",
      "     0     0 14854.5000    0   20 17687.0000 14854.5000  16.0%     -    0s\n",
      "H    0     0                    17653.000000 14879.0000  15.7%     -    0s\n",
      "H    0     0                    17602.000000 14879.0000  15.5%     -    0s\n",
      "     0     0 14879.0000    0    6 17602.0000 14879.0000  15.5%     -    0s\n",
      "     0     0 14989.0000    0    6 17602.0000 14989.0000  14.8%     -    0s\n",
      "     0     0 14989.0000    0    6 17602.0000 14989.0000  14.8%     -    0s\n",
      "H    0     0                    17579.000000 14989.0000  14.7%     -    0s\n",
      "     0     0 14989.0000    0    6 17579.0000 14989.0000  14.7%     -    0s\n",
      "     0     0 14989.0000    0    6 17579.0000 14989.0000  14.7%     -    0s\n",
      "     0     0 15020.6667    0   21 17579.0000 15020.6667  14.6%     -    0s\n",
      "H    0     0                    17366.000000 15077.0147  13.2%     -    0s\n",
      "     0     0 15154.7500    0   21 17366.0000 15154.7500  12.7%     -    0s\n",
      "     0     0 15154.7500    0   22 17366.0000 15154.7500  12.7%     -    0s\n",
      "H    0     0                    17325.000000 15154.7500  12.5%     -    0s\n",
      "     0     0 15154.7500    0   18 17325.0000 15154.7500  12.5%     -    0s\n",
      "     0     0 15187.0000    0   22 17325.0000 15187.0000  12.3%     -    0s\n",
      "     0     0 15187.0000    0   22 17325.0000 15187.0000  12.3%     -    0s\n",
      "     0     2 15372.0000    0   22 17325.0000 15372.0000  11.3%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Implied bound: 10\n",
      "  Clique: 12\n",
      "  MIR: 18\n",
      "  Zero half: 1\n",
      "  Mod-K: 6\n",
      "  RLT: 13\n",
      "\n",
      "Explored 92 nodes (1205 simplex iterations) in 0.12 seconds (0.05 work units)\n",
      "Thread count was 16 (of 16 available processors)\n",
      "\n",
      "Solution count 7: 17325 17366 17579 ... 20669\n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 1.732500000000e+04, best bound 1.732500000000e+04, gap 0.0000%\n",
      "\n",
      "Arc decisions x[k,i,j] in the optimal solution:\n",
      "x[0][0][2] = 1\n",
      "x[0][2][4] = 1\n",
      "x[0][4][0] = 1\n",
      "x[1][0][1] = 1\n",
      "x[1][1][3] = 1\n",
      "x[1][3][0] = 1\n",
      "x[2][0][6] = 1\n",
      "x[2][5][7] = 1\n",
      "x[2][6][5] = 1\n",
      "x[2][7][0] = 1\n",
      "\n",
      "Solution found with total distance = 17325.00 meters\n",
      "Vehicle 1: [0, 2, 4, 0] (demands covered: [1, 3])\n",
      "Vehicle 2: [0, 1, 3, 0] (demands covered: [3, 1])\n",
      "Vehicle 3: [0, 6, 5, 7, 0] (demands covered: [1, 2, 1])\n",
      "Map saved to vrp_route.html. Open 'vrp_route.html'.\n"
     ]
    }
   ],
   "source": [
    "# 1. Build distance matrix\n",
    "dist_matrix = get_distance_matrix(addresses, api_key)\n",
    "n = len(dist_matrix)\n",
    "\n",
    "# 2. Create demand\n",
    "random.seed(42)\n",
    "demands = [0] + [random.randint(1,3) for _ in range(n-1)]\n",
    "print(\"Demands:\", demands)\n",
    "\n",
    "# 3. Solve VRP\n",
    "vehicle_capacity = 4\n",
    "num_vehicles = 3\n",
    "routes, objVal = solve_vrp(dist_matrix, demands, vehicle_capacity, num_vehicles)\n",
    "\n",
    "print(f\"\\nSolution found with total distance = {objVal:.2f} meters\")\n",
    "for k, r in enumerate(routes):\n",
    "    print(f\"Vehicle {k+1}: {r} (demands covered: {[demands[node] for node in r if node != 0]})\")\n",
    "\n",
    "# 4. Plot the route with real roads\n",
    "coords = geocode_addresses(addresses, api_key)\n",
    "\n",
    "# 5) Plot each vehicle’s route in a different color/style\n",
    "plot_vrp_routes_on_map(addresses, coords, routes, api_key, outfile=\"vrp_route.html\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}