optimization101/bike_share/bike_rebalancing_complete.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "406c2750",
   "metadata": {},
   "source": [
    "# Bike Share Rebalancing With Mathematical Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a7e7404",
   "metadata": {},
   "source": [
    "Bike share systems have become an effective commuting method globally for everyday urban dwellers as well as tourists.\n",
    "\n",
    "Citi-Bike in NYC being the largest Bike-Sharing network had 1,588 active stations and 25,575 active bikes in July 2022.\n",
    "\n",
    "Over 3 million rides were completed July 2022 that cover NYC/Hoboken/Jersey City, with around 150,000 active annual members.\n",
    "\n",
    "During rush hours there are many bike stations that have a high demand for bikes, which means their out-flow of bikes is greater than their in-flow in these stations. \n",
    "\n",
    "Meanwhile there are stations that have a high demand for docks (riders return their bikes to these stations) which means their in-flow of bikes is greater than their out-flow.\n",
    "\n",
    "Lack of available bikes or docks in high-demand stations can cause major imbalance in the bike sharing network and result in customer dissatisfaction and lost revenue.\n",
    "\n",
    "To tackle this problem, bikes are relocated between stations to create a balance between supply and demand."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fdb0fdf",
   "metadata": {},
   "source": [
    "## Problem Statement and Solution Approach"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2142854",
   "metadata": {},
   "source": [
    "Using historical Citi-bike data in NYC and Jersey area during July 2022, we like to know:\n",
    "- What is the demand for bikes per hour at each station during the first week of August?\n",
    "- Knowing the demand, how can we minimize loss of sale?\n",
    "\n",
    "Loss of sale is caused by lack of bikes when customers demand them. So, bikes should be transferred from stations with higher in-flow of bikes to those with higher out-flow of bikes. \n",
    "\n",
    "So, first, number of bikes to be added to or removed from each station during each hour should be determined. Then, the physical transfer of bikes between stations should be scheduled. \n",
    "\n",
    "In this notebook, we'll focus on the first part and at the end, discuss how the second part can be solved.\n",
    "We'll use a mixture of Machine Learning (ML) and Mathematical Optimization (MO) to solve this problem. \n",
    "\n",
    "**Solution Approach**\n",
    "The solution approach is comprised of two steps:\n",
    "- **Step 1**: We use the historical Citi-bike data in NYC and Jersey area during July 2022 and use an ML model to predict the number of in-flow and out-flow of bikes per hour at each station for the first week of August. This is done in [predict_bike_flow](predict_bike_flow.ipynb) Notebook.\n",
    "- **Step 2**: We use an MO model to decide how many bikes should be added to or removed from each station during each hour so that the total loss of sale is minimized.\n",
    "\n",
    "To ensure that everyone can run the notebook with the Gurobi restricted license, we reduce the size of the data. To achieve that, we focus on the top 50 stations during the morning rush hours (7 am to 9 am).\n",
    "\n",
    "The top stations are chosen using the PageRank algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed4f7ee5",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Install Required Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af590656",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-16T02:10:37.955737Z",
     "start_time": "2022-11-16T02:10:35.495638Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "%pip install gurobipy\n",
    "%pip install pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8bd6f69",
   "metadata": {},
   "source": [
    "# Import Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4036b5a7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:20:41.899618Z",
     "start_time": "2022-11-18T00:20:41.431601Z"
    }
   },
   "outputs": [],
   "source": [
    "import datetime\n",
    "import gurobipy as gp\n",
    "import pandas as pd\n",
    "from gurobipy import GRB"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c165d441",
   "metadata": {},
   "source": [
    "# Optimization Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58e8f88d",
   "metadata": {},
   "source": [
    "## Problem Definition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0816abcb",
   "metadata": {},
   "source": [
    "We want to minimize the total loss of sale. Loss of sale at each station and in each hour can be defined as the difference between the total demand of bikes (number of bikes that start their trip from the station) and total supply of bikes.\n",
    "\n",
    "Total supply is comprised of number of bikes that end their trip at the station plus all the existing bikes at the station (a.k.a inventory) plus number of bikes that are added or removed from that station in that hour through some bike transfers. \n",
    "\n",
    "**Assumptions:**\n",
    "- Inventory at the beginning of first hour (in our case, hour 7) is zero.\n",
    "- At any given hour, we have access to a limited number of bikes that can be added to the stations in hope of helping reduce the imbalance without yet transferring the bikes between stations (since this analysis is during morning rush hours)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79287692",
   "metadata": {},
   "source": [
    "## Load Required Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3fd07b80",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:20:44.203583Z",
     "start_time": "2022-11-18T00:20:44.173553Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>capacity</th>\n",
       "      <th>lat</th>\n",
       "      <th>lon</th>\n",
       "      <th>region</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>station</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Cleveland Pl &amp; Spring St</th>\n",
       "      <td>33</td>\n",
       "      <td>40.722104</td>\n",
       "      <td>-73.997249</td>\n",
       "      <td>71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Grand Army Plaza &amp; Central Park S</th>\n",
       "      <td>96</td>\n",
       "      <td>40.764397</td>\n",
       "      <td>-73.973715</td>\n",
       "      <td>71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Broadway &amp; E 14 St</th>\n",
       "      <td>114</td>\n",
       "      <td>40.734546</td>\n",
       "      <td>-73.990741</td>\n",
       "      <td>71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lafayette St &amp; E 8 St</th>\n",
       "      <td>47</td>\n",
       "      <td>40.730207</td>\n",
       "      <td>-73.991026</td>\n",
       "      <td>71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Norfolk St &amp; Broome St</th>\n",
       "      <td>59</td>\n",
       "      <td>40.717227</td>\n",
       "      <td>-73.988021</td>\n",
       "      <td>71</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                   capacity        lat        lon  region\n",
       "station                                                                  \n",
       "Cleveland Pl & Spring St                 33  40.722104 -73.997249      71\n",
       "Grand Army Plaza & Central Park S        96  40.764397 -73.973715      71\n",
       "Broadway & E 14 St                      114  40.734546 -73.990741      71\n",
       "Lafayette St & E 8 St                    47  40.730207 -73.991026      71\n",
       "Norfolk St & Broome St                   59  40.717227 -73.988021      71"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stations = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/bike_share/top_stations.csv', index_col='station')\n",
    "# run locally\n",
    "# stations = pd.read_csv('top_stations.csv', index_col='station')\n",
    "stations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2d6872ea",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:20:45.797953Z",
     "start_time": "2022-11-18T00:20:45.779127Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>station</th>\n",
       "      <th>datetime</th>\n",
       "      <th>start_forecast</th>\n",
       "      <th>end_forecast</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>2022-08-01 00:00:00</td>\n",
       "      <td>4.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>2022-08-01 01:00:00</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>2022-08-01 02:00:00</td>\n",
       "      <td>4.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>2022-08-01 03:00:00</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>2022-08-01 04:00:00</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           station            datetime  start_forecast  end_forecast\n",
       "0  1 Ave & E 62 St 2022-08-01 00:00:00             4.0           5.0\n",
       "1  1 Ave & E 62 St 2022-08-01 01:00:00             2.0           3.0\n",
       "2  1 Ave & E 62 St 2022-08-01 02:00:00             4.0           2.0\n",
       "3  1 Ave & E 62 St 2022-08-01 03:00:00             2.0           2.0\n",
       "4  1 Ave & E 62 St 2022-08-01 04:00:00             2.0           2.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stations_flow = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/bike_share/stations_flow.csv')\n",
    "# run locally\n",
    "# stations_flow = pd.read_csv('stations_flow.csv')\n",
    "stations_flow['datetime'] = pd.to_datetime(stations_flow['datetime'])\n",
    "stations_flow.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a791ab30",
   "metadata": {},
   "source": [
    "The `stations_flow` data contains the prediction for the first 5 days of August 2022 and during all the hours. Our analysis is for morning rush hours, between 7 to 9 am. Also, we can run our MO model daily. For now, we only focus on the first day but at the end, we will show the full model and how it can be run daily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "190b9bce",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:20:48.030622Z",
     "start_time": "2022-11-18T00:20:48.004820Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>datetime</th>\n",
       "      <th>start_forecast</th>\n",
       "      <th>end_forecast</th>\n",
       "      <th>date</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>station</th>\n",
       "      <th>time</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">1 Ave &amp; E 62 St</th>\n",
       "      <th>7</th>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>11.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>14.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>17.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">1 Ave &amp; E 68 St</th>\n",
       "      <th>7</th>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>17.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>22.0</td>\n",
       "      <td>41.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                datetime  start_forecast  end_forecast  \\\n",
       "station         time                                                     \n",
       "1 Ave & E 62 St 7    2022-08-01 07:00:00            11.0          10.0   \n",
       "                8    2022-08-01 08:00:00            14.0          14.0   \n",
       "                9    2022-08-01 09:00:00            17.0          20.0   \n",
       "1 Ave & E 68 St 7    2022-08-01 07:00:00            17.0          32.0   \n",
       "                8    2022-08-01 08:00:00            22.0          41.0   \n",
       "\n",
       "                            date  \n",
       "station         time              \n",
       "1 Ave & E 62 St 7     2022-08-01  \n",
       "                8     2022-08-01  \n",
       "                9     2022-08-01  \n",
       "1 Ave & E 68 St 7     2022-08-01  \n",
       "                8     2022-08-01  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Pandas will give a few SettingWithCopyWarning here when new columns are created. \n",
    "# They are false alarms. So, we suppress them.\n",
    "pd.options.mode.chained_assignment = None\n",
    "morning_flow = stations_flow[stations_flow['datetime'].dt.hour.between(7, 9)]\n",
    "morning_flow['date'] = morning_flow['datetime'].dt.date\n",
    "morning_flow['time'] = morning_flow['datetime'].dt.hour\n",
    "# For now, let's run the MO model for the first date: 08/01/2022\n",
    "flow_df = morning_flow.loc[morning_flow['date'] == datetime.date(2022, 8, 1)]\n",
    "flow_df.set_index(['station', 'time'], inplace=True)\n",
    "flow_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75881443",
   "metadata": {},
   "source": [
    "## Problem Formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b8229be",
   "metadata": {},
   "source": [
    "Let's define our notations for the MO model. We want to run the model for every hour and every station. So, let's define the following two sets:\n",
    "\n",
    "**Sets**\n",
    "- $I\\quad$: Set of stations\n",
    "- $T\\quad$: Set of hours \n",
    "\n",
    "We also have the following information from `stations` and `flow_df` dataframes:\n",
    "\n",
    "**Parameters**\n",
    "- $e_{i,t}\\quad$: number of bikes that end their trip at station $i$ at hour $t$ (a.k.a. supply)\n",
    "- $s_{i,t}\\quad$: number of bikes that start their trip at station $i$ at hour $t$ (a.k.a. demand)\n",
    "- $c_{i}\\quad$: capacity of station $i$\n",
    "\n",
    "We know it's not easy to transfer bikes between stations during rush hours and heavy traffic. To mitigate loss due to unavailable bikes at high-demand stations, we assume there is a small reserve of bikes available at the beginning of an hour that we can allocate to the stations. We show that with $N$: \n",
    "- $N\\quad$: Number of bikes at hand that we can assign to stations at a given hour"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e8140631",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:20:51.414665Z",
     "start_time": "2022-11-18T00:20:51.401524Z"
    }
   },
   "outputs": [],
   "source": [
    "num_bikes = 25  # N: Number of bikes at hand that we can assign to stations at a given hour\n",
    "station_names = list(stations.index)  # set I\n",
    "time_rng = morning_flow['time'].drop_duplicates().values  # set T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "68f17fc6-abf3-470d-9848-a94d98c1f7cc",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:20:52.302162Z",
     "start_time": "2022-11-18T00:20:52.282123Z"
    }
   },
   "outputs": [],
   "source": [
    "station_time = flow_df.index  # pair of (i,t) index\n",
    "start_forecast = flow_df.start_forecast  # s\n",
    "end_forecast = flow_df.end_forecast  # e\n",
    "capacity = stations.capacity  # c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fcf3cf8",
   "metadata": {},
   "source": [
    "## Decision Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d21f6f93",
   "metadata": {},
   "source": [
    "First of all, we like to find out how many bikes should be added to or removed from each station during each hour. We show number of bikes added with $y_{i,t}$ and those removed with $z_{i,t}$:\n",
    "- $y_{i,t}\\quad$: number of bikes to be added to station $i$ at hour $t$\n",
    "- $z_{i,t}\\quad$: number of bikes to be removed from station $i$ at hour $t$\n",
    "\n",
    "Another variable is the inventory of bikes at each station and at the beginning of each hour. Since the inventory depends on the value of $y_{i,t}$ and $z_{i,t}$, it's another decision variable. We show that with $q_{i,t}$:\n",
    "- $q_{i,t}\\quad$: inventory of bikes at station $i$ at the beginning of hour $t$\n",
    "\n",
    "For simplicity, we assume that there is no inventory at the beginning of the first hour. In our case, this means no initial inventory at 7 am.\n",
    "\n",
    "The goal of the model is to reduce the lost sale at each station per hour. This value, also depends on the value of decision variables $y_{i,t}$ and $z_{i,t}$ and as a result, is a decision variable itself. We show that with $l_{i,t}$:\n",
    "- $l_{i,t}\\quad$: lost sale at station $i$ at hour $t$\n",
    "\n",
    "With these initial decision variables, we can start the MO model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "18f9b67f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:20:57.991204Z",
     "start_time": "2022-11-18T00:20:57.970149Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "mdl = gp.Model('bike_rebalancing')\n",
    "# Variables\n",
    "\n",
    "y = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='y')\n",
    "# y = mdl.addVars(station_names, time_rng, lb=0, vtype=GRB.CONTINUOUS, name='y')  # alternatively\n",
    "z = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='z')\n",
    "q = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='q')\n",
    "l = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='l')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53a33172",
   "metadata": {},
   "source": [
    "## Constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "225d8693",
   "metadata": {},
   "source": [
    "First, we set up lowerbound and upperbound values for the number of bikes that can be added to or removed from a station in each hour.\n",
    "\n",
    "In each hour, we have $s_{i,t}$ bikes that start their trip from the station and $e_{i,t}$ bikes that end their trip at the station. If $s_{i,t} \\ge e_{i,t}$, demand is exceeding supply. In this case, we may choose to add some bikes to this station to reduce the loss (we can also add nothing). One thing to remember is that we cannot add more bikes than the station's capacity.\n",
    "\n",
    "On the other hand, if $s_{i,t} \\le e_{i,t}$, we may choose to remove some of the excess bikes from that station (we can also remove nothing). However, if $e_{i,t} \\ge s_{i,t} + c_i$, then we'll have more bikes arriving to the station than even the station capacity. In that case, we must remove some bikes to avoid overflow.\n",
    "\n",
    "The above descriptions, help us to define the bounds on $y_{i,t}$ and $z_{i,t}$.\n",
    "For adding bikes, the bounds are:\n",
    "\n",
    "\\begin{align}\n",
    "&0 \\le y_{i,t} \\le c_i &\\quad \\forall i \\in I, t \\in T \\tag{1}\\\\\n",
    "\\end{align}\n",
    "\n",
    "For removing bikes, the upperbound is:\n",
    "\n",
    "\\begin{align}\n",
    "&z_{i,t} \\le \\max(0, e_{i,t} - s_{i,t}) &\\quad \\forall i \\in I, t \\in T \\tag{2}\\\\\n",
    "\\end{align}\n",
    "\n",
    "And the lowerbound is:\n",
    "\n",
    "\\begin{align}\n",
    "&z_{i,t} \\ge  e_{i,t} - s_{i,t} - c_i &\\quad \\forall i \\in I, t \\in T \\tag{3}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fe03af94-b72d-440c-b81d-e1de82eed7c6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:09.340868Z",
     "start_time": "2022-11-18T00:21:09.309277Z"
    }
   },
   "outputs": [],
   "source": [
    "# UB of y\n",
    "mdl.addConstrs((y[i, t] <= capacity[i] for i, t in station_time), 'ub_y')\n",
    "\n",
    "# LB of z\n",
    "mdl.addConstrs(\n",
    "    (z[i, t] >= (end_forecast[i, t] - start_forecast[i, t] - capacity[i])\n",
    "     for i, t in station_time), 'lb_z')\n",
    "#     for i in station_names for t in time_rng), 'lb_z')\n",
    "\n",
    "# UB of z\n",
    "mdl.addConstrs(\n",
    "        (z[i, t] <= max(0, end_forecast[i, t] - start_forecast[i, t])\n",
    "         for i, t in station_time), 'ub_z');  # add \";\" here to stop the cell to print the constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db534a31",
   "metadata": {},
   "source": [
    "Next, we set up the initial inventory (i.e. inventory at hour $t_0$) to be 0. \n",
    "\n",
    "\\begin{align}\n",
    "&q_{i,t_0} = 0 &\\quad \\forall i \\in I \\tag{4}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e0b73d5f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:11.635769Z",
     "start_time": "2022-11-18T00:21:11.621260Z"
    }
   },
   "outputs": [],
   "source": [
    "t0 = 7\n",
    "# setting up initial inventory\n",
    "mdl.addConstrs((q[i, t0] == 0 for i in stations.index), name='initial_inv');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28688fc6",
   "metadata": {},
   "source": [
    "Our next constraint is the definition of the inventory at a station at the beginning of an hour.\n",
    "The inventory is defined as the difference between all the bikes that are added to a station minus all the bikes that are removed. Since it's possible that demand of bikes at a station exceeds the total supply of bikes (which will cause loss of sale), we need to ensure that inventory at a station is a non-negative number.\n",
    "\n",
    "\\begin{align}\n",
    "&q_{i,t} = \\max(0, e_{i,t-1} + y_{i,t-1} + q_{i,t-1} - s_{i,t-1} - z_{i,t-1}) &\\quad \\forall i \\in I, t \\in T/t_0 \\tag{5}\\\\\n",
    "\\end{align}\n",
    "\n",
    "To write this constraint in Gurobi, we can take advantage of Gurobi's [`max_()`](https://www.gurobi.com/documentation/9.5/refman/py_max_.html) function. This function accepts a list of decision variables, and if desired, a constant. Our constant is 0. But $e_{i,t-1} + y_{i,t-1} + q_{i,t-1} - s_{i,t-1} - z_{i,t-1}$ is a linear expression and not a decision variable. However, the fix is very simple. We can define an auxiliary variable that is equal to the linear expression and then use that auxiliary variable in the inventory's definition. In other words, we first create $a_{i,t}$ as a new decision variable:\n",
    "\n",
    "- $a_{i,t}\\quad$: auxiliary variable\n",
    "\n",
    "Then, we define $a_{i,t}$ as:\n",
    "\n",
    "\\begin{align}\n",
    "&a_{i,t} = e_{i,t-1} + y_{i,t-1} + q_{i,t-1} - s_{i,t-1} - z_{i,t-1} &\\quad \\forall i \\in I, t \\in T/t_0 \\tag{6}\\\\\n",
    "\\end{align}\n",
    "\n",
    "Equation 5 can then be simplified to:\n",
    "\n",
    "\\begin{align}\n",
    "&q_{i,t} = \\max(0, a_{i,t}) &\\quad \\forall i \\in I, t \\in T/t_0 \\tag{7}\\\\\n",
    "\\end{align}\n",
    "\n",
    "We replace equation 5 with equations 6 and 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ba4e366d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:13.929281Z",
     "start_time": "2022-11-18T00:21:13.907145Z"
    }
   },
   "outputs": [],
   "source": [
    "a = mdl.addVars(station_time, lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS, name='a')  # auxiliary variable\n",
    "# defintion of auxiliary variable\n",
    "mdl.addConstrs((a[i, t] == end_forecast[i, t - 1] + y[i, t - 1] + q[i, t - 1]\n",
    "                - start_forecast[i, t - 1] - z[i, t - 1]\n",
    "                for i, t in station_time if t != t0), name='aux_def')\n",
    "# definition of inventory\n",
    "mdl.addConstrs((q[i, t] == gp.max_(a[i, t], 0) for i, t in station_time if t!=t0), name='inv_def');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "718a8ac0",
   "metadata": {},
   "source": [
    "The next constraint is an upperbound on the inventory at each station. We need to ensure the inventory does not exceed the capacity of the station.\n",
    "\n",
    "\\begin{align}\n",
    "&q_{i,t} \\le c_i &\\quad \\forall i \\in I, t \\in T \\tag{8}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3d1a3c0a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:15.845620Z",
     "start_time": "2022-11-18T00:21:15.835238Z"
    }
   },
   "outputs": [],
   "source": [
    "# UB of inventory\n",
    "mdl.addConstrs((q[i, t] <= capacity[i] for i, t in station_time), 'ub_inv');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdfeb12a",
   "metadata": {},
   "source": [
    "Next, we need to define how the loss of sale is calculated. Loss of sale is the difference between the demand and all the supply of bikes at a station. \n",
    "- Demand of bikes are all the bikes that leave a station, in any shape or form. So, what are the demands?\n",
    "- Supply of bikes are all the bikes that arrive at a station, in any shape or form. So, what are the supplies?\n",
    "\n",
    "Of course, if the supply is greater than the demand, there is no loss. So, we need to ensure that loss only considers non-negative values. This can be achieved by:\n",
    "\n",
    "\\begin{align}\n",
    "&l_{i,t} = \\max(0, s_{i,t} + z_{i,t} - e_{i,t} - y_{i,t} - q_{i,t}) &\\quad \\forall i \\in I, t \\in T \\tag{9}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4a467dca",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:17.785983Z",
     "start_time": "2022-11-18T00:21:17.760364Z"
    }
   },
   "outputs": [],
   "source": [
    "# loss definition\n",
    "mdl.addConstrs(\n",
    "    (l[i, t] >= start_forecast[i, t] + z[i, t]\n",
    "     - end_forecast[i, t] - y[i, t] - q[i, t]\n",
    "     for i, t in station_time), 'loss_def');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf04748f",
   "metadata": {},
   "source": [
    "We assumed that we have a small reserve of bikes at the beginning of each hour to allocate to stations. This limit is on the total number of bikes added to the stations.\n",
    "\n",
    "\\begin{align}\n",
    "&\\sum_{i} y_{i,t} \\le N &\\quad \\forall t \\in T \\tag{10}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "61f52e65",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:19.314712Z",
     "start_time": "2022-11-18T00:21:19.301113Z"
    }
   },
   "outputs": [],
   "source": [
    "# limit on number of bikes added\n",
    "mdl.addConstrs((y.sum('*', t) <= num_bikes for t in time_rng), name='total_bikes');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9e0d5df",
   "metadata": {},
   "source": [
    "## Objective"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cafeb700",
   "metadata": {},
   "source": [
    "The objective is to minimize total loss of sale. \n",
    "\n",
    "$$\\min \\sum_{i,t} l_{i,t}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a7352d11",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:21.469828Z",
     "start_time": "2022-11-18T00:21:21.461236Z"
    }
   },
   "outputs": [],
   "source": [
    "mdl.setObjective(l.sum(), GRB.MINIMIZE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93782ff9",
   "metadata": {},
   "source": [
    "We can now tell Gurobi that the model is complete and it can solve the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9b4669ff",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:22.565397Z",
     "start_time": "2022-11-18T00:21:22.500896Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 813 rows, 675 columns and 1620 nonzeros\n",
      "Model fingerprint: 0xecda0867\n",
      "Model has 90 general constraints\n",
      "Variable types: 675 continuous, 0 integer (0 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 1e+02]\n",
      "Presolve removed 750 rows and 580 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 63 rows, 95 columns, 201 nonzeros\n",
      "Variable types: 81 continuous, 14 integer (14 binary)\n",
      "Found heuristic solution: objective 40.0000000\n",
      "\n",
      "Root relaxation: objective 1.368348e+01, 62 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   13.68348    0    9   40.00000   13.68348  65.8%     -    0s\n",
      "     0     0   33.00000    0    6   40.00000   33.00000  17.5%     -    0s\n",
      "     0     0   33.00000    0    6   40.00000   33.00000  17.5%     -    0s\n",
      "     0     0   37.03846    0    2   40.00000   37.03846  7.40%     -    0s\n",
      "     0     0   39.00000    0    1   40.00000   39.00000  2.50%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 1\n",
      "  Implied bound: 3\n",
      "  MIR: 8\n",
      "  Flow cover: 14\n",
      "  Relax-and-lift: 2\n",
      "\n",
      "Explored 1 nodes (113 simplex iterations) in 0.04 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 40 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 4.000000000000e+01, best bound 4.000000000000e+01, gap 0.0000%\n"
     ]
    }
   ],
   "source": [
    "mdl.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9a1a8f4",
   "metadata": {},
   "source": [
    "## Post Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1f1b8fd1",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:26.127808Z",
     "start_time": "2022-11-18T00:21:26.037204Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total Loss : 40.0\n",
      "Total number of bikes added in each hour:\n",
      " time\n",
      "7    25\n",
      "8    25\n",
      "9    25\n",
      "Name: bikes_added, dtype: int64\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>station</th>\n",
       "      <th>time</th>\n",
       "      <th>datetime</th>\n",
       "      <th>start_forecast</th>\n",
       "      <th>end_forecast</th>\n",
       "      <th>date</th>\n",
       "      <th>capacity</th>\n",
       "      <th>bikes_added</th>\n",
       "      <th>bikes_removed</th>\n",
       "      <th>loss_sale</th>\n",
       "      <th>beginning_inventory</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>11.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>14.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>17.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>17.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0</td>\n",
       "      <td>15</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>22.0</td>\n",
       "      <td>41.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0</td>\n",
       "      <td>19</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>28.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>9.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>18.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>18.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>10 Ave &amp; W 14 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>55</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            station  time            datetime  start_forecast  end_forecast  \\\n",
       "0   1 Ave & E 62 St     7 2022-08-01 07:00:00            11.0          10.0   \n",
       "1   1 Ave & E 62 St     8 2022-08-01 08:00:00            14.0          14.0   \n",
       "2   1 Ave & E 62 St     9 2022-08-01 09:00:00            17.0          20.0   \n",
       "3   1 Ave & E 68 St     7 2022-08-01 07:00:00            17.0          32.0   \n",
       "4   1 Ave & E 68 St     8 2022-08-01 08:00:00            22.0          41.0   \n",
       "5   1 Ave & E 68 St     9 2022-08-01 09:00:00            28.0          48.0   \n",
       "6   1 Ave & E 78 St     7 2022-08-01 07:00:00             9.0           7.0   \n",
       "7   1 Ave & E 78 St     8 2022-08-01 08:00:00            18.0           5.0   \n",
       "8   1 Ave & E 78 St     9 2022-08-01 09:00:00            18.0           7.0   \n",
       "9  10 Ave & W 14 St     7 2022-08-01 07:00:00             3.0           5.0   \n",
       "\n",
       "         date  capacity  bikes_added  bikes_removed  loss_sale  \\\n",
       "0  2022-08-01        45            0              0          1   \n",
       "1  2022-08-01        45            0              0          0   \n",
       "2  2022-08-01        45            0              0          0   \n",
       "3  2022-08-01        62            0             15          0   \n",
       "4  2022-08-01        62            0             19          0   \n",
       "5  2022-08-01        62            0              0          0   \n",
       "6  2022-08-01        57            7              0          0   \n",
       "7  2022-08-01        57            8              0          0   \n",
       "8  2022-08-01        57           11              0          0   \n",
       "9  2022-08-01        55            0              2          0   \n",
       "\n",
       "   beginning_inventory  \n",
       "0                    0  \n",
       "1                    0  \n",
       "2                    0  \n",
       "3                    0  \n",
       "4                    0  \n",
       "5                    0  \n",
       "6                    0  \n",
       "7                    5  \n",
       "8                    0  \n",
       "9                    0  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame()\n",
    "if mdl.status == GRB.Status.OPTIMAL:\n",
    "    df = flow_df.copy()\n",
    "    df = df.merge(stations[['capacity']], left_on='station', right_index=True)\n",
    "    df[['bikes_added', 'bikes_removed', 'loss_sale', 'beginning_inventory']] = 0\n",
    "    for k, v in y.items():\n",
    "        df.loc[k, 'bikes_added'] = v.x\n",
    "        df.loc[k, 'bikes_removed'] = z[k].x\n",
    "        df.loc[k, 'beginning_inventory'] = q[k].x\n",
    "        df.loc[k, 'loss_sale'] = l[k].x\n",
    "    df.reset_index(inplace=True)\n",
    "    print(f'Total Loss : {mdl.objVal}')\n",
    "    total_bikes_added = df.groupby('time')['bikes_added'].sum()\n",
    "    print(f'Total number of bikes added in each hour:\\n {total_bikes_added}')\n",
    "else:\n",
    "    print('Could not find a feasible solution!')\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a926b3c9",
   "metadata": {},
   "source": [
    "## Putting it all together"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72c2cba8",
   "metadata": {},
   "source": [
    "**Sets**\n",
    "- $I\\quad$: Set of stations\n",
    "- $T\\quad$: Set of hours \n",
    "\n",
    "**Parameters**\n",
    "- $e_{i,t}\\quad$: number of bikes that end their trip at station $i$ at hour $t$ (a.k.a. supply)\n",
    "- $s_{i,t}\\quad$: number of bikes that start their trip at station $i$ at hour $t$ (a.k.a. demand)\n",
    "- $c_{i}\\quad$: capacity of station $i$\n",
    "- $N\\quad$: Number of bikes at hand that we can assign to different stations at a given hour\n",
    "\n",
    "**Variables**\n",
    "- $y_{i,t}\\quad$: number of bikes to be added to station $i$ at hour $t$\n",
    "- $z_{i,t}\\quad$: number of bikes to be removed from station $i$ at hour $t$\n",
    "- $q_{i,t}\\quad$: inventory of bikes at station $i$ at the beginning of hour $t$\n",
    "- $l_{i,t}\\quad$: lost sale at station $i$ at hour $t$\n",
    "- $a_{i,t}\\quad$: auxiliary variable needed to define inventory's relationship with other variables. We use this in the code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78b9555b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-10-19T19:10:20.288414Z",
     "start_time": "2022-10-19T19:10:20.268715Z"
    }
   },
   "source": [
    "\\begin{align}\n",
    "&\\min \\sum_{i,t} l_{i,t}&\\\\\n",
    "\\mbox{s.t: }\\\\\n",
    "&y_{i,t} \\le c_i &\\quad \\forall i \\in I, t \\in T \\tag{1}\\\\\n",
    "&z_{i,t} \\ge \\max(0, e_{i,t} - s_{i,t}) &\\quad \\forall i \\in I, t \\in T \\tag{2}\\\\\n",
    "&z_{i,t} \\le  e_{i,t} - s_{i,t} - c_i &\\quad \\forall i \\in I, t \\in T \\tag{3}\\\\\n",
    "&q_{i,t_0} = 0 &\\quad \\forall i \\in I \\tag{4}\\\\\n",
    "&a_{i,t} = e_{i,t-1} + y_{i,t-1} + q_{i,t-1} - s_{i,t-1} - z_{i,t-1} &\\quad \\forall i \\in I, t \\in T/t_0 \\tag{5}\\\\\n",
    "&q_{i,t} = \\max(0, a_{i,t}) &\\quad \\forall i \\in I, t \\in T/t_0 \\tag{6}\\\\\n",
    "&q_{i,t} \\le c_i &\\quad \\forall i \\in I, t \\in T \\tag{7}\\\\\n",
    "&l_{i,t} = \\max(0, s_{i,t} + z_{i,t} - e_{i,t} - y_{i,t} - q_{i,t}) &\\quad \\forall i \\in I, t \\in T \\tag{8}\\\\\n",
    "&\\sum_{i} y_{i,t} \\le N &\\quad \\forall t \\in T \\tag{9}\\\\\n",
    "&y_{i,t}, z_{i,t}, q_{i,t}, l_{i,t} \\ge0 &\\quad \\forall i \\in I, t \\in T \\tag{10}\\\\\n",
    "&a_{i,t} \\quad \\mbox{urs} &\\quad \\forall i \\in I, t \\in T \\tag{11}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "63a77b7e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:30.057450Z",
     "start_time": "2022-11-18T00:21:30.037195Z"
    }
   },
   "outputs": [],
   "source": [
    "def bike_rebalancing(flow_df, num_bikes):\n",
    "    station_time = flow_df.index  # pair of (i,t) index\n",
    "    start_forecast = flow_df.start_forecast  # s\n",
    "    end_forecast = flow_df.end_forecast  # e\n",
    "    capacity = stations.capacity  # c\n",
    "    mdl = gp.Model('bike_rebalancing')\n",
    "\n",
    "    # Variables\n",
    "    y = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='y')\n",
    "    z = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='z')\n",
    "    q = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='q')\n",
    "    l = mdl.addVars(station_time, lb=0, vtype=GRB.CONTINUOUS, name='l')\n",
    "    a = mdl.addVars(station_time, lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS, name='a')\n",
    "\n",
    "    # Constraints\n",
    "    mdl.addConstrs((y[i, t] <= capacity[i] for i, t in station_time), 'ub_y')\n",
    "\n",
    "    mdl.addConstrs(\n",
    "        (z[i, t] >= (end_forecast[i, t] - start_forecast[i, t] - capacity[i])\n",
    "         for i, t in station_time), 'lb_z')\n",
    "    \n",
    "    mdl.addConstrs(\n",
    "        (z[i, t] <= max(0, end_forecast[i, t] - start_forecast[i, t])\n",
    "         for i, t in station_time), 'ub_z')\n",
    "\n",
    "    t0 = 7\n",
    "    mdl.addConstrs((q[i, t0] == 0 for i in stations.index), name='initial_inv')\n",
    "\n",
    "    mdl.addConstrs((a[i, t] == end_forecast[i, t - 1] + y[i, t - 1] + q[i, t - 1]\n",
    "                    - start_forecast[i, t - 1] - z[i, t - 1]\n",
    "                    for i, t in station_time if t != t0), name='aux_def')\n",
    "\n",
    "    mdl.addConstrs((q[i, t] == gp.max_(a[i, t], 0) for i, t in station_time if t != t0), name='inv_def')\n",
    "\n",
    "    mdl.addConstrs((q[i, t] <= capacity[i] for i, t in station_time), 'ub_inv')\n",
    "\n",
    "    mdl.addConstrs((l[i, t] >= start_forecast[i, t] + z[i, t]\n",
    "                    - end_forecast[i, t] - y[i, t] - q[i, t] for i, t in station_time), 'loss_def')\n",
    "\n",
    "    mdl.addConstrs((y.sum('*', t) <= num_bikes for t in time_rng), name='total_bikes')\n",
    "\n",
    "    # Objectives\n",
    "    mdl.setObjective(l.sum(), GRB.MINIMIZE)\n",
    "    mdl.optimize()\n",
    "\n",
    "    # create output\n",
    "    df = pd.DataFrame()\n",
    "    if mdl.status == GRB.Status.OPTIMAL:\n",
    "        df = flow_df.copy()\n",
    "        df = df.merge(stations[['capacity']], left_on='station', right_index=True)\n",
    "        df[['bikes_added', 'bikes_removed', 'loss_sale', 'beginning_inventory']] = 0\n",
    "        for k, v in y.items():\n",
    "            df.loc[k, 'bikes_added'] = v.x\n",
    "            df.loc[k, 'bikes_removed'] = z[k].x\n",
    "            df.loc[k, 'beginning_inventory'] = q[k].x\n",
    "            df.loc[k, 'loss_sale'] = l[k].x\n",
    "        df.reset_index(inplace=True)\n",
    "        print(f'Total Loss : {mdl.objVal}')\n",
    "        total_bikes_added = df.groupby('time')['bikes_added'].sum()\n",
    "        print(f'Total number of bikes added in each hour:\\n {total_bikes_added}')\n",
    "    else:\n",
    "        print('Could not find a feasible solution!')\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "6ad4fea9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:40.624979Z",
     "start_time": "2022-11-18T00:21:39.730487Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 813 rows, 675 columns and 1620 nonzeros\n",
      "Model fingerprint: 0xecda0867\n",
      "Model has 90 general constraints\n",
      "Variable types: 675 continuous, 0 integer (0 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 1e+02]\n",
      "Presolve removed 750 rows and 580 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 63 rows, 95 columns, 201 nonzeros\n",
      "Variable types: 81 continuous, 14 integer (14 binary)\n",
      "Found heuristic solution: objective 40.0000000\n",
      "\n",
      "Root relaxation: objective 1.368348e+01, 62 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   13.68348    0    9   40.00000   13.68348  65.8%     -    0s\n",
      "     0     0   33.00000    0    6   40.00000   33.00000  17.5%     -    0s\n",
      "     0     0   33.00000    0    6   40.00000   33.00000  17.5%     -    0s\n",
      "     0     0   37.03846    0    2   40.00000   37.03846  7.40%     -    0s\n",
      "     0     0   39.00000    0    1   40.00000   39.00000  2.50%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 1\n",
      "  Implied bound: 3\n",
      "  MIR: 8\n",
      "  Flow cover: 14\n",
      "  Relax-and-lift: 2\n",
      "\n",
      "Explored 1 nodes (113 simplex iterations) in 0.05 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 40 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 4.000000000000e+01, best bound 4.000000000000e+01, gap 0.0000%\n",
      "Total Loss : 40.0\n",
      "Total number of bikes added in each hour:\n",
      " time\n",
      "7    25\n",
      "8    25\n",
      "9    25\n",
      "Name: bikes_added, dtype: int64\n",
      "Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 813 rows, 675 columns and 1620 nonzeros\n",
      "Model fingerprint: 0x2f4e5b06\n",
      "Model has 90 general constraints\n",
      "Variable types: 675 continuous, 0 integer (0 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 1e+02]\n",
      "Presolve removed 733 rows and 556 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 80 rows, 119 columns, 256 nonzeros\n",
      "Variable types: 101 continuous, 18 integer (18 binary)\n",
      "Found heuristic solution: objective 93.0000000\n",
      "\n",
      "Root relaxation: objective 3.704177e+01, 83 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   37.04177    0   17   93.00000   37.04177  60.2%     -    0s\n",
      "     0     0   89.15397    0    7   93.00000   89.15397  4.14%     -    0s\n",
      "     0     0     cutoff    0        93.00000   93.00000  0.00%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 2\n",
      "  Implied bound: 1\n",
      "  MIR: 13\n",
      "  Flow cover: 27\n",
      "  Relax-and-lift: 1\n",
      "\n",
      "Explored 1 nodes (138 simplex iterations) in 0.04 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 93 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 9.300000000000e+01, best bound 9.300000000000e+01, gap 0.0000%\n",
      "Total Loss : 93.0\n",
      "Total number of bikes added in each hour:\n",
      " time\n",
      "7    25.0\n",
      "8    25.0\n",
      "9    25.0\n",
      "Name: bikes_added, dtype: float64\n",
      "Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 813 rows, 675 columns and 1620 nonzeros\n",
      "Model fingerprint: 0x26913082\n",
      "Model has 90 general constraints\n",
      "Variable types: 675 continuous, 0 integer (0 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 1e+02]\n",
      "Presolve removed 743 rows and 560 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 70 rows, 115 columns, 236 nonzeros\n",
      "Variable types: 101 continuous, 14 integer (14 binary)\n",
      "Found heuristic solution: objective 83.0000000\n",
      "\n",
      "Root relaxation: objective 3.075054e+01, 70 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   30.75054    0   14   83.00000   30.75054  63.0%     -    0s\n",
      "     0     0   82.50000    0    3   83.00000   82.50000  0.60%     -    0s\n",
      "     0     0   82.50000    0    2   83.00000   82.50000  0.60%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 1\n",
      "  MIR: 6\n",
      "  Flow cover: 26\n",
      "\n",
      "Explored 1 nodes (96 simplex iterations) in 0.03 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 83 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 8.300000000000e+01, best bound 8.300000000000e+01, gap 0.0000%\n",
      "Total Loss : 83.0\n",
      "Total number of bikes added in each hour:\n",
      " time\n",
      "7    25\n",
      "8    25\n",
      "9    25\n",
      "Name: bikes_added, dtype: int64\n",
      "Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 813 rows, 675 columns and 1620 nonzeros\n",
      "Model fingerprint: 0xab6a7a4f\n",
      "Model has 90 general constraints\n",
      "Variable types: 675 continuous, 0 integer (0 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 1e+02]\n",
      "Presolve removed 751 rows and 583 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 62 rows, 92 columns, 196 nonzeros\n",
      "Variable types: 78 continuous, 14 integer (14 binary)\n",
      "Found heuristic solution: objective 45.0000000\n",
      "\n",
      "Root relaxation: objective 1.136000e+01, 56 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   11.36000    0   11   45.00000   11.36000  74.8%     -    0s\n",
      "     0     0   36.72214    0    5   45.00000   36.72214  18.4%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 2\n",
      "  Implied bound: 3\n",
      "  MIR: 7\n",
      "  Flow cover: 20\n",
      "  Relax-and-lift: 1\n",
      "\n",
      "Explored 1 nodes (89 simplex iterations) in 0.03 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 45 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 4.500000000000e+01, best bound 4.500000000000e+01, gap 0.0000%\n",
      "Total Loss : 45.0\n",
      "Total number of bikes added in each hour:\n",
      " time\n",
      "7    25\n",
      "8    25\n",
      "9    25\n",
      "Name: bikes_added, dtype: int64\n",
      "Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 813 rows, 675 columns and 1620 nonzeros\n",
      "Model fingerprint: 0xea4cbb22\n",
      "Model has 90 general constraints\n",
      "Variable types: 675 continuous, 0 integer (0 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+00, 1e+02]\n",
      "Presolve removed 732 rows and 554 columns\n",
      "Presolve time: 0.01s\n",
      "Presolved: 81 rows, 121 columns, 260 nonzeros\n",
      "Variable types: 103 continuous, 18 integer (18 binary)\n",
      "Found heuristic solution: objective 34.0000000\n",
      "\n",
      "Root relaxation: objective 3.000000e+00, 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    3.00000    0    9   34.00000    3.00000  91.2%     -    0s\n",
      "     0     0   19.58760    0   11   34.00000   19.58760  42.4%     -    0s\n",
      "     0     0     cutoff    0        34.00000   34.00000  0.00%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Implied bound: 9\n",
      "  MIR: 7\n",
      "  Flow cover: 28\n",
      "\n",
      "Explored 1 nodes (183 simplex iterations) in 0.04 seconds (0.00 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 1: 34 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 3.400000000000e+01, best bound 3.400000000000e+01, gap 0.0000%\n",
      "Total Loss : 33.99999999999997\n",
      "Total number of bikes added in each hour:\n",
      " time\n",
      "7    25.0\n",
      "8    25.0\n",
      "9    25.0\n",
      "Name: bikes_added, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "# run the model daily\n",
    "all_outputs = []\n",
    "g = morning_flow.set_index(['station', 'time']).groupby('date')\n",
    "for x in g.groups:\n",
    "    df = g.get_group(x)\n",
    "    odf = bike_rebalancing(df, num_bikes)\n",
    "    all_outputs.append(odf)\n",
    "output_df = pd.concat(all_outputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "dda9a673",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:46.965353Z",
     "start_time": "2022-11-18T00:21:46.939023Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>station</th>\n",
       "      <th>time</th>\n",
       "      <th>datetime</th>\n",
       "      <th>start_forecast</th>\n",
       "      <th>end_forecast</th>\n",
       "      <th>date</th>\n",
       "      <th>capacity</th>\n",
       "      <th>bikes_added</th>\n",
       "      <th>bikes_removed</th>\n",
       "      <th>loss_sale</th>\n",
       "      <th>beginning_inventory</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>11.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>14.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>17.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>17.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>22.0</td>\n",
       "      <td>41.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0.0</td>\n",
       "      <td>19</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>28.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>9.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>7.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>18.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>18.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>11.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>10 Ave &amp; W 14 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>55</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            station  time            datetime  start_forecast  end_forecast  \\\n",
       "0   1 Ave & E 62 St     7 2022-08-01 07:00:00            11.0          10.0   \n",
       "1   1 Ave & E 62 St     8 2022-08-01 08:00:00            14.0          14.0   \n",
       "2   1 Ave & E 62 St     9 2022-08-01 09:00:00            17.0          20.0   \n",
       "3   1 Ave & E 68 St     7 2022-08-01 07:00:00            17.0          32.0   \n",
       "4   1 Ave & E 68 St     8 2022-08-01 08:00:00            22.0          41.0   \n",
       "5   1 Ave & E 68 St     9 2022-08-01 09:00:00            28.0          48.0   \n",
       "6   1 Ave & E 78 St     7 2022-08-01 07:00:00             9.0           7.0   \n",
       "7   1 Ave & E 78 St     8 2022-08-01 08:00:00            18.0           5.0   \n",
       "8   1 Ave & E 78 St     9 2022-08-01 09:00:00            18.0           7.0   \n",
       "9  10 Ave & W 14 St     7 2022-08-01 07:00:00             3.0           5.0   \n",
       "\n",
       "         date  capacity  bikes_added  bikes_removed  loss_sale  \\\n",
       "0  2022-08-01        45          0.0              0        1.0   \n",
       "1  2022-08-01        45          0.0              0        0.0   \n",
       "2  2022-08-01        45          0.0              0        0.0   \n",
       "3  2022-08-01        62          0.0             15        0.0   \n",
       "4  2022-08-01        62          0.0             19        0.0   \n",
       "5  2022-08-01        62          0.0              0        0.0   \n",
       "6  2022-08-01        57          7.0              0        0.0   \n",
       "7  2022-08-01        57          8.0              0        0.0   \n",
       "8  2022-08-01        57         11.0              0        0.0   \n",
       "9  2022-08-01        55          0.0              2        0.0   \n",
       "\n",
       "   beginning_inventory  \n",
       "0                  0.0  \n",
       "1                  0.0  \n",
       "2                  0.0  \n",
       "3                  0.0  \n",
       "4                  0.0  \n",
       "5                  0.0  \n",
       "6                  0.0  \n",
       "7                  5.0  \n",
       "8                  0.0  \n",
       "9                  0.0  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output_df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3316684",
   "metadata": {},
   "source": [
    "# Model Enhancement"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1863587a",
   "metadata": {},
   "source": [
    "Looking at the `output_df`, you may notice stations that see a transfer of bikes (addition or removal) at every hour. Even worse are those stations where some bikes are removed from them in one hour but then bikes are added to them in the next hour. We know these additions/removals consume time and money. What can we do to avoid this situation?\n",
    "\n",
    "One way to formulate this is to introduce a fixed cost for the use of the truck (or the bicycle trailer) that transfers the bikes and then add this term to the objective function. With a cost associated with the transfer, the model is incentivized to use fewer number of transfers.\n",
    "\n",
    "So first, we should calculate the number of times a transfer has occurred. Any time bikes are added to a station or removed from a station, a transfer has happened. So, we need a way to link addition of bikes and removal of bikes to a transfer. One way to achieve this is to introduce two new binary variables as follows:\n",
    "\n",
    "- $x_{i,t}\\quad$: 1 if any bike is added to station $i$ at hour $t$; 0 otherwise\n",
    "- $w_{i,t}\\quad$: 1 if any bike is removed from station $i$ at hour $t$; 0 otherwise\n",
    "\n",
    "Next, we need to establish the relationship between $y_{i,t}$ with $x_{i,t}$, and $z_{i,t}$ with $w_{i,t}$. Basically, we want to say:\n",
    "if $y_{i,t} \\ge 0$, then $x_{i,t} = 1$ and if $z_{i,t} \\ge 0$, then $w_{i,t} = 1$. \n",
    "\n",
    "We introduce the following two constraints:\n",
    "\n",
    "\\begin{align}\n",
    "&y_{i,t} \\le M x_{i,t} &\\quad \\forall i \\in I, t \\in T \\tag{11}\\\\\n",
    "&z_{i,t} \\le M w_{i,t} &\\quad \\forall i \\in I, t \\in T \\tag{12}\\\\\n",
    "\\end{align}\n",
    "\n",
    "where $M$ is a large number.\n",
    "\n",
    "Constraint 11 ensures that if $y_{i,t}\\ge 0$, then $x_{i,t} = 1$. But by itself, this constraint cannot make $x_{i,t} =0$ if $y_{i,t} \\le 0$. The same is true with constraint 12. It ensures that $z_{i,t} \\ge 0$ make $w_{i,t} = 1$. However, it cannot force $w_{i,t} = 0$ if $z_{i,t} \\ge 0$.\n",
    "\n",
    "This can be achieved by the objective function.\n",
    "\n",
    "Total number of transfers is equal to sum of $x_{i,t}$ and $w_{i,t}$ and our goal is to minimize number of transfers. So, we add these terms to the objective function. In other words, our new objective function is:\n",
    "\n",
    "$$\\min \\sum_{i,t} (l_{i,t} + x_{i,t} + w_{i,t})$$\n",
    "\n",
    "Since minimizing total transfers is desired, the objective function tries to make both $x_{i,t}$ and $w_{i,t}$ as small as possible (or zero, in this case). Along with constraints 11 and 12, this means that for cases where $x_{i,t}$ and $w_{i,t}$ can take either 0 or 1, objective function forces them to get a value of 0. Moreover, since any extra transfer causes either $x_{i,t}$ or $w_{i,t}$ to be 1, the model is incentivized to have fewer transfers in order to minimize the objective function.\n",
    "\n",
    "Of course, you can make this model even more generic by having a coefficient for each term in the objective function (think of them as cost)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a8c272b6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:54.361099Z",
     "start_time": "2022-11-18T00:21:54.036805Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 1083 rows, 945 columns and 2160 nonzeros\n",
      "Model fingerprint: 0x361303a6\n",
      "Model has 90 general constraints\n",
      "Variable types: 675 continuous, 270 integer (270 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+03]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 1e+02]\n",
      "\n",
      "Warning: Completing partial solution with 270 unfixed non-continuous variables out of 270\n",
      "MIP start from previous solve produced solution with objective 52 (0.01s)\n",
      "MIP start from previous solve produced solution with objective 49 (0.03s)\n",
      "Loaded MIP start from previous solve with objective 49\n",
      "\n",
      "Presolve removed 974 rows and 789 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 109 rows, 156 columns, 320 nonzeros\n",
      "Variable types: 99 continuous, 57 integer (57 binary)\n",
      "\n",
      "Root relaxation: objective 1.700045e+01, 114 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   17.00045    0   34   49.00000   17.00045  65.3%     -    0s\n",
      "     0     0   42.89229    0   14   49.00000   42.89229  12.5%     -    0s\n",
      "     0     0   42.96598    0   16   49.00000   42.96598  12.3%     -    0s\n",
      "     0     0   46.50771    0   17   49.00000   46.50771  5.09%     -    0s\n",
      "     0     0   47.78795    0   11   49.00000   47.78795  2.47%     -    0s\n",
      "     0     0   47.78795    0   11   49.00000   47.78795  2.47%     -    0s\n",
      "     0     0   47.78795    0   17   49.00000   47.78795  2.47%     -    0s\n",
      "     0     0   47.78795    0   17   49.00000   47.78795  2.47%     -    0s\n",
      "     0     0   47.85591    0   14   49.00000   47.85591  2.33%     -    0s\n",
      "     0     0   48.02222    0    9   49.00000   48.02222  2.00%     -    0s\n",
      "     0     0   48.03684    0   13   49.00000   48.03684  1.97%     -    0s\n",
      "     0     0   48.06443    0   14   49.00000   48.06443  1.91%     -    0s\n",
      "     0     0   48.06443    0   13   49.00000   48.06443  1.91%     -    0s\n",
      "     0     0   48.08832    0   30   49.00000   48.08832  1.86%     -    0s\n",
      "     0     0   48.08832    0   10   49.00000   48.08832  1.86%     -    0s\n",
      "     0     0   48.08832    0   12   49.00000   48.08832  1.86%     -    0s\n",
      "     0     0   48.09318    0   12   49.00000   48.09318  1.85%     -    0s\n",
      "     0     0   48.09428    0   13   49.00000   48.09428  1.85%     -    0s\n",
      "     0     0   48.12906    0   15   49.00000   48.12906  1.78%     -    0s\n",
      "     0     0   48.13326    0   12   49.00000   48.13326  1.77%     -    0s\n",
      "     0     0   48.14788    0   13   49.00000   48.14788  1.74%     -    0s\n",
      "     0     0   48.15075    0   14   49.00000   48.15075  1.73%     -    0s\n",
      "     0     0   48.15143    0   14   49.00000   48.15143  1.73%     -    0s\n",
      "     0     0   48.15143    0   14   49.00000   48.15143  1.73%     -    0s\n",
      "     0     2   48.15143    0   14   49.00000   48.15143  1.73%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 10\n",
      "  Cover: 1\n",
      "  Implied bound: 8\n",
      "  MIR: 17\n",
      "  Flow cover: 31\n",
      "\n",
      "Explored 145 nodes (1500 simplex iterations) in 0.30 seconds (0.07 work units)\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 2: 49 52 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 4.900000000000e+01, best bound 4.900000000000e+01, gap 0.0000%\n"
     ]
    }
   ],
   "source": [
    "x = mdl.addVars(station_time, vtype=GRB.BINARY, name='x')  # 1 if y_{i,t} >= 0 \n",
    "w = mdl.addVars(station_time, vtype=GRB.BINARY, name='w')  # 1 if z_{i,t} >= 0\n",
    "\n",
    "big_m = 1000  # large number\n",
    "# relation between y and x\n",
    "mdl.addConstrs((y[i, t] <= big_m * x[i, t] for i, t in station_time), 'rel_y_x')\n",
    "# relation between z and w\n",
    "mdl.addConstrs((z[i, t] <= big_m * w[i, t] for i, t in station_time), 'rel_z_w')\n",
    "\n",
    "# new objective\n",
    "obj = l.sum() + (x.sum() + w.sum())\n",
    "mdl.setObjective(obj, GRB.MINIMIZE)\n",
    "mdl.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98a94077",
   "metadata": {},
   "source": [
    "## Post Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "bd6baf18",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-18T00:21:59.861287Z",
     "start_time": "2022-11-18T00:21:59.777976Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total Loss : 49.0\n",
      "Total number of bikes added in each hour:\n",
      " time\n",
      "7    25\n",
      "8    25\n",
      "9    25\n",
      "Name: bikes_added, dtype: int64\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>station</th>\n",
       "      <th>time</th>\n",
       "      <th>datetime</th>\n",
       "      <th>start_forecast</th>\n",
       "      <th>end_forecast</th>\n",
       "      <th>date</th>\n",
       "      <th>capacity</th>\n",
       "      <th>bikes_added</th>\n",
       "      <th>bikes_removed</th>\n",
       "      <th>loss_sale</th>\n",
       "      <th>beginning_inventory</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>11.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>14.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1 Ave &amp; E 62 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>17.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>45</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>17.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>22.0</td>\n",
       "      <td>41.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1 Ave &amp; E 68 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>28.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>62</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>9.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>8</td>\n",
       "      <td>2022-08-01 08:00:00</td>\n",
       "      <td>18.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>18</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1 Ave &amp; E 78 St</td>\n",
       "      <td>9</td>\n",
       "      <td>2022-08-01 09:00:00</td>\n",
       "      <td>18.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>57</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>10 Ave &amp; W 14 St</td>\n",
       "      <td>7</td>\n",
       "      <td>2022-08-01 07:00:00</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2022-08-01</td>\n",
       "      <td>55</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            station  time            datetime  start_forecast  end_forecast  \\\n",
       "0   1 Ave & E 62 St     7 2022-08-01 07:00:00            11.0          10.0   \n",
       "1   1 Ave & E 62 St     8 2022-08-01 08:00:00            14.0          14.0   \n",
       "2   1 Ave & E 62 St     9 2022-08-01 09:00:00            17.0          20.0   \n",
       "3   1 Ave & E 68 St     7 2022-08-01 07:00:00            17.0          32.0   \n",
       "4   1 Ave & E 68 St     8 2022-08-01 08:00:00            22.0          41.0   \n",
       "5   1 Ave & E 68 St     9 2022-08-01 09:00:00            28.0          48.0   \n",
       "6   1 Ave & E 78 St     7 2022-08-01 07:00:00             9.0           7.0   \n",
       "7   1 Ave & E 78 St     8 2022-08-01 08:00:00            18.0           5.0   \n",
       "8   1 Ave & E 78 St     9 2022-08-01 09:00:00            18.0           7.0   \n",
       "9  10 Ave & W 14 St     7 2022-08-01 07:00:00             3.0           5.0   \n",
       "\n",
       "         date  capacity  bikes_added  bikes_removed  loss_sale  \\\n",
       "0  2022-08-01        45            0              0        1.0   \n",
       "1  2022-08-01        45            0              0        0.0   \n",
       "2  2022-08-01        45            0              0        0.0   \n",
       "3  2022-08-01        62            0              0        0.0   \n",
       "4  2022-08-01        62            0              0        0.0   \n",
       "5  2022-08-01        62            0              0        0.0   \n",
       "6  2022-08-01        57            0              0        2.0   \n",
       "7  2022-08-01        57           18              0        0.0   \n",
       "8  2022-08-01        57            6              0        0.0   \n",
       "9  2022-08-01        55            0              0        0.0   \n",
       "\n",
       "   beginning_inventory  \n",
       "0                    0  \n",
       "1                    0  \n",
       "2                    0  \n",
       "3                    0  \n",
       "4                   15  \n",
       "5                   34  \n",
       "6                    0  \n",
       "7                    0  \n",
       "8                    5  \n",
       "9                    0  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame()\n",
    "if mdl.status == GRB.Status.OPTIMAL:\n",
    "    df = flow_df.copy()\n",
    "    df = df.merge(stations[['capacity']], left_on='station', right_index=True)\n",
    "    df[['bikes_added', 'bikes_removed', 'loss_sale', 'beginning_inventory']] = 0\n",
    "    for k, v in y.items():\n",
    "        df.loc[k, 'bikes_added'] = v.x\n",
    "        df.loc[k, 'bikes_removed'] = z[k].x\n",
    "        df.loc[k, 'beginning_inventory'] = q[k].x\n",
    "        df.loc[k, 'loss_sale'] = l[k].x\n",
    "    df.reset_index(inplace=True)\n",
    "    print(f'Total Loss : {mdl.objVal}')\n",
    "    total_bikes_added = df.groupby('time')['bikes_added'].sum()\n",
    "    print(f'Total number of bikes added in each hour:\\n {total_bikes_added}')\n",
    "else:\n",
    "    print('Could not find a feasible solution!')\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cf3c78f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-10-19T18:40:17.606647Z",
     "start_time": "2022-10-19T18:40:17.596562Z"
    }
   },
   "source": [
    "# Extra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6b18de9",
   "metadata": {},
   "source": [
    "## Scenario Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af26cbd6",
   "metadata": {},
   "source": [
    "An important requirements in many MO problems is scenario analysis or what-if analysis. \n",
    "Generally speaking, in what-if analysis, we're interested in knowing how the solution changes under various scenarios. \n",
    "Think about our case here.\n",
    "- How does the solution change if the number of reserved bikes increase or decrease by 10%?\n",
    "- In the enhanced model, what happens if the cost of lost sale change? How about the cost of visiting a location for adding or removing the bikes?\n",
    "- What if a new station is added close to the busiest station?\n",
    "- What if we want to ensure that every station have at least 2 available bikes at the beginning of each hour?\n",
    "\n",
    "Our model here is still a simple model. But you can imagine the value that the scenario analysis can provide. It can enable you to answer many questions by creating and comparing different scenarios and evaluating their outcomes, so that you can assess their impacts on the business goals. To learn more, you can check this [multi-scenario example](https://www.gurobi.com/documentation/9.5/examples/multiscenario_py.html) from Gurobi."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57187c2f",
   "metadata": {},
   "source": [
    "## How this problem is solved in reality?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40d32625",
   "metadata": {},
   "source": [
    "After knowing how many bikes needed in each station, the bikes need to be physically moved from one station to another. \n",
    "\n",
    "During rush hours, the traffic is already heavy. So, bicycle trailers (that can usually hold 5 bicycles) are used to move bikes around.\n",
    "\n",
    "During lighter hours (mainly at night), the bikes are transferred using trucks.\n",
    "\n",
    "In either case, some bikes should be removed from stations where there are more in-flow of bikes and should be transferred to stations where there are more out-flow of bikes to balance them out. \n",
    "This problem, where trucks need to go from one station to another and either pick up bikes or deliver them, is itself another mathematical optimization problem. \n",
    "\n",
    "In this problem, we need to ensure that all the stations that have a pickup or delivery, are visited during a certain time window and the goal can be to do this with minimum number of trucks or minimum transportation cost (for example, fuel cost plus the cost for using the truck). This problem is a variation of the famous Vehicle Routing Problem (VRP) with pickup and delivery.\n",
    "To learn more, check out [this webinar](https://www.gurobi.com/resource/how-to-synchronize-complex-routing-operations-synched-vrps-with-gurobi/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c83670ce",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}