{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiple Traveling Salesman Problem\n",
    "\n",
    "## 1. Mathematical Problem\n",
    "\n",
    "In this notebook, we show how to solve the Multiple Traveling Salesman Problem (mTSP) using CVXPY. The problem considers $m$ traveling salesmen. To solve it, I'm going to use the Miller-Tucker-Zemlin formulation, which follows:\n",
    "\n",
    "The cities are identified with the numbers $1, \\ldots, n$, with which we define:\n",
    "\n",
    "$$\n",
    "x_{ij} = \\begin{cases} 1 & \\text{the path goes from the city} \\, i \\text{ to city } j \\\\ 0 & \\text{otherwise} \\end{cases}\n",
    "$$\n",
    "\n",
    "For $i = 1, \\ldots, n$, $m$ is the number of salesmen, $u_{i}$ is an auxiliary variable and $c_{i, j}$ is the distance from city $i$ to the city $j$. Then the mTSP can be written as the following integer linear programming problem:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\min &\\sum_{i=1}^n \\sum_{j\\ne i,j=1}^n c_{i,j}x_{i,j}\\colon &&  \\\\\n",
    "     & \\sum_{i=2}^n x_{i,1} = m \\\\\n",
    "     & \\sum_{j=2}^n x_{1,j} = m \\\\\n",
    "     & \\sum_{i=2,i\\ne j}^n x_{i,j} = 1 && j=1, \\ldots, n; \\\\\n",
    "     & \\sum_{j=2,j\\ne i}^n x_{i,j} = 1 && i=1, \\ldots, n; \\\\\n",
    "     & u_{i}-u_{j} + 1 \\le (n-1)(1-x_{i,j}) && 2 \\le i \\ne j \\le n;  \\\\\n",
    "     & 0 \\le u_{i} \\le n && 2 \\le i \\le n; \\\\\n",
    "     & x_{i,j} \\in \\{0,1\\}  && i,j=1, \\ldots, n; \\\\\n",
    "     & u_{i} \\in \\mathbf{Z} && i=2, \\ldots, n. \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The objective function minimizes the distance of each salesman's routes. The first two restrictions guarantee that there is a number of salesmen $m$ who leave and return to the city of origin. The following two restrictions ensure that each city is reached from exactly one other city and that from each city there is an exit to exactly one other city. The constraints of the auxiliary variables force that no salesman passes through the same city twice.\n",
    "\n",
    "## 2. Example for a single Traveling Salesman\n",
    "\n",
    "We are going to find the optimal paths (shortest distance) for a single salesman that starts from the point of origin (O) and must pass through all points A to F and return to point O.\n",
    "\n",
    "| Point |\tCoordinates |\n",
    "| :----------: | :----------: |\n",
    "| Origin(O) | (-12.059296, -76.975893) |\n",
    "| A | (-12.079575, -77.009686) |\n",
    "| B | (-12.087303, -76.996620) |\n",
    "| C | (-12.084391, -76.975651) |\n",
    "| D | (-12.063603, -76.960483) |\n",
    "| E | (-12.056762, -77.014452) |\n",
    "| F | (-12.011531, -77.002383) |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "code_folding": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distance Matrix is:\n",
      "\n",
      "[[0.     4.3092 3.8329 2.7762 1.7441 4.2076 6.0199]\n",
      " [4.3092 0.     1.6596 3.7435 5.6408 2.5764 7.5691]\n",
      " [3.8329 1.6596 0.     2.3055 4.7278 3.8966 8.4056]\n",
      " [2.7762 3.7435 2.3055 0.     2.8311 5.2142 8.5694]\n",
      " [1.7441 5.6408 4.7278 2.8311 0.     5.9246 7.3483]\n",
      " [4.2076 2.5764 3.8966 5.2142 5.9246 0.     5.1733]\n",
      " [6.0199 7.5691 8.4056 8.5694 7.3483 5.1733 0.    ]]\n"
     ]
    }
   ],
   "source": [
    "################################################\n",
    "# Loading libraries\n",
    "################################################\n",
    "!pip install --quiet geopy\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import cvxpy as cp\n",
    "from geopy import distance # Library for geographical calculations\n",
    "\n",
    "################################################\n",
    "# Original Data\n",
    "################################################\n",
    "points = [(-12.059296, -76.975893),\n",
    "          (-12.079575, -77.009686),\n",
    "          (-12.087303, -76.996620),\n",
    "          (-12.084391, -76.975651),\n",
    "          (-12.063603, -76.960483),\n",
    "          (-12.056762, -77.014452),\n",
    "          (-12.011531, -77.002383)]\n",
    "\n",
    "################################################\n",
    "# Building distance matrix\n",
    "################################################\n",
    "n = len(points)\n",
    "C = np.zeros((n,n))\n",
    "\n",
    "for i in range(0, n):\n",
    "    for j in range(0, len(points)):\n",
    "        C[i,j] = distance.distance(points[i], points[j]).km\n",
    "\n",
    "# Showing distance matrix\n",
    "print('Distance Matrix is:\\n')\n",
    "print(np.round(C,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The path is:\n",
      "\n",
      "0 => 6 => 5 => 1 => 2 => 3 => 4 => 0\n"
     ]
    }
   ],
   "source": [
    "################################################\n",
    "# Solving the integer programming problem\n",
    "################################################\n",
    "\n",
    "# Defining variables\n",
    "X = cp.Variable(C.shape, boolean=True)\n",
    "u = cp.Variable(n, integer=True)\n",
    "ones = np.ones((n,1))\n",
    "\n",
    "# Defining the objective function\n",
    "objective = cp.Minimize(cp.sum(cp.multiply(C, X)))\n",
    "\n",
    "# Defining the constraints\n",
    "constraints = []\n",
    "constraints += [X @ ones == ones]\n",
    "constraints += [X.T @ ones == ones]\n",
    "constraints += [cp.diag(X) == 0]\n",
    "constraints += [u[1:] >= 2]\n",
    "constraints += [u[1:] <= n]\n",
    "constraints += [u[0] == 1]\n",
    "\n",
    "for i in range(1, n):\n",
    "    for j in range(1, n):\n",
    "        if i != j:\n",
    "            constraints += [ u[i] - u[j] + 1  <= (n - 1) * (1 - X[i, j]) ]\n",
    "\n",
    "# Solving the problem\n",
    "prob = cp.Problem(objective, constraints)\n",
    "prob.solve(verbose=False)\n",
    "\n",
    "# Transforming the solution to a path\n",
    "X_sol = np.argwhere(X.value==1)\n",
    "orden = X_sol[0].tolist()\n",
    "\n",
    "for i in range(1, n):\n",
    "    row = orden[-1]\n",
    "    orden.append(X_sol[row,1])\n",
    "\n",
    "# Showing the optimal path\n",
    "print('The path is:\\n')\n",
    "print( ' => '.join(map(str, orden)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc603f89730>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 36744 characters]\n",
      "text/plain": [
       "<Figure size 1008x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "################################################\n",
    "# Plotting the optimal path\n",
    "################################################\n",
    "\n",
    "# Transforming the points to the xy plane approximately\n",
    "xy_cords = np.zeros((n,2))\n",
    "\n",
    "for i in range(0, n):\n",
    "    xy_cords[i,0] = distance.distance((points[0][1],0), (points[i][1],0)).km\n",
    "    xy_cords[i,1] = distance.distance((0,points[0][0]), (0,points[i][0])).km\n",
    "\n",
    "# Plotting the points\n",
    "fig, ax = plt.subplots(figsize=(14,7))\n",
    "\n",
    "for i in range(n):\n",
    "    ax.annotate(str(i), xy=(xy_cords[i,0], xy_cords[i,1]+0.1))\n",
    "    \n",
    "ax.scatter(xy_cords[:,0],xy_cords[:,1])\n",
    "ax.plot(xy_cords[orden,0], xy_cords[orden,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The optimal distance is: 22.31 km\n"
     ]
    }
   ],
   "source": [
    "# Showing the optimal distance\n",
    "distance = np.sum(np.multiply(C, X.value))\n",
    "print('The optimal distance is:', np.round(distance,2), 'km')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Example for three Traveling Salesmen\n",
    "\n",
    "Continuing with the previous example, now we are going to find the optimal paths (shortest distance) for three salesmen that start from the point of origin (O) and between all of them must go through points A to F and return to point O.\n",
    "\n",
    "| Point |\tCoordinates |\n",
    "| :----------: | :----------: |\n",
    "| Origin(O) | (-12.059296, -76.975893) |\n",
    "| A | (-12.079575, -77.009686) |\n",
    "| B | (-12.087303, -76.996620) |\n",
    "| C | (-12.084391, -76.975651) |\n",
    "| D | (-12.063603, -76.960483) |\n",
    "| E | (-12.056762, -77.014452) |\n",
    "| F | (-12.011531, -77.002383) |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The path of Salesman_1 is:\n",
      "\n",
      "0 => 3 => 0\n",
      "\n",
      "The path of Salesman_2 is:\n",
      "\n",
      "0 => 4 => 0\n",
      "\n",
      "The path of Salesman_3 is:\n",
      "\n",
      "0 => 6 => 5 => 1 => 2 => 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "################################################\n",
    "# Solving the integer programming problem\n",
    "################################################\n",
    "\n",
    "# Defining the variables\n",
    "X = cp.Variable(C.shape, boolean=True)\n",
    "u = cp.Variable(n, integer=True)\n",
    "m = 3\n",
    "ones = np.ones((n,1))\n",
    "\n",
    "# Defining the objective function\n",
    "objective = cp.Minimize(cp.sum(cp.multiply(C, X)))\n",
    "\n",
    "# Defining the constraints\n",
    "constraints = []\n",
    "constraints += [X[0,:] @ ones == m]\n",
    "constraints += [X[:,0] @ ones == m]\n",
    "constraints += [X[1:,:] @ ones == 1]\n",
    "constraints += [X[:,1:].T @ ones == 1]\n",
    "constraints += [cp.diag(X) == 0]\n",
    "constraints += [u[1:] >= 2]\n",
    "constraints += [u[1:] <= n]\n",
    "constraints += [u[0] == 1]\n",
    "\n",
    "for i in range(1, n):\n",
    "    for j in range(1, n):\n",
    "        if i != j:\n",
    "            constraints += [ u[i] - u[j] + 1  <= (n - 1) * (1 - X[i, j]) ]\n",
    "\n",
    "# Solving the problem\n",
    "prob = cp.Problem(objective, constraints)\n",
    "prob.solve(verbose=False)\n",
    "\n",
    "# Transforming the solution to paths\n",
    "X_sol = np.argwhere(X.value==1)\n",
    "\n",
    "ruta = {}\n",
    "for i in range(0, m):\n",
    "    ruta['Salesman_' + str(i+1)] = [0]\n",
    "    j = i\n",
    "    a = 10e10\n",
    "    while a != 0:\n",
    "        a = X_sol[j,1]\n",
    "        ruta['Salesman_' + str(i+1)].append(a)\n",
    "        j = np.where(X_sol[:,0] == a)\n",
    "        j = j[0][0]\n",
    "        a = j\n",
    "\n",
    "# Showing the paths\n",
    "for i in ruta.keys():\n",
    "    print('The path of ' + i + ' is:\\n')\n",
    "    print( ' => '.join(map(str, ruta[i])))\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 43256 characters]\n",
      "text/plain": [
       "<Figure size 1008x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "################################################\n",
    "# Plotting the optimal path\n",
    "################################################\n",
    "\n",
    "# Plotting the points\n",
    "fig, ax = plt.subplots(figsize=(14,7))\n",
    "\n",
    "for i in range(n):\n",
    "    ax.annotate(str(i), xy=(xy_cords[i,0], xy_cords[i,1]+0.1))\n",
    "\n",
    "ax.scatter(xy_cords[:,0],xy_cords[:,1])\n",
    "for i in ruta.keys():\n",
    "    ax.plot(xy_cords[ruta[i],0], xy_cords[ruta[i],1], label = i)\n",
    "    ax.legend(loc='best')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The optimal distance is: 28.3 km\n"
     ]
    }
   ],
   "source": [
    "# Showing the optimal distance\n",
    "distance = np.sum(np.multiply(C, X.value))\n",
    "print('The optimal distance is:', np.round(distance,2), 'km')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be observed that depending on the location of the points that salesmen must travel, in terms of distance it may not be optimal to add more salesmen."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
