CSP-git-V1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyomo.environ import *\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import time\n",
    "import pandas as pd\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.061871</td>\n",
       "      <td>0.709058</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.985891</td>\n",
       "      <td>0.023800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.816979</td>\n",
       "      <td>0.094308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.301882</td>\n",
       "      <td>0.500483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.729545</td>\n",
       "      <td>0.863813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.848528</td>\n",
       "      <td>0.914105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.743942</td>\n",
       "      <td>0.028500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.583924</td>\n",
       "      <td>0.223909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.083004</td>\n",
       "      <td>0.596695</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.949340</td>\n",
       "      <td>0.051097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.784259</td>\n",
       "      <td>0.007328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0.030885</td>\n",
       "      <td>0.732157</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>0.190174</td>\n",
       "      <td>0.033792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.651377</td>\n",
       "      <td>0.469025</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>0.116406</td>\n",
       "      <td>0.815622</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>0.416353</td>\n",
       "      <td>0.716536</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>0.696187</td>\n",
       "      <td>0.473790</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.269826</td>\n",
       "      <td>0.977890</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>0.297694</td>\n",
       "      <td>0.464140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>0.603382</td>\n",
       "      <td>0.224822</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           x         y\n",
       "0   0.061871  0.709058\n",
       "1   0.985891  0.023800\n",
       "2   0.816979  0.094308\n",
       "3   0.301882  0.500483\n",
       "4   0.729545  0.863813\n",
       "5   0.848528  0.914105\n",
       "6   0.743942  0.028500\n",
       "7   0.583924  0.223909\n",
       "8   0.083004  0.596695\n",
       "9   0.949340  0.051097\n",
       "10  0.784259  0.007328\n",
       "11  0.030885  0.732157\n",
       "12  0.190174  0.033792\n",
       "13  0.651377  0.469025\n",
       "14  0.116406  0.815622\n",
       "15  0.416353  0.716536\n",
       "16  0.696187  0.473790\n",
       "17  0.269826  0.977890\n",
       "18  0.297694  0.464140\n",
       "19  0.603382  0.224822"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame()\n",
    "N = 20\n",
    "df['x']=[random.random() for i in range(N)]\n",
    "df['y']=[random.random() for i in range(N)]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAHSCAYAAADfUaMwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAafUlEQVR4nO3dUWxm513n8d+/zkQyLcIVKbOMEzZZFAyRSjVlCBWF3SmIdVK0ylCxUgui2goURUsRV1YzewEXXDRobhCiEEVVVCGtiFbLyIQlYCFVpmhLl7RyyTStDLOpaMfDbreAy06wlMn02YvxJJ5ZO35n5rx+7NefjzTSvOecjJ/+leSb8/j0uFprAQD6eVPvBQDAYSfGANCZGANAZ2IMAJ2JMQB0JsYA0Nkdvb7wXXfd1e69995tz7388st585vfvLcLmnBmOjwzHZ6ZDs9Mh3c7M/3c5z739dba22483i3G9957bz772c9ue255eTknT57c2wVNODMdnpkOz0yHZ6bDu52ZVtXfbnfcNjUAdCbGANCZGANAZ2IMAJ2JMQB0tmuMq+rpqvpaVX1hh/NVVb9ZVeer6oWqeufwywSAyTXKnfEnkjz0BucfTnL/5q9Hk/zO7S8LAA6PXWPcWvtUkn94g0seSfK77arPJJmpqu8caoEAMOmG+J7xbJKvbvl8YfMYADCCId7AVdsca9teWPVorm5l5+jRo1leXt72D7x06dKO57g1Zjo8Mx2emQ7PTIc3jpkOEeMLSe7Z8vnuJBe3u7C19lSSp5LkxIkTbafXiXl92/DMdHhmOjwzHZ6ZDm8cMx1im/rZJB/cfKr6XUm+0Vr7uwH+XAA4FHa9M66q30tyMsldVXUhya8mOZIkrbUnkzyX5L1Jzif55yQfGtdiAWAS7Rrj1toHdjnfkvziYCsCgEPGG7gAoLNuP8+Yyba4spYzS6u5uL6RYzPTWZify6nj/h9vANsRYwa3uLKW02fPZePylSTJ2vpGTp89lySCDLAN29QM7szS6mshvmbj8pWcWVrttCKA/U2MGdzF9Y2bOg5w2Ikxgzs2M31TxwEOOzFmcAvzc5k+MnXdsekjU1mYn+u0IoD9zQNcDO7aQ1qepgYYjRgzFqeOz4ovwIhsUwNAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnd3RewEAHAyLK2s5s7Sai+sbOTYznYX5uZw6Ptt7WRNBjAHY1eLKWk6fPZeNy1eSJGvrGzl99lySCPIAbFMDsKszS6uvhfiajctXcmZptdOKJosYA7Cri+sbN3WcmyPGAOzq2Mz0TR3n5ogxALtamJ/L9JGp645NH5nKwvxcpxVNFg9wAbCraw9peZp6PMQYgJGcOj4rvmNimxoAOhNjAOhMjAGgMzEGgM7EGAA6E2MA6EyMAaAzMQaAzsQYADoTYwDoTIwBoDMxBoDO/KAIgDFaXFnzk47YlRgDjMniylpOnz2XjctXkiRr6xs5ffZckggy17FNDTAmZ5ZWXwvxNRuXr+TM0mqnFbFfiTHAmFxc37ip4xxeYgwwJsdmpm/qOIeXGAOMycL8XKaPTF13bPrIVBbm5zqtiP3KA1wAY3LtIS1PU7MbMQYYo1PHZ8WXXdmmBoDOxBgAOhNjAOhMjAGgMzEGgM7EGAA6E2MA6EyMAaAzMQaAzsQYADoTYwDoTIwBoDM/KGIMFlfW/JQWAEYmxgNbXFnL6bPnsnH5SpJkbX0jp8+eSxJBBmBbtqkHdmZp9bUQX7Nx+UrOLK12WhEA+50YD+zi+sZNHQcAMR7YsZnpmzoOACPFuKoeqqrVqjpfVY9vc/7bquoPq+qvqurFqvrQ8Es9GBbm5zJ9ZOq6Y9NHprIwP9dpRQDsd7s+wFVVU0k+luQnklxI8nxVPdta++KWy34xyRdba/+uqt6WZLWq/nNr7ZWxrHofu/aQlqepARjVKE9TP5jkfGvtpSSpqmeSPJJka4xbkm+tqkryliT/kOTVgdd6YJw6Piu+AIxslG3q2SRf3fL5wuaxrX4ryfcluZjkXJJfbq19c5AVAsCEG+XOuLY51m74PJ/k80l+LMl3J/nTqvrz1to/XfcHVT2a5NEkOXr0aJaXl7f9gpcuXdrxHLfGTIdnpsMz0+GZ6fDGMdNRYnwhyT1bPt+dq3fAW30oyROttZbkfFV9Ocn3JvnLrRe11p5K8lSSnDhxop08eXLbL7i8vJydznFrzHR4Zjo8Mx2emQ5vHDMdZZv6+ST3V9V9VXVnkvcnefaGa76S5MeTpKqOJplL8tKQCwWASbXrnXFr7dWq+nCSpSRTSZ5urb1YVY9tnn8yya8l+URVncvVbe2PtNa+PsZ1A8DEGOnd1K2155I8d8OxJ7f8/mKSfzvs0gDgcPAGLgDoTIwBoDMxBoDOxBgAOhNjAOhMjAGgMzEGgM7EGAA6E2MA6EyMAaAzMQaAzsQYADoTYwDoTIwBoDMxBoDOxBgAOhNjAOhMjAGgMzEGgM7EGAA6u6P3Am7X4spaziyt5uL6Ro7NTGdhfi6njs/2XhYAjOxAx3hxZS2nz57LxuUrSZK19Y2cPnsuSQQZgAPjQG9Tn1lafS3E12xcvpIzS6udVgQAN+9Ax/ji+sZNHQeA/ehAx/jYzPRNHQeA/ehAx3hhfi7TR6auOzZ9ZCoL83OdVgQAN+9AP8B17SEtT1OzX3i6H7gVBzrGydUg+5cd+4Gn+4FbdaC3qWE/8XQ/cKsO/J3xrbKdyNA83Q/cqkN5Z3xtO3FtfSMtr28nLq6s9V4aB5in+4FbdShjbDuRcfB0P3CrDuU2te1ExsHT/cCtOpQxPjYznbVtwms7kdvl6X7gVhzKbWrbiQDsJ4fyzth2IgD7yaGMcWI7EYD941BuUwPAfiLGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdjRTjqnqoqlar6nxVPb7DNSer6vNV9WJV/dmwywSAyXXHbhdU1VSSjyX5iSQXkjxfVc+21r645ZqZJL+d5KHW2leq6jvGtF4AmDi7xjjJg0nOt9ZeSpKqeibJI0m+uOWan0lytrX2lSRprX1t6IX2sLiyljNLq7m4vpFjM9NZmJ/LqeOzvZcFwIQZZZt6NslXt3y+sHlsq+9J8taqWq6qz1XVB4daYC+LK2s5ffZc1tY30pKsrW/k9NlzWVxZ6700ACZMtdbe+IKqf59kvrX2C5uffy7Jg621X9pyzW8lOZHkx5NMJ/mLJD/ZWvvrG/6sR5M8miRHjx79gWeeeWbbr3np0qW85S1vudX/TYNY/V//N69c+eb/d/zOqTdl7l98a4cV3Z79MNNJY6bDM9Phmenwbmem73nPez7XWjtx4/FRtqkvJLlny+e7k1zc5pqvt9ZeTvJyVX0qyTuSXBfj1tpTSZ5KkhMnTrSTJ09u+wWXl5ez07m98qHH/yhtm42DSvLlJ07u+Xpu136Y6aQx0+GZ6fDMdHjjmOko29TPJ7m/qu6rqjuTvD/Jszdc8wdJfrSq7qiqb0nyQ0m+NOhK99ixmembOg4At2rXGLfWXk3y4SRLuRrY/9Jae7GqHquqxzav+VKSP0nyQpK/TPLx1toXxrfs8VuYn8v0kanrjk0fmcrC/FynFQEwqUbZpk5r7bkkz91w7MkbPp9Jcma4pfV17alpT1MDMG4jxfiwOnV8VnwBGDuvwwSAzsQYADoTYwDoTIwBoDMxBoDOxBgAOhNjAOhMjAGgMzEGgM7EGAA6E2MA6EyMAaAzMQaAzsQYADoTYwDoTIwBoDMxBoDOxBgAOhNjAOhMjAGgMzEGgM7EGAA6E2MA6EyMAaAzMQaAzsQYADoTYwDoTIwBoDMxBoDOxBgAOhNjAOjsjt4LgK0WV9ZyZmk1F9c3cmxmOgvzczl1fLb3sgDGSozZNxZX1nL67LlsXL6SJFlb38jps+eSRJCBiWabmn3jzNLqayG+ZuPylZxZWu20IoC9IcbsGxfXN27qOMCkEGP2jWMz0zd1HGBSiDH7xsL8XKaPTF13bPrIVBbm5zqtCGBveICLfePaQ1qepgYOGzFmXzl1fFZ8gUPHNjUAdCbGANCZGANAZ75nDDACr2plnMQYYBde1cq4iTEwsYa6m32jV7WKMUMQY2AiDXk361WtjJsHuNgXFlfW8u4nPpn7Hv+jvPuJT2ZxZa33kjjghvzBI17VyriJMd1du4NZW99Iy+t3MILM7RjybtarWhk3MaY7PzqRcRjybvbU8dl89H1vz+zMdCrJ7Mx0Pvq+t/t+MYPxPWO68/04xmFhfu667xknt3c361WtjJM7Y7rz/TjGwd0sB4k7Y7ob+g4GrnE3y0EhxnTnRycCh50Ysy+4gwEOM98zBoDOxBgAOhNjAOhMjAGgMzEGgM7EGAA6E2MA6EyMAaAzMQaAzsQYADoTYwDoTIwBoDMxBoDOxBgAOhNjAOhMjAGgMzEGgM5GinFVPVRVq1V1vqoef4PrfrCqrlTVTw+3RACYbLvGuKqmknwsycNJHkjygap6YIfrfj3J0tCLBIBJNsqd8YNJzrfWXmqtvZLkmSSPbHPdLyX5/SRfG3B9ADDxRonxbJKvbvl8YfPYa6pqNslPJXlyuKUBwOFwxwjX1DbH2g2ffyPJR1prV6q2u3zzD6p6NMmjSXL06NEsLy9ve92lS5d2PMetMdPhmenwzHR4Zjq8ccx0lBhfSHLPls93J7l4wzUnkjyzGeK7kry3ql5trS1uvai19lSSp5LkxIkT7eTJk9t+weXl5ex0jltjpsMz0+GZ6fDMdHjjmOkoMX4+yf1VdV+StSTvT/IzWy9ord137fdV9Ykk/+3GEAMA29s1xq21V6vqw7n6lPRUkqdbay9W1WOb532fGABuwyh3xmmtPZfkuRuObRvh1tp/uP1lAcDh4Q1cANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ3d0XsBANesb1zOu5/4ZC6ub+TYzHQW5udy6vhs72XB2IkxsCcWV9ZyZml1x9Aurqxl7R83srY+lSRZW9/I6bPnkkSQmXi2qYGxW1xZy+mz57K2vpGW10O7uLL22jVnllbzzdau++s2Ll/JmaXVPV4t7D0xBsbuzNJqNi5fue7YjaG9uL6x7V+703GYJGIMjN0ooT02M73tNTsdh0kixsDYjRLahfm5vKnquvPTR6ayMD831rXBfiDGwNgtzM9l+sjUdcduDO2p47OZfet0ZmemU0lmZ6bz0fe93cNbHAqepgbG7lpQ3+hp6iSZmT6S//74yQ4rhL7EGNgTp47PusuFHdimBoDORopxVT1UVatVdb6qHt/m/M9W1Qubvz5dVe8YfqkAMJl2jXFVTSX5WJKHkzyQ5ANV9cANl305yb9prX1/kl9L8tTQCwWASTXKnfGDSc631l5qrb2S5Jkkj2y9oLX26dbaP25+/EySu4ddJgBMrlFiPJvkq1s+X9g8tpOfT/LHt7MoADhMRnmaurY51rY5lqp6T67G+Ed2OP9okkeT5OjRo1leXt72C166dGnHc9waMx2emQ7PTIdnpsMbx0xHifGFJPds+Xx3kos3XlRV35/k40kebq39/XZ/UGvtqWx+P/nEiRPt5MmT237B5eXl7HSOW2OmwzPT4Znp8Mx0eOOY6Sjb1M8nub+q7quqO5O8P8mzWy+oqu9KcjbJz7XW/nrQFQLAhNv1zri19mpVfTjJUpKpJE+31l6sqsc2zz+Z5FeSfHuS366r75Z9tbV2YnzLBoDJMdIbuFprzyV57oZjT275/S8k+YVhlwYAh4M3cAFAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0BnYgwAnYkxAHQmxgDQmRgDQGdiDACdiTEAdHZH7wUAsPcWV9ZyZmk1F9c3cmxmOgvzczl1fLb3sg4tMQY4ZBZX1nL67LlsXL6SJFlb38jps+eSRJA7sU0NcMicWVp9LcTXbFy+kjNLq51WhBgDHDIX1zdu6jjjJ8YAh8yxmembOs74iTHAIbMwP5fpI1PXHZs+MpWF+blOK8IDXACHzLWHtDxNvX+IMcAhdOr4rPjuI7apAaAzMQaAzsQYADoTYwDoTIwBoDMxBoDOxBgAOhNjAOhMjAGgM2/gAmDiLK6sHajXfYoxABNlcWUtp8+ee+1nNq+tb+T02XNJsm+DbJsagIlyZmn1tRBfs3H5Ss4srXZa0e7EGICJcnF946aO7wdiDMBEOTYzfVPH9wMxBmCiLMzPZfrI1HXHpo9MZWF+rtOKducBLgAmyrWHtDxNTTcH7XF+gHE4dXz2QP27T4wnyBs9zj/TcV0AvDHfM54gB/FxfjisFlfW8u4nPpn7Hv+jvPuJT2ZxZa33kujInfEEeePH+d+8t4sBdnQQX0rBeLkzniAH8XF+OIzsYu1vPXYtxHiCHMTH+eEwOogvpTgsru1arK1vpOX1XYtxB1mMJ8ip47P56PventmZ6VSS2ZnpfPR9b7ftBfuMXaz9q9euhe8ZT5iD9jg/HEYL83PXfc84sYu1X/TatXBnDLDH7GLtX712LdwZA3RgF2t/6rVrIcYAsKnXqzTFGOAA8wrc4fXYtRBjgAPKy0Mmx0gPcFXVQ1W1WlXnq+rxbc5XVf3m5vkXquqdwy8VgK28PGRy7BrjqppK8rEkDyd5IMkHquqBGy57OMn9m78eTfI7A68TgBt4ecjkGOXO+MEk51trL7XWXknyTJJHbrjmkSS/2676TJKZqvrOgdcKwBZeHjI5RonxbJKvbvl8YfPYzV4DwIC8AndyjPIAV21zrN3CNamqR3N1GztHjx7N8vLytl/w0qVLO57j1pjp8Mx0eGZ6c2aSfPSHp/K/v3E5r1z5Zu6celOOftudmfnG32R5+W+SmOk4jGOmo8T4QpJ7tny+O8nFW7gmrbWnkjyVJCdOnGgnT57c9gsuLy9np3PcGjMdnpkOz0yHZ6bDG8dMR9mmfj7J/VV1X1XdmeT9SZ694Zpnk3xw86nqdyX5Rmvt7wZdKQBMqF3vjFtrr1bVh5MsJZlK8nRr7cWqemzz/JNJnkvy3iTnk/xzkg+Nb8kAMFlGeulHa+25XA3u1mNPbvl9S/KLwy4NAA4HP7UJADoTYwDoTIwBoDMxBoDOxBgAOhNjAOhMjAGgMzEGgM7EGAA6E2MA6EyMAaAzMQaAzurqz3jo8IWr/k+Sv93h9F1Jvr6HyzkMzHR4Zjo8Mx2emQ7vdmb6L1trb7vxYLcYv5Gq+mxr7UTvdUwSMx2emQ7PTIdnpsMbx0xtUwNAZ2IMAJ3t1xg/1XsBE8hMh2emwzPT4Znp8Aaf6b78njEAHCb79c4YAA6NrjGuqoeqarWqzlfV49ucr6r6zc3zL1TVO3us8yAZYaY/uznLF6rq01X1jh7rPEh2m+mW636wqq5U1U/v5foOolFmWlUnq+rzVfViVf3ZXq/xoBnhn/1vq6o/rKq/2pzph3qs86Coqqer6mtV9YUdzg/bp9Zal19JppL8zyT/KsmdSf4qyQM3XPPeJH+cpJK8K8n/6LXeg/BrxJn+cJK3bv7+YTO9/Zluue6TSZ5L8tO9172ff4349+lMki8m+a7Nz9/Re937+deIM/1PSX598/dvS/IPSe7svfb9+ivJv07yziRf2OH8oH3qeWf8YJLzrbWXWmuvJHkmySM3XPNIkt9tV30myUxVfedeL/QA2XWmrbVPt9b+cfPjZ5LcvcdrPGhG+fs0SX4pye8n+dpeLu6AGmWmP5PkbGvtK0nSWjPXNzbKTFuSb62qSvKWXI3xq3u7zIOjtfapXJ3RTgbtU88Yzyb56pbPFzaP3ew1vO5m5/Xzufpfduxs15lW1WySn0ry5B6u6yAb5e/T70ny1qparqrPVdUH92x1B9MoM/2tJN+X5GKSc0l+ubX2zb1Z3kQatE933PZybl1tc+zGR7tHuYbXjTyvqnpPrsb4R8a6ooNvlJn+RpKPtNauXL3pYBejzPSOJD+Q5MeTTCf5i6r6TGvtr8e9uANqlJnOJ/l8kh9L8t1J/rSq/ry19k9jXtukGrRPPWN8Ick9Wz7fnav/xXaz1/C6keZVVd+f5ONJHm6t/f0ere2gGmWmJ5I8sxniu5K8t6peba0t7skKD55R/9n/emvt5SQvV9WnkrwjiRhvb5SZfijJE+3qNzzPV9WXk3xvkr/cmyVOnEH71HOb+vkk91fVfVV1Z5L3J3n2hmueTfLBzafW3pXkG621v9vrhR4gu860qr4rydkkP+cuYyS7zrS1dl9r7d7W2r1J/muS/yjEb2iUf/b/IMmPVtUdVfUtSX4oyZf2eJ0HySgz/Uqu7jSkqo4mmUvy0p6ucrIM2qdud8attVer6sNJlnL1ScCnW2svVtVjm+efzNUnU9+b5HySf87V/7JjByPO9FeSfHuS3968k3u1eYn8jkacKTdhlJm21r5UVX+S5IUk30zy8dbatv8XE0b++/TXknyiqs7l6hbrR1prfprTDqrq95KcTHJXVV1I8qtJjiTj6ZM3cAFAZ97ABQCdiTEAdCbGANCZGANAZ2IMAJ2JMQB0JsYA0JkYA0Bn/w8GZd6SavuEXQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.scatter(df['x'],df['y'])\n",
    "plt.grid()\n",
    "plt.savefig(f\"base_nopath.png\", format='png', dpi=400)  # Save the figure as a PNG image file\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dist(n,m):\n",
    "    return np.sqrt( (df.loc[n-1,'x']-df.loc[m-1,'x'])**2+(df.loc[n-1,'y']-df.loc[m-1,'y'])**2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 12),\n",
       " (2, 10),\n",
       " (4, 19),\n",
       " (7, 11),\n",
       " (8, 20),\n",
       " (10, 2),\n",
       " (11, 7),\n",
       " (12, 1),\n",
       " (14, 17),\n",
       " (17, 14),\n",
       " (19, 4),\n",
       " (20, 8)]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes = [i for i in range(1,1+len(df))]\n",
    "R= 0.07\n",
    "dic={(n,m):dist(n,m) for n in nodes for m in nodes if dist(n,m)>0}\n",
    "cover = []\n",
    "for (n,m),d in dic.items():\n",
    "    if d<R and d>0:\n",
    "        cover+=[(n,m)]\n",
    "cover"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(df)\n",
    "model = AbstractModel()\n",
    "model.i = RangeSet(N)\n",
    "model.j = Set(initialize=model.i)\n",
    "\n",
    "model.U = Var(model.i,model.j, within = Binary)\n",
    "model.Tour = Var(model.i, within = Binary)\n",
    "model.M = Var(model.i,model.j,bounds=(0,N), within = NonNegativeReals)\n",
    "\n",
    "#model.Neighbor = Var(model.i, within = Binary)\n",
    "model.Order = Var(model.i, bounds=(0,N), within = NonNegativeReals)\n",
    "\n",
    "def rule_C0(model,i):\n",
    "    return  sum(model.U[j,i] for j in model.j if (j,i) in dic) == model.Tour[i] \n",
    "model.C0 = Constraint(model.i, rule=rule_C0)\n",
    "\n",
    "def rule_C1(model,i):\n",
    "    return  sum(model.U[i,j] for j in model.j if (i,j) in dic) == model.Tour[i] \n",
    "model.C1 = Constraint(model.i, rule=rule_C1)\n",
    "\n",
    "# both i,j should be on tour path \n",
    "def rule_C1A(model,i,j):\n",
    "    if (i,j) in dic:\n",
    "        return  model.U[i,j] <= model.Tour[i] \n",
    "    else:\n",
    "        return Constraint.Skip\n",
    "model.C1A = Constraint(model.i,model.j, rule=rule_C1A)\n",
    "def rule_C1B(model,i,j):\n",
    "    if (i,j) in dic:\n",
    "        return  model.U[i,j] <= model.Tour[j] \n",
    "    else:\n",
    "        return Constraint.Skip\n",
    "model.C1B = Constraint(model.i,model.j, rule=rule_C1B)\n",
    "\n",
    "def rule_C3(model,i):\n",
    "    return  1- model.Tour[i] <= sum(model.Tour[j] for j in model.j if (i,j) in cover) \n",
    "model.C3 = Constraint(model.i, rule=rule_C3)\n",
    "\n",
    "def rule_C4(model,i,j):\n",
    "    if j >1 and i>0 and (i,j) in dic:\n",
    "        return  model.Order[i]- model.Order[j]+ model.M[i,j] <= sum(model.Tour[n] for n in model.i)-1 +10*N*(2-model.Tour[i]-model.Tour[j])\n",
    "    else:\n",
    "        return Constraint.Skip\n",
    "model.C4 = Constraint(model.i,model.j, rule=rule_C4)\n",
    "\n",
    "def rule_C5(model,i,j):\n",
    "    if (i,j) in dic:\n",
    "        return  model.M[i,j] <= N*model.U[i,j]\n",
    "    else:\n",
    "        return Constraint.Skip\n",
    "model.C5 = Constraint(model.i,model.j, rule=rule_C5)\n",
    "\n",
    "def rule_C6(model,i,j):\n",
    "    if (i,j) in dic:\n",
    "        return  model.M[i,j]  <= sum(model.Tour[n] for n in model.i)\n",
    "    else:\n",
    "        return Constraint.Skip\n",
    "model.C6 = Constraint(model.i,model.j, rule=rule_C6)\n",
    "\n",
    "def rule_C7(model,i,j):\n",
    "    if (i,j) in dic:\n",
    "        return  model.M[i,j] >= sum(model.Tour[n] for n in model.i) - (1-model.U[i,j])*N\n",
    "    else:\n",
    "        return Constraint.Skip\n",
    "model.C7 = Constraint(model.i,model.j, rule=rule_C7)\n",
    "\n",
    "def rule_C8(model,i):\n",
    "    return  model.Order[i]  <= N*model.Tour[i]\n",
    "model.C8 = Constraint(model.i, rule=rule_C8)\n",
    "def rule_C8B(model,i):\n",
    "    return  model.Order[i]  >= model.Tour[i]\n",
    "model.C8B = Constraint(model.i, rule=rule_C8B)\n",
    "\n",
    "def rule_C9(model,i,j):\n",
    "    if (i,j) in dic and i>j:\n",
    "        return  model.U[i,j]+model.U[j,i] <=1\n",
    "    else:\n",
    "        return Constraint.Skip\n",
    "model.C9 = Constraint(model.i,model.j, rule=rule_C9)\n",
    "\n",
    "def rule_of(model):\n",
    "    return sum(model.U[i,j]*dic[i,j] for i in model.i for j in model.j if (i,j) in dic)\n",
    "model.obj = Objective(rule=rule_of, sense=minimize)\n",
    "instance = model.create_instance()\n",
    "\n",
    "instance.Order[1].fix(1)\n",
    "instance.Tour[1].fix(1)\n",
    "    \n",
    "opt = SolverFactory('cbc')\n",
    "opt.options = {'sec': 900, 'threads': 4, 'ratio': 0.05}\n",
    "results = opt.solve(instance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(value(instance.obj))\n",
    "\n",
    "def circ(x,y,r):\n",
    "    theta = np.linspace(-np.pi,np.pi,100)\n",
    "    X= [x+r*cos(th) for th in theta]\n",
    "    Y= [y+r*sin(th) for th in theta]\n",
    "    plt.plot(X,Y,'--k', alpha=0.3)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9,9))\n",
    "#plt.scatter(df['x'],df['y'],s= 100)\n",
    "plt.grid()\n",
    "for j in instance.j:\n",
    "    #print(j, value(instance.Order[j]) ,value(instance.Tour[j]) )\n",
    "    for i in instance.i:\n",
    "        if (i,j) in dic and value(instance.U[i,j])>0:        \n",
    "            x1,y1 = df.loc[i-1,'x'], df.loc[i-1,'y']\n",
    "            x2,y2 = df.loc[j-1,'x'], df.loc[j-1,'y']\n",
    "            plt.plot([x1,x2],[y1,y2])\n",
    "            plt.scatter(x1,y1,s=100,c='k')\n",
    "            circ(x1,y1,R)\n",
    "            print(i,j, value(instance.Order[i]) , value(instance.Order[j]))\n",
    "for i in instance.i:\n",
    "    if value(instance.Tour[i]) ==0:\n",
    "        x1,y1 = df.loc[i-1,'x'], df.loc[i-1,'y']\n",
    "        plt.scatter(x1,y1,s=20,c='r')\n",
    "        \n",
    "\n",
    "plt.savefig(f\"path.png\", format='png', dpi=400)  # Save the figure as a PNG image file\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.write()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}