p-median.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a136407c",
   "metadata": {},
   "source": [
    "$$\n",
    "z_{ij} = \n",
    "\\begin{Bmatrix}\n",
    "1, \\text{if the ith facility is assigned to the jth center} \\\\\n",
    "0, \\text{otherwise}\n",
    "\\end{Bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66d9156d",
   "metadata": {},
   "source": [
    "$$\n",
    "y_j \n",
    "= \n",
    "\\begin{Bmatrix}\n",
    "1, \\text{if jth location is selected as the location for a depot} \\\\\n",
    "0, \\text{otherwise}\n",
    "\\end{Bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d230caf4",
   "metadata": {},
   "source": [
    "$$\n",
    "d_{ij} = \\text{Distance between ith and jth location}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79fa518c",
   "metadata": {},
   "source": [
    "$$\n",
    "a_{i} = \\text{Weight (demand) of the ith location}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f73ac33",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\min z = & \\sum_{i=1}^{n} \\sum_{j=1}^n a_i d_{ij} z_{ij} \\\\\n",
    "\\text{Subject to:} & \\\\\n",
    "& \\sum_{j = 1}^{n} z_{ij = 1}  \\;\\;\\;, \\forall_i, \\;\\;\\; i, j = 1,2,\\ldots, n \\\\\n",
    "& z_{ij} \\le y_j  \\;\\;\\;, \\forall_{i, j}\\;\\;\\; i, j = 1,2,\\ldots, n \\\\\n",
    "& \\sum_{j=1}^n y_j = \\text{number_of_depots} \\\\\n",
    "& z_{ij}, y_j \\in \\{0,1\\} \\;\\;\\; i, j = 1,2,\\ldots, n \\\\\n",
    "& a_{ij} \\ge 0 \\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "55a59845",
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuMP, HiGHS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "843eae19",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "number_of_depots = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f0b94343",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×2 Matrix{Float64}:\n",
       " 1.0  5.0\n",
       " 2.0  4.0\n",
       " 4.0  3.0\n",
       " 2.0  1.0\n",
       " 7.0  8.0\n",
       " 8.0  7.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coords = Float64[\n",
    "    1 5;\n",
    "    2 4;\n",
    "    4 3;\n",
    "    2 1;\n",
    "    7 8;\n",
    "    8 7\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9c87f422",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "euclidean (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function euclidean(u, v)\n",
    "    return ((u .- v) .* (u .- v)) |> sum |> sqrt\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f7182d72",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 2)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n, p = size(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3da2d1ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×6 Matrix{Float64}:\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances = zeros(Float64, n, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bb0bf35e",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in 1:n\n",
    "    for j in i:n\n",
    "        d = euclidean(coords[i, :], coords[j, :])\n",
    "        distances[i, j] = d\n",
    "        distances[j, i] = d\n",
    "    end \n",
    "end "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e4c53f77",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×6 Matrix{Float64}:\n",
       " 0.0      1.41421  3.60555  4.12311  6.7082   7.28011\n",
       " 1.41421  0.0      2.23607  3.0      6.40312  6.7082\n",
       " 3.60555  2.23607  0.0      2.82843  5.83095  5.65685\n",
       " 4.12311  3.0      2.82843  0.0      8.60233  8.48528\n",
       " 6.7082   6.40312  5.83095  8.60233  0.0      1.41421\n",
       " 7.28011  6.7082   5.65685  8.48528  1.41421  0.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "7bc95e7c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running HiGHS 1.5.1 [date: 1970-01-01, git hash: 93f1876e4]\n",
      "Copyright (c) 2023 HiGHS under MIT licence terms\n",
      "Presolving model\n",
      "43 rows, 42 cols, 114 nonzeros\n",
      "43 rows, 42 cols, 114 nonzeros\n",
      "\n",
      "Solving MIP model with:\n",
      "   43 rows\n",
      "   42 cols (42 binary, 0 integer, 0 implied int., 0 continuous)\n",
      "   114 nonzeros\n",
      "\n",
      "        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      \n",
      "     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time\n",
      "\n",
      "         0       0         0   0.00%   0               inf                  inf        0      0      0         0     0.0s\n",
      " T       0       0         0   0.00%   0               8.064495102      100.00%        0      0      0        17     0.0s\n",
      "\n",
      "Solving report\n",
      "  Status            Optimal\n",
      "  Primal bound      8.06449510225\n",
      "  Dual bound        8.06449510225\n",
      "  Gap               0% (tolerance: 0.01%)\n",
      "  Solution status   feasible\n",
      "                    8.06449510225 (objective)\n",
      "                    0 (bound viol.)\n",
      "                    0 (int. viol.)\n",
      "                    0 (row viol.)\n",
      "  Timing            0.00 (total)\n",
      "                    0.00 (presolve)\n",
      "                    0.00 (postsolve)\n",
      "  Nodes             1\n",
      "  LP iterations     17 (total)\n",
      "                    0 (strong br.)\n",
      "                    0 (separation)\n",
      "                    0 (heuristics)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6-element Vector{Float64}:\n",
       " -0.0\n",
       "  1.0\n",
       "  0.0\n",
       " -0.0\n",
       "  0.0\n",
       "  1.0"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = Model(HiGHS.Optimizer)\n",
    "\n",
    "@variable(m, z[1:n, 1:n], Bin)\n",
    "@variable(m, y[1:n], Bin)\n",
    "\n",
    "@constraint(m, sum(y) == number_of_depots)\n",
    "\n",
    "for i in 1:n\n",
    "    for j in 1:n\n",
    "        @constraint(m, z[i, j] .<= y[j])\n",
    "    end \n",
    "end\n",
    "\n",
    "for i in 1:n\n",
    "    @constraint(m, sum(z[i, :]) == 1)\n",
    "end \n",
    "\n",
    "@objective(m, Min, sum(distances[1:n, 1:n] .* z[1:n, 1:n]))\n",
    "\n",
    "optimize!(m)\n",
    "\n",
    "value.(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "f76c3643",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×6 Matrix{Float64}:\n",
       " -0.0  1.0   0.0   0.0   0.0  0.0\n",
       " -0.0  1.0   0.0   0.0   0.0  0.0\n",
       "  0.0  1.0  -0.0  -0.0   0.0  0.0\n",
       "  0.0  1.0  -0.0  -0.0   0.0  0.0\n",
       "  0.0  0.0   0.0   0.0  -0.0  1.0\n",
       "  0.0  0.0   0.0   0.0  -0.0  1.0"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value.(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "aad18146",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Int64}:\n",
       " 2\n",
       " 6"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "location_of_depots = findall(x -> x == 1, value.(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97b8edaf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.0-rc1",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}