regression-multiline-V4.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 random \n",
    "import operator\n",
    "from pyomo import environ as pe\n",
    "import os\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "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)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAI4CAYAAABndZP2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvpUlEQVR4nO3df3Rc533f+fcXPwhRoOKmlEmHckiRxzWdmEkMmqbqjdGQ65OsnJ71gnFrl1vmVJZJZQ/tRtzj0812z25iL9MT97RN2dphXFF1ZIc+YZIygrNdh45jk1tYqSOLhGS7UribUIIZsiIV+EcMGAYxmGf/mAsZAUlwAGLuzDzzfp1zz8XceQb44hHE+cxzn+feSCkhSZKUk65mFyBJkrTSDDiSJCk7BhxJkpQdA44kScqOAUeSJGWnp9kFNMLdd9+d7r333maXUbfJyUn6+/ubXUa27N/Gs48by/5tPPu4sRrZv2fPnv3LlNIrFx7PMuDce++9PPXUU80uo25nzpxh165dzS4jW/Zv49nHjWX/Np593FiN7N+IGLvRcU9RSZKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlx4AjSZKyY8CRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOwYcCRJUnZKCzgRcSgivhwRsxGRIuIDt2h/R0R8OCKuRsRURDwREfeVVK4kSWpjZY7gvBH4OnCxzvZHgPcBV4Bh4M3AZyPi7kYUJ0mS8lFawEkp/WxKaRfw9K3aRsQ64EGgCrw1pbQX+CRwF7XQI0mSdFM9zS7gJl4P9AIvpJSuFseeAvYBb7jRCyLiIeAhgPXr13PmzJnGV7lCJiYm2qredmP/Np593Fj2b+PZx43VjP5t1YCzvthPzDs2WexfdaMXpJQeAR4B2LFjR9q1a1fDiltpZ86coZ3qbTf2b+PZx41l/zaefdxYzejfVl1FdaXYr5l3bO7rF0uuRZIktZmWCDgR8YqIeF1E3FscehaYATZGxNxozpuK/TNl1ydJktpLaaeoImI/8BZge3FoqAg0w8DfAH6DWnh5Q0rpSkQ8BhwAPhcRXwXeSe2U1UfKqlmSJLWnMufgvAX4R/Me/1ixvVBsCz1MbRTnncBrgC8C708pvdTQKiVJUtsrLeCklB4AHlikyWML2k8B7y02SZKkurXEHBxJkqSVZMCRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOyUebNNSZJuaWx8kmMjFxgevczkdIX+vh6GBjZwYHALm9b2N7s8tQkDjiSpZZw+f5WDx88xM1ulUk0ATExXOPHkRU6evcTRfdvZvXVdk6tUO/AUlSSpJYyNT3Lw+DmmZmZfDjdzKtXE1MwsB4+fY2x8skkVqp0YcCRJLeHYyAVmZquLtpmZrfLoyPMlVaR2ZsCRJLWE4dHL143cLFSpJh4fvVRSRWpnBhxJUkuYnK7U1+5afe3U2ZxkLElaskasdOrv62GijpDTv8q3Lt2aIziSpCU5ff4q9x8Z4cSTF5mYrpD43kqn+4+McPr81WV936GBDfR0xaJterqCPQP3LOv7q7MYcCRJdWvkSqcDg1vo7V78bam3u4v9g5uX/L3VeQw4kqS6NXKl06a1/Rzdt53Vvd3XjeT0dAWre7s5um+7F/tTXQw4kqS6NXql0+6t6zh1aJC9Ozeypq+HCFjT18PenRs5dWjQi/ypbs7UkiTVrYyVTpvW9nN4aBuHh7Yt+3tIjuBIkurW31ff52JXOqnZDDiSpLq50kntwoAjSaqbK53ULgw4kqS6udJJ7cKAI0laElc6qR04C0yStGSudFKrcwRHkiRlx4AjSZKyY8CRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOwYcCRJUnYMOJIkKTsGHEmSlB0DjiRJyo4BR5IkZceAI0mSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlp6fZBUgrbWx8kmMjFxgevczkdIV/8qOz/NHwVzgwuIVNa/ubXZ4kqQQGHGXl9PmrHDx+jpnZKpVqAmA2JU48eZGTZy9xdN92dm9d1+QqJUmN5ikqZWNsfJKDx88xNTP7criZU6kmpmZmOXj8HGPjk02qUJJUFgOOsnFs5AIzs9VF28zMVnl05PmSKpIkNYsBR9kYHr183cjNQpVq4vHRSyVVJElqFgOOsjE5Xamv3bX62kmS2pcBR9no76tvznz/KufWS1LuDDjKxtDABnq6YtE2PV3BnoF7SqpIktQsBhxl48DgFnq7F/+T7u3uYv/g5pIqkiQ1iwFH2di0tp+j+7azurf7upGcnq5gdW83R/dt92J/ktQBDDjKyu6t6zh1aJC9Ozeypq+HCOiOYO/OjZw6NOhF/iSpQzjbUtnZtLafw0PbODy0DYAzZ85wcNe2JlclSSqTIziSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdlxmbgkSW1ubHySYyMXGB69zOR0hf6+HoYGNnBgcEvHXtzUgCNJUhs7ff4qB4+fY2a2SqWaAJiYrnDiyYucPHuJo/u2d+RFTj1FJUlSmxobn+Tg8XNMzcy+HG7mVKqJqZlZDh4/x9j4ZJMqbB4DjiRJberYyAVmZquLtpmZrfLoyPMlVdQ6DDiSJLWp4dHL143cLFSpJh4fvVRSRa3DgCNJUpuanK7U1+5afe1yYsCRJKlN9ffVt1aof1XnrSky4EiS1KaGBjbQ0xWLtunpCvYM3FNSRa2jtIATEXdExIcj4mpETEXEExFx3yLt3xARn4mI8Yj4TkQ8GxEHy6pXkqRWd2BwC73di7+V93Z3sX9wc0kVtY4yR3COAO8DrgDDwJuBz0bE3TdpPwz8FHAJ+DTwOuDXImJ3owuVJKkdbFrbz9F921nd233dSE5PV7C6t5uj+7Z35MX+SjkpFxHrgAeBKvDWlNLViKgA+6iFng8saN8L/GDx8H9MKX01Ip4C3gjcW0bNkqTm8Kq8S7N76zpOHRrk0ZHneXz0EpPXKvSv6mHPwD3sH9zcsX0WKS2+vGxFfkht1OXzwAsppc3FsYepjep8KqU0dIPX/CrwPwNfBv4/4GeKr3ellL55g/YPAQ8BrF+//o0nTpxoxK/SEBMTE6xZs6bZZWTL/m08+7ixOql/vz1d4Wvj3yElSHzv/SkIImDj2ju5q86JtUvRSX3cDI3s3927d59NKe1YeLysadXri/3EvGNzl1V81U1eMwzsAX602GaKY9++UeOU0iPAIwA7duxIu3btup16S3XmzBnaqd52Y/82nn3cWJ3Sv2Pjk9x/ZISpme6btlndW+HUoTev+KhEp/RxszSjf8uag3Ol2M+Pb3Nfv7iwcUSsBf6A2umoQeBvAk8DvwT8XKOKlCQ1j1fl1UoqK+A8S20EZmNEzI3mvKnYPxMRr4iI10XEvcWxzcCdxWu+lFL6BvBc8dwPlVSzJKlEXpVXK6mUgJNSugI8Vvy8z0XECWAvtVNWH6F2Kuo5aqegKL7+OtBbtP940R7gC2XULEkql1fl1Uoqc5n4w8BRavNxhoAvAj+VUnppYcOU0iTw08AfURux+fvAnwGHUkq/XVbBkqTyeFVeraTS/kpSSlPAe4tt4XOPURvhmX/sT4CfLKM2SVLzDQ1s4MSTFxc9TdWpV+XV0nmrBklSS/CqvFpJBhxJUkvwqrxaSQYcSVLLmLsq796dG1nT10MErOnrYe/OjZw6NMjureuaXaLahDO1JEktZdPafg4PbePw0LZml6I25giOJEnKjgFHkiRlx4AjSZKyY8CRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOwYcCRJUnYMOJIkKTsGHEmSlB0DjiRJyo4BR5IkZceAI0mSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlx4AjSZKyY8CRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdnpaXYBqhkbn+TYyAWGRy8zOV2hv6+HoYENHBjcwqa1/c0uT5KktmLAaQHfnq5w/5ERZmarVKoJgInpCieevMjJs5c4um87u7eua3KVkiS1D09RNdnY+CRfG/8OUzOzL4ebOZVqYmpmloPHzzE2PtmkCiVJaj8GnCY7NnKBlBZvMzNb5dGR58spSJKkDBhwmmx49DKJxRNOpZp4fPRSSRVJktT+DDhNNjldqa/dtfraSZIkA07T9ffVN8+7f5XzwSVJqpcBp8mGBjYQxKJterqCPQP3lFSRJEntz4DTZAcGtxCL5xt6u7vYP7i5nIIkScqAAafJNq3tZ+PaO1nd201P119POj1dwerebo7u2+7F/iRJWgIDTgu4q6+HU4cG2btzI2v6eoiANX097N25kVOHBr3InyRJS+TM1RaxaW0/h4e2cXhoW7NLkSSp7TmCI0mSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlxwv9SVKbGhuf5NjIBYZHLzM5XaG/r4ehgQ0cGNzi7V3U8Qw4ktSGTp+/ysHj55iZrVKpJgAmpiucePIiJ89e4ui+7d7mRR3NU1SS1GbGxic5ePwcUzOzL4ebOZVqYmpmloPHzzE2PtmkCqXmM+BIUps5NnKBmdnqom1mZqs8OvJ8SRVJrceAI0ltZnj08nUjNwtVqonHRy+VVJHUegw4ktRmJqcr9bW7Vl87KUcGHElqM/199a0P6V/lOhJ1LgOOJLWZoYEN9HTFom16uoI9A/eUVJHUegw4ktRmDgxuobd78X++e7u72D+4uaSKpNZjwJGkNrNpbT9H921ndW/3dSM5PV3B6t5uju7b7sX+1NEMOJLUhnZvXcepQ4Ps3bmRNX09RMCavh727tzIqUODXuRPHc8ZaJK0Qsq+dcKmtf0cHtrG4aFtK/69pXZnwJGkFeCtE6TW4ikqSbpN3jpBaj0GHEm6Td46QWo9BhxJuk3eOkFqPQYcSbpN3jpBaj0GHEm6Td46QWo9BhxJuk3eOkFqPQYcSbpN3jpBaj2lBZyIuCMiPhwRVyNiKiKeiIj7bvGaPRHxpaL9tyLiCxHx/WXVLEn18NYJUuspcwTnCPA+4AowDLwZ+GxE3H2jxhGxF/g94EeATwG/C9wF3FlCrZK0JN46QWotpcx4i4h1wINAFXhrSulqRFSAfdRCzwcWtA/gnxcP708pnSmjTkm6Hd46QWodkdLi125YkR8SsRv4PPBCSmlzcexhaqM6n0opDS1o/1rgPDAFnAH+DvAi8K9TSr92k5/xEPAQwPr169944sSJRvwqDTExMcGaNWuaXUa27N/Gs48by/5tPPu4sRrZv7t37z6bUtqx8HhZaxbXF/uJecfmrln+qhu0nztttRrYAvwOsBf4SERcSikNL3xBSukR4BGAHTt2pF27dt1+1SU5c+YM7VRvu7F/G88+biz7t/Hs48ZqRv+WNQfnSrGfH9/mvn7xBu1fmvf1z6aUHgQ+Vjx++wrXJkmSMlNWwHkWmAE2RsTcaM6biv0zEfGKiHhdRNxbHBsD/mrB95hbmjCBJEnSIkoJOCmlK8Bjxc/7XEScoHbKaQL4CLAHeI7a6ipSSteozc8B+EREfAx4NzALfLKMmiVJUvsq87rhD1MbxXkn8Brgi8D7U0ov1RZNXecwsAp4AHgX8FXgl1JKf1JKtZIkqW2VFnBSSlPAe4tt4XOPURvhmX+sAvzTYpMkSaqbt2qQJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOwYcCRJUnYMOJIkKTsGHEmSlB0DjiRJyo4BR5IkZceAI0mSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpO8sKOBHx4xGxbaWLkSRJWgl1BZyI+GhEPB01x4H/BDwTEQcbW54kSdLS1TuC81PAXwB3Au8CvgJ8G/j5BtUlSZK0bPUGnB8AxoAfLl7zs8BvAxsbVJckSdKy1RtwJoA3AHuBGeBPgVXAdxtTliRJ0vL11NnuDPAO4G8Dn0kpzUTEjwDnG1WYJEnSctUbcA5QCzM9wJGIWAX8PnCuUYVJkiQtV10BJ6X0TeB/j4g7gC0ppReB/7ORhUmSJC1XvcvEV0fEv6c2F+fLEXFvRPx5RPzjxpYnSZK0dPVOMj4MvLtoHymlF4BvAP+wQXVJkiQtW70B5+8BZ4Ffn3fsGeC1K16RJEnSbao34Hw/8Cxwbd6x1UCseEWSJEm3qd6A86fA/cAPAUTEe4Eh4LnGlCVJkrR89Qacfwa8EvjJ4vG/Be4APtSIoiRJkm5HXQEnpfT7wNuAU9RGbf4A+OniuCRJUkup90J/pJQ+A3ymgbVIkiStiLoCTkR87CZPpZTSe1awHkmSpNtW7wjOA0Die6um5r5OgAFHkiS1lHoDzgfnfd0N/AjwduBmIzuSJElNU++9qD648FhEPALcs+IVSZIk3aZ65+BsXHDo+4CtwI+teEWSJEm3qd5TVM/f5PhXVqoQSZKklVLvhf5iwTYF/DG1yceSJEktpd45OPUGIUmSpKa7acC5wbybG0opfW3lypEkSbp9i43g3GzezXzpFt9D6khj45McG7nA8OhlJqcr9Pf1MDSwgQODW9i0tr/Z5UlS9hYLJ7HIc0tpI3WU0+evcvD4OWZmq1SqCYCJ6QonnrzIybOXOLpvO7u3rmtylZKUt5vOrUkpddWzlVms1OrGxic5ePwcUzOzL4ebOZVqYmpmloPHzzE2PtmkCiWpMyw5oETEXRHxfXNbI4qS2tWxkQvMzFYXbTMzW+XRkXrOAEuSlquugBMRr46IT0fEd4BvAt8otq83sDap7QyPXr5u5GahSjXx+OilkiqSpM5U7wThXwPub2QhUg4mpyv1tbtWXztJ0vLUe4rq7wDngE9QWzl1H/DnwIMNqktqS/199X1m6F/l4kNJaqR6A86dwH8GrhSPvwL8P8D/1oiipHY1NLCBnq7FFxf2dAV7BrxPrSQ1Ur0B5+vAGuBFakvDTwJ7gA0NqktqSwcGt9Dbvfj/Vr3dXewf3FxSRZLUmeoNOF8F3gj8X0AVeBvw/cAfNqguqS1tWtvP0X3bWd3bfd1ITk9XsLq3m6P7tnuxP0lqsHonArwN6EopXYuI3cDPUDtd9ZGGVSa1qd1b13Hq0CCPjjzP46OXmLxWoX9VD3sG7mH/4GbDjSSVoN6A8/PAceBqSmkEGGlcSVL727S2n8ND2zg8tK3ZpUhSR6r3FNW/BP4iIv7viHhnRPQ1sihJkqTbUW/A+TjwbWqnqn4LeDEijkXEYMMqkyRJWqa6Ak5K6d3AeuC/B36zOPwe4HSD6pIkSVq2uu9FlVKqUFtN9SzwteKwdxOXJEktp65JxhFxCHgXsHPuEPAlaqeuJEmSWkq9q6h+tdhfpLaa6hMppfONKUmSJOn21BtwPlFsp1NKi98qWZIkqcnqCjgppQcaXIf014yNT3Js5ALDo5eZnK7Q39fD0MAGDgxu8UJ5kqRb8pbGajmnz1/l4PFzzMxWqVRrA4YT0xVOPHmRk2cvcXTfdnZvXdfkKiVJrazuVVRSGcbGJzl4/BxTM7Mvh5s5lWpiamaWg8fPMTY+2aQKJUntwICjlnJs5AIzs9VF28zMVnl05PmSKpIktaNFA05EDEfED5RVjDQ8evm6kZuFKtXE46OXSqpIktSObjWC83bg2Yj4uTKKkSanK/W1u1ZfO0lSZ7pVwBkBXgEcjYj/FBGvLaEmdbD+vvrmvfevcn68JOnmFg04KaWfAN4NvAS8BXgmIj4VEb9XbCfLKFKdY2hgAz1di98BpKcr2DNwT0kVSZLa0S0/BqeUPh4RF4DPAn3Ubrj58tONKkyd6cDgFk6evUSlOnvTNr3dXewf3FxiVZKkdrNowImIfuCDwD8GeoFLwG8AN3/3kW7DprX9HN23/brr4EBt5Ka3u4uj+7Z7sT9J0qJuNYLzHDB3LuCjwC+klL7d2JLU6XZvXcepQ4M8OvI8j49eYvJahf5VPewZuIf9g5sNN5KkW7pVwHk18KfAgZTSEyXUIwG1kZzDQ9s4PLSt2aVIktrQrVZR/TLwBsONJElqJ4uO4KSUfrGsQiRJklaKt2qQJEnZMeBIkqTslBZwIuKOiPhwRFyNiKmIeCIi7qvjdXsjIhXbkRJKlSRJba7MEZwjwPuAK8Aw8GbgsxFx981eEBGvBo4C3nhIkiTVrZSAExHrgAeBKvDWlNJe4JPAXdRCz41eE8DHgcuAt4SQJEl1K+uOha+ndiXkF1JKV4tjTwH7gDfc5DWHqN3/6r7i60VFxEPAQwDr16/nzJkzt1NvqSYmJtqq3nZj/zaefdxY9m/j2ceN1Yz+LSvgrC/2E/OOTRb7Vy1sHBHbgF8BfjGl9HRtMGdxKaVHgEcAduzYkXbt2nU79ZbqzJkztFO97cb+bTz7uLHs38azjxurGf1bVsC5UuzXzDs29/WLN2j/DmAV8BMRMQj8WHH87RExlVL6p40pU5JUr7HxSY6NXGB49DKT0xX6+3oYGtjAgcEt3lJFTVdWwHkWmAE2RsT6lNIV4E3Fc89ExCuAHwC+m1J6AYhie9uC77OZ2uRkSVITnT5/9bqb4k5MVzjx5EVOnr3E0X3b2b11XZOrVCcrZZJxEWgeK37e5yLiBLCX2imrjwB7qN3Yc7ho/4GUUsxt1CYbA/yblNKuMmqWJN3Y2PgkB4+fY2pm9uVwM6dSTUzNzHLw+DnGxidv8h2kxitzmfjD1JZ8rweGgC8CP5VSeqnEGiRJt+nYyAVmZquLtpmZrfLoyPMlVSRdr6xTVKSUpoD3FtvC5x6jNsJzs9c+ADzQmMp0uzwPL3WW4dHL143cLFSpJh4fvcThoW0lVSX9daUFHOXJ8/BS55mcru/aq5PXvEarmsd7UWnZPA8vdab+vvo+G/ev8jO0mseAo2XzPLzUmYYGNtDTtfj1yXq6gj0D95RUkXQ9A46WbSnn4SXl48DgFnq7F3/76O3uYv/g5pIqkq5nwNGyeR5e6kyb1vZzdN92Vvd2XzeS09MVrO7t5ui+7S4yUFMZcLRsnoeXOtfures4dWiQvTs3sqavhwhY09fD3p0bOXVo0MUFajrfebRsQwMbOPHkxUVPU3keXsrXprX9HB7a5lJwtSQDjpbtwOAWTp69RKU6e9M2nodXJ/LaUFLzeYpKy+Z5eOl6p89f5f4jI5x48iIT0xUS37s21P1HRjh9/mqzS5Q6ggFHt8Xz8NL3eG0oqXV4ikq3zfPwUs1Srg3l/y9SYzmCI0krxGtDSa3DERy1DSduqtV5bSipdTiCo7bgxE21A68NJbUOA45anhM31S68R5PUOgw4anne1FPtwns0Sa3DgKOW58RNtQuvDSW1DgOOWp4TN9VOvDaU1Bqc6aaW19/Xw0QdIceJm2oVXhtKaj5HcNTynLgpSVoqA45anhM3JUlLZcBRy3PipiRpqQw4agtO3JQkLYWzMtU2nLgpSaqXIziSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdlxmbgkdbix8UmOjVxgePQyk9MV+vt6GBrYwIHBLV5AU23LgCNJHez0+ascPH6OmdkqlWoCYGK6woknL3Ly7CWO7tvuhTTVljxFJUkdamx8koPHzzE1M/tyuJlTqSamZmY5ePwcY+OTTapQWj4DjiR1qGMjF5iZrS7aZma2yqMjz5dUkbRyDDiS1KGGRy9fN3KzUKWaeHz0UkkVSSvHgCNJHWpyulJfu2v1tZNaiQFHkjpUf19960z6V7keRe3HgCNJHWpoYAM9XbFom56uYM/APSVVJK0cA44kdagDg1vo7V78baC3u4v9g5tLqkhaOQYcSepQm9b2c3Tfdlb3dl83ktPTFazu7ebovu1e7E9tyYAjSR1s99Z1nDo0yN6dG1nT10MErOnrYe/OjZw6NOhF/tS2nDkmSR1u09p+Dg9t4/DQtmaXIq0YR3AkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOwYcCRJUnYMOJIkKTveqkGSpFsYG5/k2MgFhkcvMzldob+vh6GBDRwY3OLNSFuUAUeSpEWcPn+Vg8fPMTNbpVJNAExMVzjx5EVOnr3E0X3bvSlpCzLgqGP5iUzSrYyNT3Lw+DmmZmave65STVSqsxw8fo5Thwb9d6PFOAdHHen0+avcf2SEE09eZGK6QuJ7n8juPzLC6fNXm12ipBZwbOQCM7PVRdvMzFZ5dOT5kipSvQw46jjzP5HNDTfPqVQTUzO1T2Rj45NNqlBSqxgevXzdvxMLVaqJx0cvlVSR6mXAUcfxE5mkek1OV+prd62+diqPAUcdx09kkurV31ffVNX+VU5pbTUGHHUcP5FJqtfQwAZ6umLRNj1dwZ6Be0qqSPUycqrj9Pf1MFFHyOmkT2SuKJNu7MDgFk6evUSlev0qqjm93V3sH9xcYlWqhyM46jh+IvvrXFEm3dymtf0c3bed1b3d1/270dMVrO7t5ui+7X4QaEEGHHWcA4Nb6O1e/E+/Uz6RuaJMurXdW9dx6tAge3duZE1fDxGwpq+HvTs3curQoBf5a1GdMwYvFeY+kS28MinUPpH1dnd1zCeypawoOzy0raSqpNazaW0/h4e2+f9BG3EERx3JT2Q1riiTlCtHcNSx/ETmijJJ+XIER+pgXuNDUq4MOFIHc0WZpFwZcKQO5ooySbky4EgdzGt8SMqVAUfqcK4ok5QjZw5KckWZpOw4giNJkrJjwJEkSdkx4EiSpOwYcCRJUnYMOJIkKTsGHEmSlB0DjiRJyo4BR5IkZae0gBMRd0TEhyPiakRMRcQTEXHfIu2PRcSzETEREeMR8emIeH1Z9UqSpPZV5gjOEeB9wBVgGHgz8NmIuPsm7fcDfwX8VrF/G/CZiLij4ZVKkqS2VkrAiYh1wINAFXhrSmkv8EngLmqh50Z+PKX0t1NKB4DdxbF7gB9udL2SJKm9RUqp8T8kYjfweeCFlNLm4tjD1EZ1PpVSGrrF618LnKcWkF6dUvqvN2jzEPAQwPr169944sSJlfwVGmpiYoI1a9Y0u4xs2b+NZx83lv3bePZxYzWyf3fv3n02pbRj4fGybra5vthPzDs2WexftdgLI2IN8Fjx8F/dKNwApJQeAR4B2LFjR9q1a9dyay3dmTNnaKd6243923j2cWPZv41nHzdWM/q3rDk4V4r9/Pg29/WLN3tRMT/n89Tm6xwDfqEh1UmSpKyUFXCeBWaAjRExN5rzpmL/TES8IiJeFxH3zr0gIjYBTxTtPpRSeiiVcT5NkiS1vVICTkrpCrXTTF3A5yLiBLCX2imrjwB7gOeora6a88fAa4GvAasj4kix7SyjZkmS1L7KmoMD8DC1UZx3Aq8Bvgi8P6X0UkTcqP2GYr+xeO2cp4EnG1emJElqd6UFnJTSFPDeYlv43GN8byLx3LEbph5JkqRb8VYNkiQpOwYcSZKUHQOOJEnKjgFHkiRlx4AjSZKyY8CRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOwYcCRJUnYMOJIkKTsGHEmSlB0DjiRJyo4BR5IkZceAI0mSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlx4AjSZKyY8CRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOwYcCRJUnYMOJIkKTsGHEmSlB0DjiRJyo4BR5IkZceAI0mSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlp6fZBShPY+OTHBu5wPDoZSanK/T39TA0sIEDg1vYtLa/2eVJkjJnwNGKO33+KgePn2NmtkqlmgCYmK5w4smLnDx7iaP7trN767omVylJypmnqLSixsYnOXj8HFMzsy+HmzmVamJqZpaDx88xNj7ZpAolSZ3AgKMVdWzkAjOz1UXbzMxWeXTk+ZIqkiR1IgOOVtTw6OXrRm4WqlQTj49eKqkiSVInMuBoRU1OV+prd62+dpIkLYcBRyuqv6++eev9q5zfLklqHAOOVtTQwAZ6umLRNj1dwZ6Be0qqSJLUiQw4WlEHBrfQ2734n1Vvdxf7BzeXVJEkqRMZcLSiNq3t5+i+7azu7b5uJKenK1jd283Rfdu92J8kqaEMOFpxu7eu49ShQfbu3Miavh4iYE1fD3t3buTUoUEv8idJajhneqohNq3t5/DQNg4PbWt2KZKkDuQIjiRJyo4BR5IkZceAI0mSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjlcyliRpEWPjkxwbucDw6GUmpyv09/UwNLCBA4NbvK9eCzPgSJJ0E6fPX+Xg8XPMzFapVBMAE9MVTjx5kZNnL3F033bvr9eiPEUlSdINjI1PcvD4OaZmZl8ON3Mq1cTUzCwHj59jbHyySRVqMQYcSZJu4NjIBWZmq4u2mZmt8ujI8yVVpKUw4EiSdAPDo5evG7lZqFJNPD56qaSKtBQGHEmSbmByulJfu2v1tVO5DDiSJN1Af19963D6V7lepxWVFnAi4o6I+HBEXI2IqYh4IiLuW6n2kiStpKGBDfR0xaJterqCPQP3lFSRlqLMEZwjwPuAK8Aw8GbgsxFx9wq1lyRpxRwY3EJv9+Jvk73dXewf3FxSRVqKUgJORKwDHgSqwFtTSnuBTwJ3UQsxt9VekqSVtmltP0f3bWd1b/d1Izk9XcHq3m6O7tvuxf5aVKS0+AzxFfkhEbuBzwMvpJQ2F8cepjZK86mU0tDttC+efwh4CGD9+vVvPHHiRIN+m5U3MTHBmjVrml1GtuzfxrOPG8v+bbzF+vjabJW//PY1vvmda8ymRHcEf+POVdx91ypW3WKERzWN/BvevXv32ZTSjoXHy5oZtb7YT8w7NndlpFetQHtSSo8AjwDs2LEj7dq1a1mFNsOZM2dop3rbjf3bePZxY9m/jWcfN1Yz+res6Hml2M+Pb3Nfv7gC7SVJkl5WVsB5FpgBNkbE3OjMm4r9MxHxioh4XUTcW0/7MgqWJEntq5SAk1K6AjxW/LzPRcQJYC+1U1AfAfYAz1FbLVVPe0mSpJsq8+pED1MblXkn8Brgi8D7U0ovRdzwOgM3bV9OuZIkqV2VFnBSSlPAe4tt4XOPURuxqau9JEnSYlzfJkmSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlx4AjSZKyY8CRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUeSJGXHgCNJkrJjwJEkSdkx4EiSpOwYcCRJUnYMOJIkKTsGHEmSlB0DjiRJyo4BR5IkZceAI0mSsmPAkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlx4AjSZKyY8CRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpSdSCk1u4YVFxEvAWPNrmMJ7gb+stlFZMz+bTz7uLHs38azjxurkf27KaX0yoUHsww47SYinkop7Wh2HbmyfxvPPm4s+7fx7OPGakb/eopKkiRlx4AjSZKyY8BpDY80u4DM2b+NZx83lv3bePZxY5Xev87BkSRJ2XEER5IkZceAI0mSsmPAkSRJ2THgNFhE3BERH46IqxExFRFPRMR9K9W+0y2jf49FxLMRMRER4xHx6Yh4fZk1t5vl/k1GxN6ISMV2pIRS29Jy+jci9kTEl4r234qIL0TE95dVc7tZxr8Tb4iIzxT/Rnyn+DfjYJk1t5OIOBQRX46I2eL/9w/con0573MpJbcGbsBHgQR8BfgtoAr8FXD3SrTv9G0Z/ZuALwLHgOeLx38B3NHs36VVt+X8TQKvBr4BzBSvPdLs36NVt2X8De8t2n8XOAE8CjwD3NPs36VVt2X08QtF+y8D/6Fon4Ddzf5dWnEDfhM4M6/fPrCS/z2WXVezOybnDVgHXANmgXXz/hBu+Aew1Padvi2nv4D/Zt7X9xZtE7C92b9PK27L7OMAPgf8l+IN2ICzQv1b9O3Xiud3Nbv+dtiW0ce9RdsEbCuOPVU8fnezf59W3oDhW71flfk+5ymqxno9tf9ZvpZSuloce6rYv2EF2ne6JfdXSumP5z1cVeyrwH9tRIEZWM7f5CHgLcA/pDbKoJtbav/+LeAHgSngfylOtf5ZRLy34ZW2ryX1cUppBvg3xcNPRsR/ALZTGyV7vLGldoTS3ucMOI21vthPzDs2WexftQLtO92y+ysi1gCPFQ//VUrJgHNjS+rjiNgG/ArwiymlpxtbWhaW+jd8d7FfDWwBfge4B/hIRAw1osAMLOffiWFqp1t+FHgHUCmOfXvFq+s8pb3PGXAa60qxXzPv2NzXL65A+063rP6KiLuBzwNvpjYX5xcaUl0eltrH76A2MvYTEfEfgbcWx98eEb/SmBLb2lL796V5X/9sSulB4GPF47evcG25WFIfR8Ra4A+oncIeBP4m8DTwS8DPNarIDlLa+5wBp7GepTbJcmNEzKXWNxX7ZyLiFRHxuoi4t572ZRTcZpbav0TEJuCJot2HUkoPpeIksG5oqX0cxfY24O9Sm2wMsJlaoNRft9T+HaM2GXO+KPYT6EaW2sebgTuL13wppfQN4LniuR8qqeZsNPV9rtmTknLfqN1/IwFfpTbhskptmPOVwAPFc0/X077Zv0srbsvo30vFsTHgyLxtZ7N/l1bdltrHC177GE4yXtH+BT5YHHuO2ujNFLVTKPc1+3dp1W0pfQz0A+PFsS8AH6c2KTYB72r279KKG7C/+H99bgL808XjoWa+z/WgRnuYWlp9J/AaakuU359SeikiltS+nHLbzlL7d0Ox31i8ds7TwJONK7OtLbWPtTRL7d/D1E4DPgC8i9qbxC+llP6klGrbU919nFKajIifBn6Z2uTi7cCfAf8upfTbpVbdPt4C/KN5j3+s2F4otoVKeZ/zZpuSJCk7zsGRJEnZMeBIkqTsGHAkSVJ2DDiSJCk7BhxJkpQdA44kScqOAUdSdiLigYhIEXGm2bVIag4DjqSGiYh7I+LbRdi4f97xf1kc+1JEdDezRkl5MuBIapiU0gt872amj0TE90XETuAQtcvfvzulNNuk8iRlzIAjqdF+HTgN/CDwb6ndP6kb+OWU0lcXNo6INRExERGViFhXHLtz3rFXRcRPRsRoRHwrImYiYiwiPnizAm50yioiXiiO7Zr3Mz4UEX8WEZMRcS4ihua1/8mIOFs8963i+Z9ZiQ6StPIMOJIaKtXuB/MeYJLa/WpeT+3eXx+6SfsJ4HFqIegdxeG/S+0miH+YUnoRuAf4S2o36vtN4C7gFyPiH9xGqf+e2mjTt4CT1ALZ780FIOA3qN1f52SxVYFtt/HzJDWQAUdSw6WUngf+47xDH00pzSzykk8U+3cV+3cuOP4J4F8DF6kFkj8vjv+3y6kvIl4J/ANqoeWPga8D/wUI4H8qmvUC3wV+H/jnwE7gny3n50lqPO8mLqnhIuI+4O8DiVpo+D8i4kRK6Vs3ecnngL8ABiPiNcBPUwsyw8Xzvw48dIPXvXIJZc2f3Hxvse8C3reg3WuK/c8B/wL43eLxeNH2xBJ+pqSSOIIjqaEioo/a6Z0u4APA56mdYvrVm70mpVQFPlm85mPAncDvppS+WzSZG9l5gFpQ+fW5H3eTbzlZ7L+vqGkt8Kp5z79Q7K8Br0wpRUopgFXAnuK5P0gp/S3gbuDvAWtxBEdqWQYcSY32AeCHgGeAXwH2UwscD0bEf7fI6z5e7AeL/SfmPXel2P88tTk4D9yihmeojR69ISJ+Dfg080awU0ovAb9DLdD8SUR8NCJ+l9opsPcUzUYj4tPUQs2DxbFv3uLnSmoSA46khomINwL/BKhQWxI+U8zH+V+LJo9ExF03em1K6TngqeLh88AX5j29H/hT4IepTTD+d4vVkVL6f4ufOQ78D8AfAl9b0Ow91CY+V6kFph8H/jNwqnj+j4Ct1CZKvwU4U9QhqQVFbYGDJElSPhzBkSRJ2THgSJKk7BhwJElSdgw4kiQpOwYcSZKUHQOOJEnKjgFHkiRlx4AjSZKy8/8DlAz2Rq9MlBgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.figure(figsize=(8,8))\n",
    "plt.scatter(df['x'], df['y'], s= 100)\n",
    "plt.xlim(-0.1,1.05)\n",
    "plt.ylim(-0.1,1.05)\n",
    "plt.xlabel('  X values ', fontsize=12, fontweight='bold')\n",
    "plt.ylabel('  Y values ', fontsize=12, fontweight='bold')\n",
    "plt.xticks( fontsize=12, fontweight='bold')\n",
    "plt.yticks( fontsize=12, fontweight='bold')\n",
    "plt.grid()\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.savefig('Dots '+  '.jpg', format='jpg', dpi=400)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BM= 4\n",
    "K= 3\n",
    "M = 3\n",
    "plt.figure(figsize=(8,8))\n",
    "counter =0 \n",
    "for M in [3]:\n",
    "    for K in [3]:\n",
    "        counter +=1\n",
    "        model = AbstractModel()\n",
    "        model.i = RangeSet(N)\n",
    "        model.m = RangeSet(M)\n",
    "        model.k = RangeSet(K) # number of lines \n",
    "        model.target = Param(model.i, initialize=0, mutable=True)\n",
    "        model.r = Var(model.i,bounds=(0, 1000) , within = Reals)\n",
    "        model.c = Var(model.k, model.m, bounds=(-100, 100) , within = Reals)\n",
    "        model.U = Var(model.i, model.k,bounds=(0, 1) , within = Binary)\n",
    "        def rule_C1(model,i,k):\n",
    "            return sum(model.c[k,m]*df.loc[i-1,'x']**(m-1) for m in model.m ) - df.loc[i-1,'y'] <= model.r[i] + BM*(1-model.U[i,k])\n",
    "        model.C1 = Constraint(model.i, model.k, rule=rule_C1)\n",
    "        def rule_C2(model,i,k):\n",
    "            return sum(model.c[k,m]*df.loc[i-1,'x']**(m-1) for m in model.m ) - df.loc[i-1,'y'] >= -model.r[i] - BM*(1-model.U[i,k])\n",
    "        model.C2 = Constraint(model.i,model.k,  rule=rule_C2)\n",
    "        def rule_C3(model,i):\n",
    "            return sum(model.U[i,k] for k in model.k) == 1\n",
    "        model.C3 = Constraint(model.i,  rule=rule_C3)\n",
    "        def rule_C4(model,k):\n",
    "            if k>1:\n",
    "                return sum(model.U[i,k-1] for i in model.i) <= sum(model.U[i,k] for i in model.i)\n",
    "            else:\n",
    "                return Constraint.Skip\n",
    "        model.C4 = Constraint(model.k,  rule=rule_C4)\n",
    "        \n",
    "        def rule_C5(model,k):\n",
    "            return sum(model.U[i,k] for i in model.i) >= 1\n",
    "        model.C5 = Constraint(model.k,  rule=rule_C5)\n",
    "        \n",
    "        def rule_of(model):\n",
    "            return sum(model.r[i] for i in model.i)\n",
    "        model.obj = Objective(rule=rule_of, sense=minimize)\n",
    "\n",
    "        instance = model.create_instance()\n",
    "        #solver_manager = pe.SolverManagerFactory('neos')\n",
    "        #results = solver_manager.solve(instance, solver='cplex', options={'mipgap':0.5})\n",
    "        \n",
    "        opt = SolverFactory('cbc')\n",
    "        opt.options = {'sec': 360, 'threads': 2, 'ratio': 0.5}\n",
    "        results = opt.solve(instance)\n",
    "\n",
    "\n",
    "        print('OF= ',value(instance.obj))\n",
    "        plt.scatter(df['x'], df['y'], s= 100)\n",
    "\n",
    "        for k in instance.k:\n",
    "            X=np.linspace(0,1,100)\n",
    "            Y = [sum(value(instance.c[k,m])*x**(m-1) for m in instance.m ) for x in X ]\n",
    "            plt.plot(X,Y)\n",
    "            xdata = [df.loc[i-1,'x'] for i in instance.i if value(instance.U[i,k])>0]\n",
    "            ydata = [df.loc[i-1,'y'] for i in instance.i if value(instance.U[i,k])>0]\n",
    "            plt.scatter(xdata,ydata, s= 100, zorder =2)\n",
    "\n",
    "\n",
    "        plt.xlim(-0.1,1.05)\n",
    "        plt.ylim(-0.1,1.05)\n",
    "        plt.xlabel('  X values ', fontsize=12, fontweight='bold')\n",
    "        plt.ylabel('  Y values ', fontsize=12, fontweight='bold')\n",
    "        plt.xticks( fontsize=12, fontweight='bold')\n",
    "        plt.yticks( fontsize=12, fontweight='bold')\n",
    "        plt.grid()\n",
    "        plt.tight_layout()\n",
    "        plt.savefig('Multi-line ' + str(K+(M-1)*3)+  '.jpg', format='jpg', dpi=400)\n",
    "        \n",
    "\n",
    "plt.savefig('Multi-line '+  '.jpg', format='jpg', dpi=400)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "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
}