milk_collection/milk_collection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Milk Collection Problem\n",
    "\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "In this example, you’ll discover how mathematical optimization can be leveraged to solve a capacitated vehicle routing problem: the Milk Collection Problem. With only one tanker truck with limited capacity, you will need to determine the best possible route for the tanker to take to collect milk every day from a set of farms. It’s a complicated problem to solve, but mathematical optimization will help show you the way!\n",
    "\n",
    "This model is example 23 from the fifth edition of Model Building in Mathematical Programming by H. Paul Williams on pages 278-281 and 336-337.\n",
    "\n",
    "This modeling example is at the advanced level, where we assume that you know Python and the Gurobi Python API and you have advanced knowledge of building mathematical optimization models. Typically, the objective function and/or constraints of these examples are complex or require advanced features of the Gurobi Python API\n",
    "\n",
    "**Download the Repository** <br />\n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Problem Description\n",
    "\n",
    "A small milk processing company is committed to collecting milk from 20 farms and taking it back to the depot for processing. The company has one tanker lorry with the capacity to carry 80 000 liters of milk. Eleven of the farms are small and need a collection only every other day. The other nine farms need a collection every day. The positions of the farms in relation to the depot (numbered 1) are given in the following table together with their collection requirements.\n",
    "\n",
    "![farmCoordinates](farmCoordinates.PNG)\n",
    "\n",
    "The goal is to find the optimal route for the tanker lorry on each day, bearing in mind that it has to: \n",
    "1. Visit all the ‘every day’ farms. \n",
    "2. Visit some of the ‘every other day’ farms. \n",
    "3. Work within its capacity. \n",
    "\n",
    "On alternate days, it must again visit the ‘every day’ farms and also visit the ‘every other day’ farms not visited on the\n",
    "previous day."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Model Formulation\n",
    "\n",
    "### Sets and Indices\n",
    "$i, j \\in \\text{Farms} = \\{0,1,2, ..,20 \\}$: Indices and set of farms. The depot index number is 0.\n",
    "\n",
    "$\\text{everyDay} = \\{0,1,2, ..,9 \\} \\subset \\text{Farms}$: Farms that need to visit every day.\n",
    "\n",
    "$\\text{otherDay} = \\{10,11, 12, ..,20 \\} \\subset \\text{Farms}$: Farms that need to visit every other day.\n",
    "\n",
    "$k \\in K = \\{1,2 \\} $: Day type for farms that are visited every other day.\n",
    "\n",
    "$\\text{Edges}= \\{(i,j) \\in Farms \\times Farms \\}$: Set of allowed Edges\n",
    "\n",
    "$S_k \\in S  $: Subtour in route of day $k$.\n",
    "\n",
    "$G = (\\text{Farms} , \\text{Edges})$: A graph where the set $\\text{Farms}$ defines the set of nodes and the set $\\text{Edges}$ defines the set of edges. \n",
    "\n",
    "### Parameters \n",
    "\n",
    "$d_{i, j} \\in \\mathbb{R}^+$: Distance from farm $i$ to farm $j$, for all $(i, j) \\in \\text{Edges}$. \n",
    "\n",
    "Notice that the distance from farm $i$ to farm $j$ is the same as the distance from farm $j$ to farm $i$, i.e. $d_{i, j} = d_{j, i}$.\n",
    "\n",
    "$C \\in \\mathbb{R}^+$: The capacity of the tanker lorry.\n",
    "\n",
    "$R_i \\in \\mathbb{R}^+$: Milk collection requirements of farm $i$.\n",
    "\n",
    "### Decision Variables\n",
    "$x_{i, j, k} \\in \\{0, 1\\}$: This variable is equal to 1, if the tour on day $k$ connects directly farm $i$ with farm $j$. Otherwise, the decision variable is equal to zero.\n",
    "\n",
    "$y_{i, k} \\in \\{0, 1\\}$: This variable is equal to 1, if farm $i \\in otherDay$ is visited on the tour of day $k \\in K$. \n",
    "\n",
    "### Objective Function\n",
    "- **Shortest Route**. Minimize the total distance of both routes. \n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Min} \\quad Z = \\sum_{k \\in K} \\sum_{(i,j) \\in \\text{Edges}} \\frac{1}{2} d_{i,j} \\cdot x_{i,j,k}\n",
    "\\tag{0}\n",
    "\\end{equation}\n",
    "\n",
    "### Constraints \n",
    "- **Symmetry Constraints**. For each edge $(i,j)$, ensure that the farm $i$ and $j$ are connected, if the former is visited immediately before or after visiting the latter.\n",
    "\n",
    "\\begin{equation}\n",
    "x_{i, j, k} = x_{j, i, k} \\quad \\forall k \\in dayType, \\; (i, j) \\in Edges\n",
    "\\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "- **Entering and leaving an every day farm**. For each farm $i$, ensure that this farm is connected to two other farms. \n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{j: (i,j) \\in \\text{Edges}} x_{i,j,k} = 2 \\quad \\forall  i \\in everyDay, \\; k \\in dayType \n",
    "\\tag{2}\n",
    "\\end{equation}\n",
    "\n",
    "- **Entering and leaving an every other day farm**. For each farm $i$, ensure that this farm is connected to two other farms. \n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{j: (i,j) \\in \\text{Edges}} x_{i,j,k}  = 2 \\cdot y_{i, k} \\quad \\forall  i \\in otherDay, \\; k \\in dayType \n",
    "\\tag{3}\n",
    "\\end{equation}\n",
    "\n",
    "- **Tanker capacity**. The tanker has limited capacity.\n",
    "\\begin{equation}\n",
    "\\sum_{i \\in \\text{otherDay}} R_{i} \\cdot y_{i,k} \\leq C -\\sum_{i \\in everyDay} R_{i} \\quad \\forall  k \\in K \n",
    "\\tag{4}\n",
    "\\end{equation}\n",
    "\n",
    "- **Farms visited**. Limit on visiting some farms only every other day.\n",
    "\n",
    "\\begin{equation}\n",
    "y_{i,1} + y_{i,2}  = 1 \\quad \\forall  i \\in \\text{otherDay}\n",
    "\\tag{5}\n",
    "\\end{equation}\n",
    "\n",
    "- **Subtour elimination**. These constraints ensure that for each route there is no cycle. \n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{(i,j) \\in S_k}x_{i,j,k} \\leq |S_k|-1 \\quad \\forall  k \\in K, \\;   S_k \\in S\n",
    "\\tag{6}\n",
    "\\end{equation}\n",
    "\n",
    "Where the subset $S$ is the set of farms visited by the tour, and this tour is defined by the values of the decision variables in the LHS of the inequality."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module and other Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "import math\n",
    "from itertools import combinations\n",
    "\n",
    "# tested with Python 3.7.0 & Gurobi 9.1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input data  \n",
    "We define all the input data for the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dictionary to capture the farm coordinates (10 miles) and collection requirements (1,000).\n",
    "\n",
    "farms, coordinates, collect  = gp.multidict({\n",
    "    0: [(0,0),0],\n",
    "    1: [(-3,3),5],\n",
    "    2: [(1,11),4],\n",
    "    3: [(4,7),3],\n",
    "    4: [(-5,9),6],\n",
    "    5: [(-5,-2),7],\n",
    "    6: [(-4,-7),3],\n",
    "    7: [(6,0),4],\n",
    "    8: [(3,-6),6],\n",
    "    9: [(-1,-3),5],\n",
    "    10: [(0,-6),4],\n",
    "    11: [(6,4),7], \n",
    "    12: [(2,5),3],\n",
    "    13: [(-2,8),4],\n",
    "    14: [(6,10),5],\n",
    "    15: [(1,8),6],\n",
    "    16: [(-3,1),8],\n",
    "    17: [(-6,5),5],\n",
    "    18: [(2,9),7],\n",
    "    19: [(-6,-5),6],\n",
    "    20: [(5,-4),6]\n",
    "})\n",
    "\n",
    "# List of farm  including the depot\n",
    "farms = [*range(0,21)]\n",
    "\n",
    "# List of farms that are visited everyday\n",
    "everyDay = [*range(0,10)]\n",
    "\n",
    "# List of farms that are visited each other day\n",
    "otherDay = [*range(10,21)]\n",
    "\n",
    "# List of day types\n",
    "dayType = [1,2]\n",
    "\n",
    "# Tanker lorry capacity (1,000)\n",
    "tankerCap = 80\n",
    "\n",
    "# Every day farms requirements\n",
    "everyDayReq = 0\n",
    "for i in everyDay:\n",
    "    everyDayReq += collect[i]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data processing\n",
    "Here, we calculate the distance for each pair of farms and store it in a Python dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute pairwise distance matrix\n",
    "# numpy linalg norm = euclidean n=2\n",
    "\n",
    "def distance(city1, city2):\n",
    "    c1 = coordinates[city1]\n",
    "    c2 = coordinates[city2]\n",
    "    diff = (c1[0]-c2[0], c1[1]-c2[1])\n",
    "    return math.sqrt(diff[0]*diff[0]+diff[1]*diff[1])\n",
    "\n",
    "dist = {(c1, c2): distance(c1, c2) for c1, c2 in combinations(farms, 2)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Deployment\n",
    "We now determine the capacitated vehicle routing model for the milk collection problem by defining decision variables, constraints, and objective function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the model object m\n",
    "m = gp.Model('MilkCollection.lp')\n",
    "\n",
    "# Decision variables: \n",
    "\n",
    "# Edge variables = 1, if farm 'i' is adjacent to farm 'j' on the tour of day type 'k'.\n",
    "vars = m.addVars(dist, dayType, vtype=GRB.BINARY, name='x')\n",
    "\n",
    "# Other day variables = 1, if farm 'i' is visited on the tour of day type 'k'.\n",
    "other_var = m.addVars(otherDay, dayType, vtype=GRB.BINARY, name='y') \n",
    "\n",
    "# Symmetry constraints: copy the object (not a constraint)\n",
    "vars.update({(j,i,k):vars[i,j,k] for i,j,k in vars.keys()})\n",
    "\n",
    "# Every day constraints: two edges incident to an every day farm on tour of day type 'k'. \n",
    "m.addConstrs((vars.sum(i,'*',k) == 2 for i in everyDay for k in dayType  ), name='everyDay')\n",
    "\n",
    "# Other day constraints: two edges incident to an other day farm on tour of day type 'k'.\n",
    "m.addConstrs((vars.sum(i,'*',k) == 2*other_var[i,k] for i in otherDay for k in dayType ), name='otherDay')\n",
    "\n",
    "# Tanker capacity constraint.\n",
    "m.addConstrs(( gp.quicksum(collect[i]*other_var[i,k] for i in otherDay ) <= tankerCap-everyDayReq for k in dayType ),\n",
    "             name='tankerCap')\n",
    "\n",
    "# Other day farms are visited on day type 1 or 2.\n",
    "otherDayFarms = m.addConstrs((other_var.sum(i, '*') == 1 for i in otherDay), name='visited')\n",
    "\n",
    "# Avoid symmetric alternative solutions\n",
    "other_var[11,1].lb = 1\n",
    "\n",
    "# Objective function: minimize total distance travel\n",
    "m.setObjective(gp.quicksum(dist[i,j]*vars[i,j,k] for i,j in dist.keys() for k in dayType), GRB.MINIMIZE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding A Cycle\n",
    "The following function determines a cycle not connected to the depot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the edges from solution values, as a tuplelist for each dayType\n",
    "def selected(vals):\n",
    "    s = {k:gp.tuplelist() for k in dayType}\n",
    "    for i, j, k in vals.keys():\n",
    "        if vals[i,j,k] > 0.5:\n",
    "            s[k].append((i,j))\n",
    "    return s\n",
    "# Alternately, using comprehension syntax:\n",
    "#    return {k:gp.tuplelist((i, j) for i, j, k in vals.keys().select('*','*',k) if vals[i,j,k] > 0.5) for k in dayType}\n",
    "            \n",
    "# Given a tuplelist of edges, find the shortest subtour\n",
    "def subtour(edges):\n",
    "    nodes = set(i for e in edges for i in e)\n",
    "    unvisited = list(nodes)\n",
    "    cycle = list(nodes)\n",
    "    while unvisited:  \n",
    "        thiscycle = []\n",
    "        neighbors = unvisited\n",
    "        while neighbors:\n",
    "            current = neighbors[0]\n",
    "            thiscycle.append(current)\n",
    "            unvisited.remove(current)\n",
    "            neighbors = [j for i, j in edges.select(current, '*')\n",
    "                         if j in unvisited]\n",
    "        if len(thiscycle) <= len(cycle): # even if it's the same, we reuse it so that we get the final tour in order\n",
    "            cycle = thiscycle # New shortest subtour\n",
    "    return cycle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Callback Definition\n",
    "Subtour constraints prevent multiple loops in a tour. Because there are an exponential number of these constraints, we don't want to add them all to the model. Instead, we use a callback function to find violated subtour constraints and add them to the model as lazy constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Callback - use lazy constraints to eliminate sub-tours\n",
    "def subtourelim(model, where):\n",
    "    if where == GRB.Callback.MIPSOL:\n",
    "        # get edges selected in the current solution\n",
    "        vals = model.cbGetSolution(model._vars)\n",
    "        edges = selected(vals)\n",
    "        for k in dayType:\n",
    "            tour = subtour(edges[k])\n",
    "            if len(tour) < 0.5*len(edges[k]): # 0.5 due to symmetry: there are both edges i,j and j,i\n",
    "                # add subtour elimination constr. for farms visited that day\n",
    "                model.cbLazy(gp.quicksum(model._vars[i, j, k] for i, j in combinations(tour, 2))\n",
    "                             <= len(tour)-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve\n",
    "We can now optimize the model with the lazy subtour constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimize the model\n",
    "\n",
    "m.reset()\n",
    "m._vars = vars\n",
    "m.Params.lazyConstraints = 1\n",
    "m.optimize(subtourelim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "We print the optimal total distance traveled and the associated optimal routes for each day type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print optimal routes and distance traveled\n",
    "\n",
    "print(f\"The optimal distance traveled is: {10*round(m.objVal)} miles.\")\n",
    "\n",
    "vals = m.getAttr('X', vars)\n",
    "edges = selected(vals)\n",
    "\n",
    "for k in dayType:\n",
    "    tour = subtour(edges[k])\n",
    "    tour.append(0) # return to depot\n",
    "    print (\"Route for day type %i: %s\" % (k, \" -> \".join(map(str, tour))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}