DABP_project_v1.5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gurobipy import *\n",
    "import numpy as np\n",
    "import csv\n",
    "import os\n",
    "import pandas as pd\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "file1 = \"ZipInfoCleaned.csv\"\n",
    "file2 = \"POD SitesNew.csv\"\n",
    "#file3 = \"Allegheny_population.csv\"\n",
    "file3 = \"HospitalSites.csv\"\n",
    "file4 = \"Scenarios.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "zipcode = pd.read_csv(file1,delimiter=\",\")\n",
    "podsites = pd.read_csv(file2,delimiter=\",\")\n",
    "#population = pd.read_csv(file3,delimiter=\",\")\n",
    "hospitals = pd.read_csv(file3,delimiter=\",\")\n",
    "scenarios = pd.read_csv(file4,delimiter=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "zipcode = zipcode[['ZIP', 'LAT', 'LON', 'Population']].dropna()\n",
    "#zipcode_np = zipcode[['LAT', 'LON']].to_numpy()\n",
    "#zipcode_np\n",
    "#merged_zipcode = population.merge(zipcode, left_on = \"ZIP\", right_on = \"ZIP\")\n",
    "zipcode_np = zipcode[['LAT', 'LON', 'Population']].to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "podsites1 = podsites[[\"SCHOOL/FACILITY NAME\", \"Latitude\", \"Longitude\"]].dropna()\n",
    "podsites_np = podsites1[[\"Latitude\",\"Longitude\"]].to_numpy()\n",
    "\n",
    "HospitalSites = hospitals[[\"Address for the hospital\", \"Latitude_hospital\", \"Longitude_hospital\"]].dropna()\n",
    "HospitalSites_np = HospitalSites[[\"Latitude_hospital\",\"Longitude_hospital\"]].to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenarios_merged = scenarios.merge(zipcode, left_on = \"ZIP\", right_on = \"ZIP\")\n",
    "scenarios_merged = scenarios_merged.drop([\"LAT\",\"LON\",\"Population\"], axis=1)\n",
    "scenarios_np = scenarios_merged.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(47, 92)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#calculate empty list of lists rows = podsites, cols = zips\n",
    "distance_zip = [[0 for i in range(len(zipcode_np))] for j in range(len(podsites_np))]\n",
    "np.shape(distance_zip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(47, 14)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#calculate empty list of lists rows =podsites, cols = hospitals\n",
    "distance_hos = [[0 for i in range(len(HospitalSites_np))] for j in range(len(podsites_np))]\n",
    "np.shape(distance_hos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Indices for POD sites and zipcodes\n",
    "i = len(podsites_np) # Potential POD sites\n",
    "j = len(zipcode_np) # All the zipcodes\n",
    "sites = range(i) # Equivalently, sites = [0,1,2,3,4]\n",
    "zips= range(j) # Equivalently, zips = [0,1,2,3]\n",
    "scenarios = scenarios_np\n",
    "scenarios = range(9)\n",
    "\n",
    "# calculate haversine distance\n",
    "for i in sites:\n",
    "    for j in zips:\n",
    "        #distance between sites (cols) and zips (rows)\n",
    "        \n",
    "        #Haversine Distance\n",
    "        #earth radius in meters\n",
    "        R = 6373000.\n",
    "\n",
    "        lat1 = np.radians(podsites_np[i,0])\n",
    "        lon1 = np.radians(podsites_np[i,1])\n",
    "        lat2 = np.radians(zipcode_np[j,0])\n",
    "        lon2 = np.radians(zipcode_np[j,1])\n",
    "\n",
    "        phi1 = np.radians(lat1)\n",
    "        phi2 = np.radians(lat2)\n",
    "        delta_phi = np.radians(lat2 - lat1)\n",
    "        delta_lambda = np.radians(lon2 - lon1)\n",
    "        a = np.sin(delta_phi / 2)**2 + np.cos(phi1) * np.cos(phi2) *   np.sin(delta_lambda / 2)**2\n",
    "        res = R * (2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)))\n",
    "        distance_zip[i][j] = np.round(res, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Indices for POD sites and Hospitals\n",
    "i = len(podsites_np) # Potential POD sites\n",
    "j = len(HospitalSites_np) # All the zipcodes\n",
    "sites = range(i) # Equivalently, sites = [0,1,2,3,4]\n",
    "hospitals= range(j) # Equivalently, hospital = [0,1,2,3]\n",
    "\n",
    "# calculate haversine distance\n",
    "for i in sites:\n",
    "    for j in hospitals:\n",
    "        #distance between sites (cols) and zips (rows)\n",
    "        \n",
    "        #Haversine Distance\n",
    "        #earth radius in meters\n",
    "        R = 6373000.\n",
    "\n",
    "        lat1 = np.radians(podsites_np[i,0])\n",
    "        lon1 = np.radians(podsites_np[i,1])\n",
    "        lat2 = np.radians(HospitalSites_np[j,0])\n",
    "        lon2 = np.radians(HospitalSites_np[j,1])\n",
    "\n",
    "        phi1 = np.radians(lat1)\n",
    "        phi2 = np.radians(lat2)\n",
    "        delta_phi = np.radians(lat2 - lat1)\n",
    "        delta_lambda = np.radians(lon2 - lon1)\n",
    "        a = np.sin(delta_phi / 2)**2 + np.cos(phi1) * np.cos(phi2) *   np.sin(delta_lambda / 2)**2\n",
    "        res = R * (2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)))\n",
    "        distance_hos[i][j] = np.round(res, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "DIST1 = np.array(distance_zip)\n",
    "DIST2 = np.array(distance_hos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = scenarios_np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = [[1/9],[1/9],[1/9],[1/9],[1/9],[1/9],[1/9],[1/9],[1/9]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "wage = 14.37\n",
    "cost1 = wage *8*260/20\n",
    "B = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([15017.  ,  4053.25,     0.  ,     0.  ,     0.  ,     0.  ,\n",
       "       16213.  ,  4053.25,     0.  , 12159.75])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi903\\gurobi.lic\n",
      "Academic license - for non-commercial use only\n"
     ]
    }
   ],
   "source": [
    "# Setting up model object\n",
    "m = Model(\"POD\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sites/zips/hospitals/scenarios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Decision Variable\n",
    "\n",
    "# 1st stage\n",
    "x = m.addVars(sites, vtype=GRB.BINARY)\n",
    "st = m.addVars(sites, lb=0.0, vtype=GRB.INTEGER)\n",
    "z = m.addVars(hospitals,sites, vtype=GRB.BINARY)\n",
    "\n",
    "# 2nd stage\n",
    "y = m.addVars(sites,zips,scenarios,vtype = GRB.BINARY)\n",
    "b = m.addVars(hospitals, sites, scenarios, lb=0.0, vtype=GRB.INTEGER)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Objective Function\n",
    "\n",
    "# Objective 1\n",
    "obj1_1 = LinExpr()\n",
    "obj1_2 = LinExpr()\n",
    "for s in scenarios:\n",
    "    for i in sites:\n",
    "        for j in zips:\n",
    "            obj1_1 += DIST1[i,j]*D[j,s]*y[i,j,s]*1/9\n",
    "\n",
    "for s in scenarios:\n",
    "    for i in sites:\n",
    "        for h in hospitals:\n",
    "            obj1_2 += DIST2[i,h] * b[h,i,s] *1/9\n",
    "m.setObjective(obj1_1+obj1_2, GRB.MINIMIZE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 100000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Constrains\n",
    "\n",
    "# Demand constraint: The total resources at a POD (initial stock+surge) should meet the demand of all neighborhoods that are assigned to it.\n",
    "for s in scenarios:\n",
    "    for i in sites:\n",
    "        m.addConstr((sum(b[h,i,s] for h in hospitals) + st[i]) >= sum(D[j,s] * y[i,j,s] for j in zips))\n",
    "        \n",
    "# Budget constraint\n",
    "m.addConstr(sum(st[i]*cost1 for i in sites) <= 0.4*B)\n",
    "\n",
    "            \n",
    "# One neighborhood (per zip code) can have (at most?) only one POD       \n",
    "for j in zips:\n",
    "    for s in scenarios:\n",
    "        m.addConstr(sum(y[i,j,s] for i in sites) == 1) ### DISCUSS <=/==\n",
    "        \n",
    "# Neighborhoods cannot be assigned to a POD that is not constructed.\n",
    "for i in sites:\n",
    "    for j in zips:\n",
    "        for s in scenarios:\n",
    "            m.addConstr(y[i,j,s] <= x[i]) \n",
    "\n",
    "# Only if a POD is built and it is assigned to a hospital, supplies can be transported from that hospital to that POD\n",
    "for i in sites:\n",
    "    for h in hospitals:\n",
    "        for s in scenarios:\n",
    "            m.addConstr(b[h,i,s] <= M*\n",
    "                        x[i]) \n",
    "for i in sites:\n",
    "    for h in hospitals:\n",
    "        for s in scenarios:\n",
    "            m.addConstr(b[h,i,s] <= M*z[h,i])\n",
    "            \n",
    "#Only if a POD is built there can be initial storage\n",
    "for i in sites:\n",
    "    m.addConstr(st[i] <= M*x[i])\n",
    "            \n",
    "        \n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.0.3 build v9.0.3rc0 (win64)\n",
      "Optimize a model with 52059 rows, 45590 columns and 166239 nonzeros\n",
      "Model fingerprint: 0x0669e203\n",
      "Variable types: 0 continuous, 45590 integer (39621 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+05]\n",
      "  Objective range  [2e+00, 5e+06]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Presolve removed 51920 rows and 40608 columns\n",
      "Presolve time: 0.28s\n",
      "Presolved: 139 rows, 4982 columns, 9306 nonzeros\n",
      "Variable types: 0 continuous, 4982 integer (4324 binary)\n",
      "Found heuristic solution: objective 8.863512e+07\n",
      "\n",
      "Root relaxation: objective 5.771705e+07, 217 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 5.7717e+07    0   12 8.8635e+07 5.7717e+07  34.9%     -    0s\n",
      "H    0     0                    5.786899e+07 5.7717e+07  0.26%     -    0s\n",
      "H    0     0                    5.786886e+07 5.7717e+07  0.26%     -    0s\n",
      "*    0     0               0    5.784777e+07 5.7848e+07  0.00%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 2\n",
      "  MIR: 3\n",
      "  RLT: 4\n",
      "\n",
      "Explored 1 nodes (229 simplex iterations) in 0.45 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 4: 5.78478e+07 5.78689e+07 5.7869e+07 8.86351e+07 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 5.784776892972e+07, best bound 5.784776892972e+07, gap 0.0000%\n"
     ]
    }
   ],
   "source": [
    "m.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "92\n",
      "92\n",
      "92\n",
      "92\n",
      "92\n",
      "92\n",
      "92\n",
      "92\n",
      "92\n",
      "47\n"
     ]
    }
   ],
   "source": [
    "#check if it worked\n",
    "totalzip_assigned = 0\n",
    "total_Sites =0\n",
    "for s in scenarios:\n",
    "    totalzip_assigned = 0\n",
    "    for i in sites:\n",
    "        for j in zips:\n",
    "\n",
    "            if y[i,j,s].x == 1:\n",
    "                totalzip_assigned +=1\n",
    "    print(totalzip_assigned)\n",
    "\n",
    "\n",
    "for i in sites:\n",
    "    if x[i].x == 1:\n",
    "        total_Sites+=1\n",
    "\n",
    "print(total_Sites)  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "zip_index =[]\n",
    "pod_index =[]\n",
    "pod_zip_pair =[]\n",
    "for i in sites:\n",
    "    for j in zips:\n",
    "        if y[i,j,2].x == 1:\n",
    "            pod_index.append(i)\n",
    "            zip_index.append(j)\n",
    "            pod_zip_pair.append((i,j))\n",
    "\n",
    "model_distance =[]            \n",
    "for pair in pod_zip_pair:\n",
    "    model_distance.append(distance_zip[pair[0]][pair[1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(model_distance)\n",
    "plt.title(\"Model distacne (in mtrs) Distribution\")\n",
    "plt.xlabel(\"distance\")\n",
    "plt.ylabel(\"Number of pairs\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = m.addVars(sites, vtype=GRB.BINARY)\n",
    "y = m.addVars(sites,zips, vtype = GRB.BINARY)\n",
    "st = m.addVar(lb=0.0, vtype=GRB.INTEGER)\n",
    "\n",
    "z = m.addVars(hospitals,sites, vtype=GRB.BINARY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.setObjective (sum(sum((distance[i][j]*y[i,j]*zipcode_np[j,2]) for i in sites) for j in zips))\n",
    "m.modelSense = GRB.MINIMIZE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#capacity constraint\n",
    "for i in sites:\n",
    "    m.addConstr(sum(zipcode_np[j,2]*y[i,j] for j in zips ) <= z[i]*x[i] )\n",
    "    \n",
    "# number of pod\n",
    "m.addConstr(sum(x[i] for i in sites) <= 10)# number of pod?\n",
    "\n",
    "#one zip area to one pod\n",
    "for j in zips:\n",
    "    m.addConstr(sum(y[i,j] for i in sites) == 1)\n",
    "for j in zips:\n",
    "    for i in sites:\n",
    "        m.addConstr(y[i,j] <= x[i])   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zip_index =[]\n",
    "pod_index =[]\n",
    "pod_zip_pair =[]\n",
    "for i in sites:\n",
    "    for j in zips:\n",
    "        if y[i,j].x == 1:\n",
    "            pod_index.append(i)\n",
    "            zip_index.append(j)\n",
    "            pod_zip_pair.append((i,j))\n",
    "\n",
    "model_distance =[]            \n",
    "for pair in pod_zip_pair:\n",
    "    model_distance.append(distance[pair[0]][pair[1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_distance\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(model_distance)\n",
    "plt.title(\"Model distacne (in mtrs) Distribution\")\n",
    "plt.xlabel(\"distance\")\n",
    "plt.ylabel(\"Numbe of pairs\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.7.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}