aemarkov/bomj-lqr.ipynb

## bomj-lqr.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse import csc_matrix\n",
    "from scipy.sparse.linalg import spsolve, lsqr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ТАУ для долбоебов, часть I \n",
    "# Линейно-квадратичный регулятор (LQR)\n",
    "\n",
    "По мотивам статьи [Математика на пальцах: линейно-квадратичный регулятор](https://habr.com/post/277671/). Я не понял, поэтому попытался понять.\n",
    "\n",
    "---\n",
    "```\n",
    "Задача, в которой динамическая система описывается линейными дифференциальными уравнениями, а показатель качества представляет - квадратными, называется задачей линейно-квадратичного управления. \n",
    "```\n",
    "---\n",
    "\n",
    "Дискретная система, описанная в пространстве состояний:\n",
    "$$\n",
    "\\begin{equation}\n",
    "    x_{k+1} = Ax_k + Bu_k\\\\\n",
    "    y_{k+1} = Cx_k + Du_k\\\\\n",
    "\\end{equation}\n",
    "$$\n",
    "$\n",
    "\\begin{align*}\n",
    "    \\text{где } & x - \\text{вектор состояния,}\\\\\n",
    "               & u - \\text{вектор управления,}\\\\\n",
    "               & y - \\text{вектор выходов.}\n",
    "\\end{align*}\n",
    "$\n",
    "\n",
    "Критерий оптимальности (perfomance index):\n",
    "$$\n",
    "\\begin{equation}\n",
    "    J = \\sum_{k=0}^{\\infty}{(x_k^T Q_{x_k} + u_k^T R_{u_k})}\n",
    "\\end{equation}\n",
    "$$\n",
    "$\n",
    "\\begin{align*}\n",
    "    \\text{где } & Q, R - \\text{матрицы коэффициентов.}                \n",
    "\\end{align*}\n",
    "$\n",
    "\n",
    "Минимизируем критерий оптимальности, и тогда оптимальное управление будет:\n",
    "$$\n",
    "\\begin{equation}\n",
    "    u_k = -Fx_k\n",
    "\\end{equation}\n",
    "$$\n",
    "$\n",
    "\\begin{align*}\n",
    "    \\text{где } & F - \\text{матрица коэффициентов, находимая путем решения уравнений Риккати}\n",
    "\\end{align*}\n",
    "$\n",
    "\n",
    "К уравнениям Риккати мы еще вернемся (может быть), но пока зайдем с другой стороны. \n",
    "\n",
    "Рассмотрим уравнение линейной регрессии (внимание: буквы не имеют отношения к тем, что были раньше):\n",
    "$$\n",
    "\\begin{equation}\n",
    "    y = Xb + \\varepsilon\n",
    "\\end{equation}\n",
    "$$\n",
    "$\n",
    "\\begin{align*}\n",
    "    \\text{где } & y - \\text{вектор наблюдений (выходы),} \\\\\n",
    "                & X - \\text{матрица факторов (входы),} \\\\\n",
    "                & b - \\text{искомые коэффициенты,} \\\\\n",
    "                & \\varepsilon - \\text{ошибка}.\n",
    "\\end{align*}\n",
    "$\n",
    "\n",
    "Это неиллюзорно похоже на формулу оптимального управления. Поэтому попробуем в начале найти оптимальное управление с помощью метода наименьших квадратов.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "### Метод наименьших квадратов на Python\n",
    "\n",
    "Простейший пример метода наименьших квадратов и решение его с помощью `scipy`.\n",
    "\n",
    "Надо найти прямую, проходящую через несколько точек (линейная регрессия). Простое уравнение прямой (не может быть вертикальной, зато просто):\n",
    "$$\n",
    "y = kx + b\n",
    "$$\n",
    "\n",
    "Дано несколько точек $(x_1, y_1), (x_2, y_2) ... (x_n, y_n)$. Составим систему:\n",
    "$$\n",
    "\\begin{cases}\n",
    "    x_1 k + b = y_1 \\\\\n",
    "    x_2 k + b = y_2 \\\\\n",
    "    ...             \\\\\n",
    "    x_n k + b = y_n\n",
    "\\end{cases}r\n",
    "$$\n",
    "\n",
    "В матричной форме:\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "    x_1& 1 \\\\\n",
    "    x_2& 1 \\\\\n",
    "    x_n& 1\n",
    "\\end{bmatrix}\n",
    "\\cdot\n",
    "\\begin{bmatrix}\n",
    "    k\\\\\n",
    "    b\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "    y_1\\\\\n",
    "    y_2\\\\\n",
    "    y_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "$$\n",
    "AX=B\n",
    "$$\n",
    "\n",
    "И решим эту систему методом наименьших квадратов с помощью `scipy.sparse.linalg.lsqr`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VFX+//HXAUISQgDpJYRQElpoITRxEVcRRKVYVhBFbHGLrquu+90vQUDKqmtbu0ZEyrJi+QpEAUERBUWUIphCCYQAoRNqCAkp5/fHRH6BTSCQSW5m5v18PHhk7r0ncz43Ce+5c++dc4y1FhER8S1VnC5AREQqnsJfRMQHKfxFRHyQwl9ExAcp/EVEfJDCX0TEByn8RUR8kMJfRMQHKfxFRHxQNacLKEn9+vVtWFiY02WIiHiUdevWHbbWNrhYu0ob/mFhYaxdu9bpMkREPIoxZmdp2um0j4iID1L4i4j4IIW/iIgPqrTn/IuTm5tLeno62dnZTpfilQICAggJCcHPz8/pUkSknHlU+KenpxMcHExYWBjGGKfL8SrWWjIyMkhPT6dly5ZOlyMi5cyjTvtkZ2dTr149BX85MMZQr149vasS8REeFf6Agr8c6Wcr4js8LvxFRLzZkqT9fLhmV7n3o/Avo4kTJ/LCCy+UuH3+/PkkJydXYEUi4onSj2bxwMy1PDR7HR+tTaegoHznV/fu8J8zB8LCoEoV19c5cyq8BIW/iFxIbn4BcSu2M+ClFXy/7TBjB7djbkxvqlQp39Ow3hv+c+ZATAzs3AnWur7GxLjlBWDq1KlERERw1VVXsWXLFgDeffddevToQZcuXbj11lvJyspi1apVxMfH8+STT9K1a1e2b99ebDsR8U3rdh7l5te+4x+LNtO3TX2+euJqYvq1xq9q+Uez94Z/bCycH6xZWa71ZbBu3Trmzp3Lhg0bWLRoEWvWrAHglltuYc2aNWzcuJH27dvz3nvvceWVVzJkyBCef/55NmzYQOvWrYttJyK+5XhWLmPnJXDb26s4fjqXd+7uzrR7omlWJ7DCavCo+/wvya4SLpiUtL6UVq5cyfDhw6lRowYAQ4YMASAxMZFx48Zx7NgxMjMzGThwYLHfX9p2IuJ9rLUs2LCXKQuTOZqVy/19W/LYgAiC/Cs+ir03/ENDXad6iltfDsaMGcP8+fPp0qULM2bM4JtvvilTOxHxLjsOn+Kp+Yl8t+0wXZrXYeZ9kXRsWtuxetxy2scYM90Yc9AYk1jC9lHGmF+MMQnGmFXGmC7u6PeCpk6FwqPzs2rUcK0vg379+jF//nxOnz7NyZMn+eyzzwA4efIkTZo0ITc3lzlFrisEBwdz8uTJs8sltRMR75STl8+/vtrKwH+tYGP6MSYPi+TTP1xZfPBX4E0q7jrynwG8DswqYfsO4Gpr7VFjzA1AHNDLTX0Xb9Qo19fYWNepntBQV/D/uv4yRUVFcccdd9ClSxcaNmxIjx49AJg8eTK9evWiQYMG9OrV62zgjxgxggcffJBXX32VTz75pMR2IuJ9Vm07zLj5iaQePsXNXZry1E3taRgcUHzjX29S+fVa5a83qUCZc6s4xlr33EtqjAkDPrfWRl6k3RVAorW22YXaRUdH2/Mnc9m0aRPt27cvY6VyIfoZi5Td4cwcpi7cxLyf99CiXg0mD42kX8RFJtcKCyv+VHWLFpCWVuq+jTHrrLXRF2vnxDn/+4HFxW0wxsQAMQCh5XRuXkSkvBQUWD5cu5tnF28m60wej/y2DX+6pg0BflUv/s3ldJNKSSo0/I0x1+AK/6uK226tjcN1Sojo6Ojy/XibiIgbbd5/gth5iazbeZReLesydXgn2jSsWfonqOCbVCos/I0xnYFpwA3W2oyK6ldEpDxlncnjlWUpvLdyB7UC/Xjh9i7cGtXs0gdKnDr13HP+4JabVEpSIeFvjAkFPgXuttZurYg+RUTK27JNBxi/IIk9x05zR3Rz/n5DO64Iqn55T1ZON6mUxC3hb4z5AOgP1DfGpAMTAD8Aa+3bwHigHvBm4athXmkuSIiIVEb7jp/m6fhkvkjaT3jDmnz0UB96tqxb9iceNarcwv58bgl/a+3Ii2x/AHjAHX2JiDglL7+AmT/s5KWlW8i3lr8NassDV7WiejXPGynHez/h60PGjx9Pv379uO6665wuRcRrbdx9jNj5CSTuOUH/tg2YPDSS5nVrXPwbKymFfxlYa7HWUqXK5b3q5+XlUa1a2X8FkyZNKvNziEjxTmTn8uKSLcxavZMGNf15484oBndq7PEz33neexWHpaWl0bZtW0aPHk1kZCSzZ8+mT58+REVFcfvtt5OZmQnAokWLaNeuHd27d+fPf/4zN910E+Ca/OXuu++mb9++3H333eTn5/Pkk0/So0cPOnfuzDvvvAPAvn376NevH127diUyMpKVK1eSn5/PmDFjiIyMpFOnTrz88suAa7ygTz75BIBly5bRrVs3OnXqxH333UdOTg4AYWFhTJgwgaioKDp16sTmzZsr+kcn4lGstXz+y16ue/FbZq3eyT19wlj2xNXc2LmJxwc/ePCR/9OfJZG894Rbn7ND01pMuLnjRdulpKQwc+ZM2rRpwy233MJXX31FUFAQzz33HC+99BJ/+9vfeOihh1ixYgUtW7Zk5MhzL4kkJyfz3XffERgYSFxcHLVr12bNmjXk5OTQt29frr/+ej799FMGDhxIbGws+fn5ZGVlsWHDBvbs2UNiomsIpWPHjp3zvNnZ2YwZM4Zly5YRERHB6NGjeeutt/jLX/4CQP369Vm/fj1vvvkmL7zwAtOmTXPTT07Eu+zKyOKpBYl8u/UQkc1qMe2eaDqH1HG6LLfSkf9laNGiBb1792b16tUkJyfTt29funbtysyZM9m5cyebN2+mVatWtGzZEuC/wn/IkCEEBrrG7V66dCmzZs2ia9eu9OrVi4yMDFJSUujRowfvv/8+EydOJCEhgeDgYFq1akVqaiqPPPIIX3zxBbVq1Trnebds2ULLli2JiIgA4J577mHFihVnt99yyy0AdO/enbRL+Li4iK84k1fAG8u3MeDlb1m38ygTbu7Agj9d5XXBDx585F+aI/TyEhQUBLjeFg4YMIAPPvjgnO0bNmwo1ff/+hyvvfZaseP6r1ixgoULFzJmzBgef/xxRo8ezcaNG1myZAlvv/02H330EdOnTy913f7+/gBUrVqVvLy8Un+fiC/4accRYuclkHIwkxsiGzPh5o40rl3CIGxeQEf+ZdC7d2++//57tm3bBsCpU6fYunUrbdu2JTU19ezR9YcffljicwwcOJC33nqL3NxcALZu3cqpU6fYuXMnjRo14sEHH+SBBx5g/fr1HD58mIKCAm699VamTJnC+vXrz3mutm3bkpaWdrae2bNnc/XVV5fDnot4jyOnzvDkxxv53Ts/cDo3n+ljonnrru5eHfzgwUf+lUGDBg2YMWMGI0eOPHthdcqUKURERPDmm28yaNAggoKCzg77XJwHHniAtLQ0oqKisNbSoEED5s+fzzfffMPzzz+Pn58fNWvWZNasWezZs4d7772XgoICAJ555plznisgIID333+f22+/nby8PHr06MHvf//78vsBiHgway2frEvnH4s2cTI7j99f3ZpHrw0nsHopBmHzAm4b0tndPH1I58zMTGrWrIm1lj/96U+Eh4fz2GOPOV3WRXnSz1jkcm07eJKx8xL5accRoltcwdThnWjbONjpstyiMg/p7BPeffddZs6cyZkzZ+jWrRsPPfSQ0yWJ+Lzs3Hxe/3ob76zYTo3q1Xj2lk78Lro5Vap4/q2bl0rhX04ee+wxjzjSF/EV3249xFPzE9l1JItbopoRO7g99Wr6O12WYzwu/K21XvEBi8qosp4CFCmLgyeymfR5Mp//so9WDYL4z4O9uLJ1fafLcpxHhX9AQAAZGRnUq1dPLwBuZq0lIyODgADvvsNBfEd+gWXOjzt5/ost5OQX8PiACB66uhX+1Xzjgu7FeFT4h4SEkJ6ezqFDh5wuxSsFBAQQEhLidBkiZZa45zix8xLYmH6cq9rUZ/KwSFrWD7r4N/oQjwp/Pz+/s5+aFRE5X2ZOHi8t3cqMVTuoG+TPKyO6MqRLU50pKIZHhb+ISHGstSxJOsDTnyWx/0Q2d/YM5W+D2lE70M/p0iothb+IeLT0o1lMWJDEss0Hadc4mDdGRREVeoXTZVV6Cn8R8Ui5+QVM/24H//oqBYDYwe25t28Y1apq1JrScNccvtOBm4CD1trIYrYb4BVgMJAFjLHWrj+/nYhIaazbeZTYeQls3n+S69o34umhHWlWJ9DpsjyKu478ZwCvA7NK2H4DEF74rxfwVuFXEZFSO56Vy7NfbOaDn3bRtHYAcXd35/qOjZ0uyyO5awL3FcaYsAs0GQrMsq5PEa02xtQxxjSx1u5zR/8i4t2stSzYsJcpC5M5mpXLA1e15LEBEQT568z15aqon1wzYHeR5fTCdeeEvzEmBogBCA0NraDSRKQySz2UyVMLEvl+WwZdm9dh5n2RdGxa2+myPF6letm01sYBceAa1dPhckTEQdm5+bz97XbeXL4df78qTB4WyZ09Q6nqg4OwlYeKCv89QPMiyyGF60RE/suqbYcZNz+R1MOnGNKlKeNuak/DYA094k4VFf7xwMPGmLm4LvQe1/l+ETnf4cwcpi7cxLyf99CiXg1m3deTfhENnC7LK7nrVs8PgP5AfWNMOjAB8AOw1r4NLMJ1m+c2XLd63uuOfkXEOxQUWOau2c2zizdxOjefP/+2DX+8pg0BfhqErby4626fkRfZboE/uaMvEfEum/efYOynCazfdYzereoyZVgn2jSs6XRZXq9SXfAVEd+RdSaPV75KYdp3O6gd6MeLt3fhlqhmGoStgij8RaTCLdt0gPELkthz7DQjejTnfwa144qg6k6X5VMU/iJSYfYdP83E+CSWJB0golFNPv59H3qE1XW6LJ+k8BeRcpeXX8DMH3by0tIt5FvL3wa15YGrWlG9mgZhc4rCX0TK1Ybdx4idl0DS3hP0b9uAyUMjaV63htNl+TyFv4iUixPZubywZAuzV++kYbA/b46K4obIxrqgW0ko/EXEray1fP7LPiZ9nkxGZg739AnjiesjCA7QrFqVicJfRNxmV0YW4xYksmLrITo1q81790TTOaSO02VJMRT+IlJmZ/IKeHdlKq8uS8GvahUm3NyB0X3CNAhbJabwF5Ey+TE1g9j5iWw7mMngTo0Zf1NHGtfWIGyVncJfRC7LkVNneGbRJj5el07IFYG8P6YH17Rr6HRZUkoKfxG5JNZaPl6XzjOLNnEyO48/9G/Nn38bTmB1DcLmSRT+IlJqKQdOEjs/kZ92HCG6xRVMHd6Jto2DnS5LLoPCX0QuKjs3n9e+TiFuRSpB/tV47tZO3N69OVV0QddjKfxF5IK+2XKQ8QuS2HUki1ujQhg7uB31avo7XZaUkcJfRIp18EQ2T3+ezMJf9tGqQRAfPNibPq3rOV2WuInCX0TOkV9gmfPjTp7/Ygs5+QU8MSCCmKtb4V9NF3S9ibumcRwEvAJUBaZZa589b3soMBOoU9jm79baRe7oW0TcJ3HPcWLnJbAx/Ti/Ca/P5KGRhNUPcrosKQdlDn9jTFXgDWAAkA6sMcbEW2uTizQbB3xkrX3LGNMB15y+YWXtW0TcIzMnj5eWbmXGqh3UDfLnlRFdGdKlqQZh82LuOPLvCWyz1qYCGGPmAkOBouFvgVqFj2sDe93Qr4iUkbWWJUn7mRifzIGT2YzqFcqTA9tRO1CDsHk7d4R/M2B3keV0oNd5bSYCS40xjwBBwHVu6FdEymD3kSwmxiexbPNB2jepxZt3RREVeoXTZUkFqagLviOBGdbaF40xfYDZxphIa21B0UbGmBggBiA0NLSCShPxLbn5Bbz33Q5e+SoFYyB2cHvu7RtGtaqaVcuXuCP89wDNiyyHFK4r6n5gEIC19gdjTABQHzhYtJG1Ng6IA4iOjrZuqE1Eili38whjP01ky4GTDOjQiIlDOtKsTqDTZYkD3BH+a4BwY0xLXKE/ArjzvDa7gGuBGcaY9kAAcMgNfYtIKRzLOsNzX2zmg59207R2AHF3d+f6jo2dLkscVObwt9bmGWMeBpbguo1zurU2yRgzCVhrrY0HngDeNcY8huvi7xhrrY7sRcqZtZZ5P+9h6sJNHDudy4O/aclfrosgyF8f8fF1bvkLKLxnf9F568YXeZwM9HVHXyJSOqmHMhk3P5FV2zPo2rwOs4d3okPTWhf/RvEJevkX8TLZufm89c123vpmO/5+VZgyLJI7e4ZqEDY5h8JfxIt8v+0w4+YnsuPwKYZ2bUrsje1pGKxZteS/KfxFvMChkzlMXZjM/A17CatXg9n39+Q34Q2cLksqMYW/iAcrKLDMXbObZxdvIju3gD9fG84f+7cmwE+DsMmFKfxFPNSmfSeInZfA+l3H6N2qLlOGdaJNw5pOlyUeQuEv4mGyzuTxylcpTPtuB7UD/Xjpd10Y3q2ZBmGTS6LPc4tcijlzICwMqlRxfZ0zp0K7/yr5AANeWsE7K1K5vXsIXz9xNbdEhSj45ZLpyF+ktObMgZgYyMpyLe/c6VoGGDWqXLvee+w0T3+WxJKkA0Q0qsnHv+9Dj7C65dqneDdTWT9oGx0dbdeuXet0GSL/X1iYK/DP16IFpKWVS5d5+QXMWJXGy19uJd9aHr02gvuvakn1anrTLsUzxqyz1kZfrJ2O/EVKa9euS1tfRht2H2Pspwkk7zvBNW0bMGloJM3r1iiXvsT3KPxFSis0tPgjfzcPP34iO5fnv9jCv3/cScNgf94cFcUNkY11Xl/cSuEvUlpTp557zh+gRg3Xejew1vLZL/uY/HkyGZk53NMnjCeujyA4QLNqifsp/EVK69eLurGxrlM9oaGu4HfDxd6dGacYNz+RlSmH6dSsNtPv6UGnkNplfl6Rkij8RS7FqFFuvbMnJy+fuG9TeX35NvyqVmHizR24u08YVTUIm5Qzhb+IQ1anZhA7L4Hth05xY6cmPHVTBxrX1iBsUjEU/iIV7MipM0xduIn/W59OyBWBvD+mB9e0a+h0WeJjFP4iFcRay8dr0/nH4k1kZufxx/6teeS34QRW1yBsUvEU/iIVIOXASWLnJfJT2hF6hF3B1OGdiGgU7HRZ4sPcEv7GmEHAK7jm8J1mrX22mDa/AybimsN3o7X2/EneRbzO6TP5vL48hbgVqQT5V+Oft3bmtu4hmlVLHFfm8DfGVAXeAAYA6cAaY0x84by9v7YJB/4X6GutPWqM0QlO8XrfbDnIUwsS2X3kNLdGhTB2cDvq1fR3uiwRwD1H/j2BbdbaVABjzFxgKJBcpM2DwBvW2qMA1tqDbuhXpFI6cCKbSZ8ns/CXfbRuEMQHD/amT+t6Tpclcg53hH8zYHeR5XSg13ltIgCMMd/jOjU00Vr7hRv6Fqk08gss/169kxeWbCEnv4AnBkQQc3Ur/Kvpgq5UPhV1wbcaEA70B0KAFcaYTtbaY0UbGWNigBiAUDePlyJSnhL3HGfsvAR+ST/Ob8LrM3loJGH1g5wuS6RE7gj/PUDzIsshheuKSgd+tNbmAjuMMVtxvRisKdrIWhsHxIFrSGc31CZSrjJz8nhx6RZmrkqjbpA/r47sxs2dm2gQNqn03BH+a4BwY0xLXKE/Ajj/Tp75wEjgfWNMfVyngVLd0LeII6y1fJG4n6c/S+bAyWxG9QrlyYHtqB2oQdjEM5Q5/K21ecaYh4EluM7nT7fWJhljJgFrrbXxhduuN8YkA/nAk9bajLL2LeKE3UeymBCfxNebD9K+SS3euiuKbqFXOF2WyCXRTF4ipZSbX8C0lTt4ZdlWqhjD4wMiGHNlGNWqalYtqTw0k5eIG61NO0LsvES2HDjJ9R0aMWFIR5rVCXS6LJHLpvAXuYBjWWd4dvFm5q7ZTdPaAbw7OpoBHRo5XZZImSn8RYphrWXez3uYunATx07nEtOvFY9eG06Qv/7LiHfQX7LIebYfymTcvER+SM2gW2gdZg/rRIemtZwuS8StFP4ihbJz83nzm+28/c12AvyqMHV4JCN7hGoQNvFKCn8R4LuUw4ybn0BaRhZDuzZl3I0daBCsQdjEeyn8xacdOpnDlIXJLNiwl7B6Nfj3/b24Kry+02WJlDuFv/ikggLLB2t28ezizeTkFvDna8P5Y//WBPhpEDbxDQp/8Tmb9p1g7LwEft51jD6t6jFleCStG9R0uiyRCqXwF5+RdSaPf32Vwnvf7aBOoB8v/a4Lw7s10yBs4pMU/uITvkw+wMT4JPYcO83Ins35n0HtqFOjutNliThG4S9ebe+x00yMT2Jp8gHaNgrmk9/3ITqsrtNliThO4S9eKS+/gBmr0njpy60UWMv/DGrHA79piZ8GYRMBFP7ihTbsPsbYTxNI3neCa9o2YNLQSJrXreF0WSKVisJfvMbx07k8v2Qzc37cRcNgf94aFcWgyMa6oCtSDIW/eDxrLZ/9so/JnyeTkZnDmCvDeHxABMEBmlVLpCQKf/FoaYdP8dSCRFamHKZzSG2m39ODTiG1nS5LpNJT+ItHysnLJ+7bVF5bvo3qVavw9JCO3NW7BVU1CJtIqbgl/I0xg4BXcM3hO81a+2wJ7W4FPgF6WGs1R6Nclh+2ZzBufgLbD53ixk5NGH9zBxrVCnC6LBGPUubwN8ZUBd4ABgDpwBpjTLy1Nvm8dsHAo8CPZe1TfFNGZg7/WLSZ/1ufTvO6gbx/bw+uadvQ6bJEPJI7jvx7AtustakAxpi5wFAg+bx2k4HngCfd0Kf4kIICy8frdvPM4s1kZufxx/6teeS34QRW1yBsIpfLHeHfDNhdZDkd6FW0gTEmCmhurV1ojFH4S6ltPXCS2HkJrEk7Ss+wukwZHklEo2CnyxLxeOV+wdcYUwV4CRhTirYxQAxAaGho+RYmldrpM/m8+nUK765IpWZANf55a2du6x6iWbVE3MQd4b8HaF5kOaRw3a+CgUjgm8IP2zQG4o0xQ86/6GutjQPiAKKjo60bahMPtHzLQcYvSGT3kdPc1j2EsYPbUzdIg7CJuJM7wn8NEG6MaYkr9EcAd/660Vp7HDg7NZIx5hvgr7rbR8534EQ2T3+WxKKE/bRuEMTcmN70blXP6bJEvFKZw99am2eMeRhYgutWz+nW2iRjzCRgrbU2vqx9iHfLL7DM/iGNF5ZuJTe/gL9eH8GD/VrhX00XdEXKi1vO+VtrFwGLzls3voS2/d3Rp3iHxD3HGTsvgV/Sj/Ob8PpMGRZJi3pBTpcl4vX0CV9xxMnsXF5cupVZP6RRN8ifV0d24+bOTTQIm0gFUfhLhbLW8kXifiZ+lsTBkznc1asFfx3YltqBGoRNpCIp/KXC7D6SxfgFiSzfcogOTWrx9l3d6RZ6hdNlifgkhb+Uu9z8Aqat3MEry7ZSxRjG3dieMVeGUU2zaok4RuEv5Wpt2hHGzktg64FMru/QiIlDOtK0TqDTZYn4PIW/lItjWWd4dvFm5q7ZTbM6gbw7OpoBHRo5XZaIFFL4i1tZa/l0/R6mLtrE8dO5xPRrxaPXhhPkrz81kcpE/yPFbbYfymTcvER+SM2gW2gd/jG8E+2b1HK6LBEphsJfyiw7N583l2/j7W9TCfCrwtThkYzsEapB2EQqMYW/lMnKlEM8NT+RtIwshnVtSuyNHWgQ7O90WSJyEQp/uSwHT2Yz5fNNxG/cS8v6Qfz7/l5cFV7/4t8oIpWCwl8uSUGB5T8/7eK5LzaTk1vAo9eG84f+rQnw0yBsIp5E4S+llrz3BGPnJbBh9zGubF2PycMiad2gptNlichlUPjLRZ3KyeNfX21l+vdp1An04+U7ujCsazMNwibiwRT+ckFLk/YzMT6JvcezGdmzOf8zqB11amhWLRFPp/CXYu09dpoJ8Ul8mXyAto2C+WRkN6LD6jpdloi4icJfzpGXX8CMVWm89OVWCqzl7ze04/6rWuKnQdhEvIrCX876eddRxs5LZNO+E/y2XUOeHtKR5nVrOF2WiJQDt4S/MWYQ8AquOXynWWufPW/748ADQB5wCLjPWrvTHX1L2R0/ncvzSzYz58ddNAoO4O27ohjYsbEu6Ip4sTKHvzGmKvAGMABIB9YYY+KttclFmv0MRFtrs4wxfwD+CdxR1r6lbKy1xG/cy+TPN3HkVA5jrgzjievbUlODsIl4PXf8L+8JbLPWpgIYY+YCQ4Gz4W+tXV6k/WrgLjf0K2WQdvgUTy1IZGXKYTqH1GbGvT2IbFbb6bJEpIK4I/ybAbuLLKcDvS7Q/n5gsRv6lcuQk5fPO9+m8vrybVSvWoWnh3Tkrt4tqKpB2ER8SoW+vzfG3AVEA1eXsD0GiAEIDQ2twMp8ww/bM4idn0DqoVPc2LkJ42/qQKNaAU6XJSIOcEf47wGaF1kOKVx3DmPMdUAscLW1Nqe4J7LWxgFxANHR0dYNtQmQkZnD1EWb+HT9HprXDWTGvT3o37ah02WJiIPcEf5rgHBjTEtcoT8CuLNoA2NMN+AdYJC19qAb+pRSKCiwfLxuN88s3sypnDz+dE1rHr4mnMDqGoRNxNeVOfyttXnGmIeBJbhu9ZxurU0yxkwC1lpr44HngZrAx4W3D+6y1g4pa99Ssq0HThI7L4E1aUfpGVaXqcMjCW8U7HRZIlJJuOWcv7V2EbDovHXjizy+zh39yMWdPpPPq1+n8O6KVIIDqvHP2zpze/cQ3bMvIufQDd1eZPnmgzy1IJH0o6e5rXsIYwe3p26QBmETkf+m8PcC+49nM+nzJBYl7KdNw5rMjelN71b1nC5LRCoxhb8Hyy+wzPohjReXbiU3v4C/Xh9BTL/WVK+mQdhE5MIU/h7ql/RjxM5LJGHPcX4TXp8pwyJpUS/I6bJExEMo/D3MyexcXly6lVk/pFGvpj+vjezGTZ2b6IKuiFwShb+HsNayOHE/T3+WxMGTOdzVqwV/HdiW2oF+TpcmIh5I4e8Bdh/JYvyCRJZvOUSHJrV45+5oujav43RZIuLBFP6VWG5+Ae+uTOXVZSlUMYZxN7ZnzJVhVNPgzg+EAAAK0ElEQVSsWiJSRgr/SmpN2hFi5yWw9UAmAzs2YsLNHWlaJ9DpskTESyj8K5mjp87w7OLNfLh2N83qBDJtdDTXdWjkdFki4mUU/pWEtZZP1+9h6qJNHD+dy0P9WvHodeHUqK5fkYi4n5KlEth2MJNx8xNYnXqEqNA6TB3eifZNajldloh4MYW/g7Jz83lz+Tbe+nY7gX5V+cfwTozo0ZwqmlVLRMqZwt8hK1MOMW5+IjszshjerRljB7enQbC/02WJiI9Q+FewgyezmfL5JuI37qVl/SDmPNCLvm3qO12WiPgYhX8FKSiwzPlpF//8YjM5uQU8em04f+jfmgA/zaolIhVP4V8BkvYeJ3ZeIht2H+PK1vWYMiySVg1qOl2WiPgwhX85OpWTx8tfbuX9VWnUCfTj5Tu6MKxrMw3CJiKOc0v4G2MGAa/gmsN3mrX22fO2+wOzgO5ABnCHtTbNHX1XVkuS9jMxPol9x7MZ2TOUvw9qR+0aGoRNRCqHMoe/MaYq8AYwAEgH1hhj4q21yUWa3Q8ctda2McaMAJ4D7ihr35XRnmOnmbAgia82HaBto2Bev7Mb3VvUdbosEZFzuOPIvyewzVqbCmCMmQsMBYqG/1BgYuHjT4DXjTHGWmvd0H+lkJtfwPvf7+DlL1OwWP5+Qzvuv6olfhqETUQqIXeEfzNgd5HldKBXSW2stXnGmONAPeCwG/p33PpdRxn7aQKb95/k2nYNmTikI83r1nC6LBGRElWqC77GmBggBiA0NNThai7ueFYu/1yymf/8tItGwQG8fVcUAzs21gVdEan03BH+e4DmRZZDCtcV1ybdGFMNqI3rwu85rLVxQBxAdHR0pT0lZK0lfuNeJn+ezJFTZ7j3ypY8fn0ENf0r1WupiEiJ3JFWa4BwY0xLXCE/ArjzvDbxwD3AD8BtwNeeer4/7fApnlqQyMqUw3QOqc2Me3sS2ay202WJiFySMod/4Tn8h4EluG71nG6tTTLGTALWWmvjgfeA2caYbcARXC8QHiUnL593vk3l9eXb8K9ahUlDOzKqVwuqahA2EfFAbjlPYa1dBCw6b934Io+zgdvd0ZcTVm0/zLj5iaQeOsVNnZvw1E0daFQrwOmyREQum05SX0BGZg5TF23i0/V7CK1bgxn39qB/24ZOlyUiUmYK/2IUFFg+WrubZxZvJutMHg9f04aHf9tGg7CJiNdQ+J9ny/6TxM5LYO3Oo/RsWZepwyIJbxTsdFkiIm6l8C90+kw+ryxLYdrKVIIDqvHP2zpze/cQ3bMvIl5J4Q98vfkA4xckkX70NLd3D+F/B7enblB1p8sSESk3Ph3++49n8/RnSSxO3E+bhjX5MKY3vVrVc7osEZFy55Phn19gmbkqjReXbiGvwPLkwLY8+JtWVK+mQdhExDf4XPj/kn6MsfMSSNxzgn4RDZg8tCMt6gU5XZaISIXymfA/kZ3Li0u2MGv1TurX9Oe1kd24qXMTXdAVEZ/k9eFvrWVRwn6e/iyJQ5k53N27BX8d2JZaAZpVS0R8l1eH/66MLMbHJ/LNlkN0aFKLuNHRdG1ex+myREQc55XhfyavgHdXpvLqshSqVTE8dVMH7unTgmqaVUtEBPDC8N99JIv7Zqwh5WAmgzo2ZsKQDjSpHeh0WSIilYrXhX/DWv40qRPI329ox7XtGzldjohIpeR14e9frSqz7uvpdBkiIpWaToKLiPgghb+IiA9S+IuI+KAyhb8xpq4x5ktjTErh1yuKadPVGPODMSbJGPOLMeaOsvQpIiJlV9Yj/78Dy6y14cCywuXzZQGjrbUdgUHAv4wx5fdJqzlzICwMqlRxfZ0zp9y6EhHxVGUN/6HAzMLHM4Fh5zew1m611qYUPt4LHAQalLHf4s2ZAzExsHMnWOv6GhOjFwARkfOUNfwbWWv3FT7eD1zwxnpjTE+gOrC9jP0WLzYWsrLOXZeV5VovIiJnXfQ+f2PMV0DjYjadk6jWWmuMsRd4nibAbOAea21BCW1igBiA0NDQi5X233bturT1IiI+6qLhb629rqRtxpgDxpgm1tp9heF+sIR2tYCFQKy1dvUF+ooD4gCio6NLfCEpUWio61RPcetFROSssp72iQfuKXx8D7Dg/AbGmOrAPGCWtfaTMvZ3YVOnQo0a566rUcO1XkREzipr+D8LDDDGpADXFS5jjIk2xkwrbPM7oB8wxhizofBf1zL2W7xRoyAuDlq0AGNcX+PiXOtFROQsY+2ln12pCNHR0Xbt2rVOlyEi4lGMMeustdEXa6dP+IqI+CCFv4iID1L4i4j4IIW/iIgPUviLiPigSnu3jzHmEFDMJ7ZKrT5w2E3leApf22df21/QPvuKsuxzC2vtRcdPq7ThX1bGmLWlud3Jm/jaPvva/oL22VdUxD7rtI+IiA9S+IuI+CBvDv84pwtwgK/ts6/tL2iffUW577PXnvMXEZGSefORv4iIlMCjw98YM8gYs8UYs80Y81/zBxtj/I0xHxZu/9EYE1bxVbpXKfb5cWNMsjHmF2PMMmNMCyfqdKeL7XORdrcaY6wxxuPvDCnNPhtjflf4u04yxvynomt0t1L8bYcaY5YbY34u/Pse7ESd7mKMmW6MOWiMSSxhuzHGvFr48/jFGBPl1gKstR75D6iKazrIVrimhtwIdDivzR+BtwsfjwA+dLruCtjna4AahY//4Av7XNguGFgBrAaina67An7P4cDPwBWFyw2drrsC9jkO+EPh4w5AmtN1l3Gf+wFRQGIJ2wcDiwED9AZ+dGf/nnzk3xPYZq1NtdaeAebimlC+qKITzH8CXGuMMRVYo7tddJ+ttcuttb9OZLwaCKngGt2tNL9ngMnAc0B2RRZXTkqzzw8Cb1hrjwJYa4udRc+DlGafLVCr8HFtYG8F1ud21toVwJELNBmKaxIsa10zINYpnDHRLTw5/JsBu4sspxeuK7aNtTYPOA7Uq5Dqykdp9rmo+3EdOXiyi+5z4dvh5tbahRVZWDkqze85AogwxnxvjFltjBlUYdWVj9Ls80TgLmNMOrAIeKRiSnPMpf5/vyQXncNXPJMx5i4gGrja6VrKkzGmCvASMMbhUipaNVynfvrjene3whjTyVp7zNGqytdIYIa19kVjTB9gtjEm0lpb4HRhnsiTj/z3AM2LLIcUriu2jTGmGq63ihkVUl35KM0+Y4y5DogFhlhrcyqotvJysX0OBiKBb4wxabjOjcZ7+EXf0vye04F4a22utXYHsBXXi4GnKs0+3w98BGCt/QEIwDUGjrcq1f/3y+XJ4b8GCDfGtCycJH4Ergnliyo6wfxtwNe28EqKh7roPhtjugHv4Ap+Tz8PDBfZZ2vtcWttfWttmLU2DNd1jiHWWk+eA7Q0f9vzcR31Y4ypj+s0UGpFFulmpdnnXcC1AMaY9rjC/1CFVlmx4oHRhXf99AaOW2v3uevJPfa0j7U2zxjzMLAE150C0621ScaYScBaa2088B6ut4bbcF1YGeFcxWVXyn1+HqgJfFx4bXuXtXaIY0WXUSn32auUcp+XANcbY5KBfOBJa63Hvqst5T4/AbxrjHkM18XfMZ58MGeM+QDXC3j9wusYEwA/AGvt27iuawwGtgFZwL1u7d+Df3YiInKZPPm0j4iIXCaFv4iID1L4i4j4IIW/iIgPUviLiPgghb+IiA9S+IuI+CCFv4iID/p/6M9y05/bFyEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[0,   1],\n",
    "              [0.5, 1],\n",
    "              [1,   1]])\n",
    "B = np.array([[-0.2],\n",
    "              [0.7],\n",
    "              [1.2]])\n",
    "\n",
    "X = lsqr(A, B)[0]\n",
    "\n",
    "plt.plot([0, 0.5, 1], [-0.2, 0.7, 1.2], 'ro', label='data')\n",
    "plt.plot([0, 1], [X[0]*0+X[1], X[0]*1+X[1]], label='regression')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Решение в лоб для идиотов\n",
    "Рассмотрим задачу: машинка едет с некой начальной скоростью $x_s$ и ускоряется до некой заданной $x_e$. Описание системы:\n",
    "$$\n",
    "x_{k+1} = x_k + u_k\n",
    "$$\n",
    "$\n",
    "\\begin{align*}\n",
    "    \\text{где } & x - \\text{скорость,}\\\\\n",
    "                & u - \\text{ускорение.}\n",
    "\\end{align*}\n",
    "$\n",
    "\n",
    "Решим с помощью метода наименьших квадратов. Входными данными будет скорость, выходными - ускорение. Но, в отличие от простейшего примера с прямой, у нас нет просто матрицы входов и вектора выходов, потому что каждый следующий выход $x_{k+1}$ зависит от ускорения $u_k$. Поэтому будем считать состояние системы (машинки) итеративно (пошагово).\n",
    "\n",
    "1. Итеративно промоделируем систему и найдем входные ($x$) и выходные ($u$) данные\n",
    "2. Используя полученный на предыдущем шаге датасет, применем метод наименьших кваратов, найдем искомый коэффициент $F$ для уравнения: $u=-Fx$\n",
    "\n",
    "### 1. Итеративно моделируем\n",
    "\n",
    "При этом должны соблюдаться ограничения:\n",
    " - скорость и ускорение должны удовлетворять уравнению системы\n",
    " - начальная скорость $x_s$\n",
    " - конечная скорость $x_e$\n",
    "\n",
    "Составим систему уравнений, описывающую машинку пошагово:\n",
    "$$\n",
    "\\begin{cases}\n",
    "    x_1 = x_0 + u_0 \\\\\n",
    "    x_2 = x_1 + u_1 \\\\\n",
    "    ... \\\\\n",
    "    x_n = x_{n-1} + u_{n-1}\n",
    "\\end{cases}\n",
    "$$\n",
    "Добавим в нее ограничения на начальную и конечную скорость:\n",
    "$$\n",
    "\\begin{cases}\n",
    "    x_1 = x_0 + u_0 \\\\\n",
    "    x_2 = x_1 + u_1 \\\\\n",
    "    ... \\\\\n",
    "    x_n = x_{n-1} + u_{n-1} \\\\\n",
    "    -------\\\\\n",
    "    x_0 = x_s \\\\\n",
    "    x_n = x_e\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "В таком виде систему не представить в необходимой нам матиричной форме $AX=B$. Перепишем ее:\n",
    "$$\n",
    "\\begin{cases}\n",
    "    1 x_0 - 1 x_1 + 0 x_2 + ... + 0 x_{n-1} + 0 x_n + 1 u_0 + 0 u_1 + ... + 0 u_n = 0 \\\\\n",
    "    0 x_0 + 1 x_1 - 1 x_2 + ... + 0 x_{n-1} + 0 x_n + 0 u_0 + 1 u_1 + ... + 0 u_n = 0 \\\\\n",
    "    ... \\\\\n",
    "    0 x_0 + 0 x_1 + 0 x_2 + ... + 1 x_{n-1} - 1 x_n + 0 u_0 + 0 u_0 + ... + 1 u_n = 0 \\\\\n",
    "    ---------------------------\\\\\n",
    "    1 x_0 + 0 x_1 + 0 x_2 + ... + 0 x_{n-1} + 0 x_n + 0 u_0 + 0 u_1 + ... + 0 u_n = x_b \\\\\n",
    "    0 x_0 + 0 x_1 + 0 x_2 + ... + 0 x_{n-1} + 1 x_n + 0 u_0 + 0 u_1 + ... + 0 u_n = x_e \\\\\n",
    "\\end{cases}\n",
    "$$\n",
    "В матричном виде:\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "1& -1&  0& ...& 0&  0& 1& 0& ...& 0& \\\\\n",
    "0&  1& -1& ...& 0&  0& 0& 1& ...& 0& \\\\\n",
    "... \\\\\n",
    "0&  0&  0& ...& 1& -1& 0& 0& ...& 1& \\\\\n",
    "1&  0&  0& ...& 0&  0& 0& 0& ...& 0& \\\\\n",
    "0&  0&  0& ...& 1&  0& 0& 0& ...& 0&\n",
    "\\end{bmatrix}\n",
    "\\cdot\n",
    "\\begin{bmatrix}\n",
    "x_0 \\\\\n",
    "x_1 \\\\\n",
    "x_2 \\\\\n",
    "... \\\\\n",
    "x_{n-1} \\\\\n",
    "x_n \\\\\n",
    "u_0 \\\\\n",
    "u_1 \\\\\n",
    "... \\\\\n",
    "u_n \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "... \\\\\n",
    "0 \\\\\n",
    "x_b \\\\\n",
    "x_e \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "$$\n",
    "AX=B\n",
    "$$\n",
    "\n",
    "Получилась вот такая странная и разреженная система.\n",
    "\n",
    "Размерность матрицы получается:\n",
    " - n уравнений шагов + 2 уравнения ограничений = $n + 2$ уравнений\n",
    " - $x_0$ ... $x_n$ $\\rightarrow$ n+1 переменная, $u_0$ ... $u_{n-1}$ $\\rightarrow$ n переменная, итого $2n+1$ переменных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Упрощает запись сложных систем для МНК\n",
    "class LSQR_Helper:\n",
    "    \n",
    "    # Инициализация\n",
    "    # h - высота матрицы (кол-во уравнений)\n",
    "    # w - ширина матрицы (кол-во переменных)\n",
    "    def __init__(self, h, w):\n",
    "        self.a = csc_matrix((h, w))\n",
    "        self.b = np.zeros((h,1))\n",
    "        self._row = 0\n",
    "    \n",
    "    # Добавляет уравнение в систему\n",
    "    # a - строка матрицы\n",
    "    # b - свободный член\n",
    "    def add_row(self, a, b):\n",
    "        for i in range(0, len(a)):\n",
    "            self.a[self._row, i] = a[i]\n",
    "        \n",
    "        self.b[self._row] = b\n",
    "        self._row+=1\n",
    "        \n",
    "    # Записывает коэффициент в матрицу в текущую строку (уравнение)\n",
    "    # i  - номер столбца (переменной)\n",
    "    # ai - значение\n",
    "    def set_a(self, i, ai):\n",
    "        self.a[self._row, i] = ai\n",
    "        \n",
    "    # Записывает свободный член в текущую строку (уравнение)\n",
    "    # b - значение\n",
    "    def set_b(self, b):\n",
    "        self.b[self._row] = b\n",
    "        \n",
    "    # Начинает новую строку уравнения\n",
    "    def begin_new_row(self):\n",
    "        self._row+=1\n",
    "        \n",
    "# Решает систему и рисует графики скорости и ускорения\n",
    "def solve_and_plot(eq):\n",
    "    sol = lsqr(eq.a, eq.b)[0]\n",
    "    v = sol[0:n+1]\n",
    "    u = sol[n+1:]\n",
    "\n",
    "    plt.plot(v, label='v')\n",
    "    plt.plot(u, label='u')\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Решим получившуюся систему с помощью метода наименьших квадратов (`scipy.sprase.linalg.lsqr`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/scipy/sparse/compressed.py:746: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.\n",
      "  SparseEfficiencyWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt83HWd7/HXJ/febwnpJYUWaIG2XFpLRRFFQWg5So+3BZRzUFk4exSRh57dR4+uleXsHmRFz8PVri4oirqKXdS1YktBrt6ABlpoegFKKXSapEnTNL3knvmcP2YmTJKZZNrO5TeT9/PxyCMzv983M5/HL9N3f/l+v7/vz9wdEREpLEW5LkBERNJP4S4iUoAU7iIiBUjhLiJSgBTuIiIFSOEuIlKAFO4iIgVI4S4iUoAU7iIiBagkV29cWVnpc+bMydXbi4jkpeeff/6Au1eN1C5n4T5nzhxqa2tz9fYiInnJzN5IpZ26ZURECpDCXUSkACncRUQKUM763BPp6ekhFArR2dmZ61KSqqiooKamhtLS0lyXIiKSVKDCPRQKMWHCBObMmYOZ5bqcIdydlpYWQqEQc+fOzXU5IiJJjdgtY2b3mVmTmdUl2W9m9i9mtsvMXjKzJSdaTGdnJ9OmTQtksAOYGdOmTQv0XxYiIpBan/uPgOXD7F8BzIt+3Qx892QKCmqwxwS9PhERSCHc3f1p4OAwTVYCP/aIZ4DJZjYjXQWKiBSK7t4wd27YwYt7D2X8vdIxW2YWsDfueSi6bQgzu9nMas2strm5OQ1vLSKSPxrbOvm3p3bz8v4jGX+vrE6FdPd73H2puy+tqhrx6lkRkYISam0HoGbKmIy/VzrCfR8wO+55TXRb3lm1ahVr1qzpf3777bdz991357AiESkkodYOAGZPGZvx90rHVMh1wC1m9gDwdqDN3RtO9kX/4bfb2F5/+KSLi7dg5kS++sGFSfdfc8013HbbbXz2s58FYO3atWzcuDGtNYjI6BVqbafIYPqkioy/14jhbmY/By4FKs0sBHwVKAVw9+8B64GrgF1AO/CpTBWbaYsXL6apqYn6+nqam5uZMmUKs2fPHvkHRURSEGrtYMakMZQWZ75HfMRwd/frRtjvwGfTVlHUcGfYmfSxj32MBx98kMbGRq655pqc1CAihSnU2sGsLPS3Q8CuUA2Ca665hptuuokDBw7w1FNP5bocESkgodZ2LjpjWlbeSwuHDbJw4UKOHDnCrFmzmDFD0/VFJD26e8M0Hu6kJguDqaAz94S2bt2a6xJEpMA0tnUSdqiZnJ1uGZ25i4hkQehQ9ua4g8JdRCQrYnPcs9Uto3AXEcmCUGtH1ua4g8JdRCQrQq3tTJ9YQVlJdmJX4S4ikgWh1o6sdcmAwl1EJCv2tXZkbTAVFO4iIhnX0xemoU3hLiJSUGJz3LO19AAo3IfYs2cPixYt6n9+9913c/vtt+euIBHJe3v713HPXp97cK9Q3bAKGtN8pej0c2HF19L7miIiI9jXP8ddZ+4iIgUj1NqBGcyYlL1wD+6Ze47OsEtKSgiHw/3POzs7c1KHiBSOUGtHVue4g87ch6iurqapqYmWlha6urp46KGHcl2SiOS5UGt7VrtkIMhn7jlSWlrK6tWrWbZsGbNmzeLss8/OdUkikudCrR0smzs1q++pcE/g1ltv5dZbb811GSJSAHr7Yuu4Z/fMXd0yIiIZ1NDWSV/YFe4iIoUkttTvrMnZm+MOAQz3yP22gyvo9YlIsOw7lP057hCwcK+oqKClpSWwAerutLS0UFGRnfWYRST/hVrbI3PcJ2c3NwI1oFpTU0MoFKK5uTnXpSRVUVFBTU1NrssQkTwRau2gekIF5SXFWX3fQIV7aWkpc+fOzXUZIiJpk4s57hCwbhkRkUITyvI67jEKdxGRDOntC9PQ1pnV1SBjFO4iIhnSeDg3c9xB4S4ikjGh/qV+deYuIlIw+i9g0pm7iEjhiN2kY2aW57hDiuFuZsvN7GUz22VmqxLsP9XMnjCzzWb2kpldlf5SRUTyS6i1neqJ5Vmf4w4phLuZFQNrgBXAAuA6M1swqNnfA2vdfTFwLfCv6S5URCTfRKZBZr+/HVI7c18G7HL33e7eDTwArBzUxoGJ0ceTgPr0lSgikp9Ch3JzAROkFu6zgL1xz0PRbfFuB643sxCwHvhcWqoTEclTvX1hGg5lfx33mHQNqF4H/Mjda4CrgJ+Y2ZDXNrObzazWzGqDvH6MiMjJ2n+ki96wB7pbZh8wO+55TXRbvBuBtQDu/hegAqgc/ELufo+7L3X3pVVVVSdWsYhIHggdbAeyv9RvTCrhvgmYZ2ZzzayMyIDpukFt3gQuAzCzc4iEu07NRWTUyuUFTJBCuLt7L3ALsBHYQWRWzDYzu8PMro42+yJwk5m9CPwc+KQHdVF2EZEsiIX7jEm5uf9DSkv+uvt6IgOl8dtWxz3eDlyc3tJERPLXvkPtnDKhnIrS7M9xB12hKiKSEbla6jdG4S4ikgG5vIAJFO4iImnXF3bqD+nMXUSkoOw/3JnTOe6gcBcRSbu3pkHqzF1EpGCEWnN7ARMo3EVE0i7Uv467wl1EpGCEWtupyuEcd1C4i4ik3b4cz5QBhbuISNrleo47KNxFRNIqCHPcQeEuIpJWTUc66elzhbuISCHJ9VK/MQp3EZE0CsIcd1C4i4ikVehg5Mx9Vg7nuIPCXUQkrUKtHVSOz+0cd1C4i4ikVehQe867ZEDhLiKSVvtyfJOOGIW7iEiahMMevTo1tzNlQOEuIpI2TUe6AjHHHRTuIiJpE5RpkKBwFxFJm6BcwAQKdxGRtNGZu4hIAQrKHHdQuIuIpE2otYNZAThrB4W7iEjahFqDcQETKNxFRNIiHHbqD3Uq3EVECknz0S66+8KBmCkDCncRkbQI0kwZULiLiByf3m7o6RyyOTbHfbbCXUQkD73+NNw5C0LPD9gcC/dZk/OoW8bMlpvZy2a2y8xWJWnzV2a23cy2mdnP0lumiEhANGyGcC9Unjlgc6i1ncrxZYwpy/0cd4CSkRqYWTGwBng/EAI2mdk6d98e12Ye8L+Bi9291cxOyVTBIiI5Vb8Fpp0JFZMGbA61duT87kvxUjlzXwbscvfd7t4NPACsHNTmJmCNu7cCuHtTessUEQmI+i0wc/GQzaHWYCz1G5NKuM8C9sY9D0W3xZsPzDezP5nZM2a2PNELmdnNZlZrZrXNzc0nVrGISK4cbYLDIZhxwYDN4bAH5iYdMekaUC0B5gGXAtcB95rZ5MGN3P0ed1/q7kurqqrS9NYiIllSvyXyfdCZ+4H+Oe75Fe77gNlxz2ui2+KFgHXu3uPurwOvEAl7EZHC0bAFMJhx3oDNewO01G9MKuG+CZhnZnPNrAy4Flg3qM1/Ejlrx8wqiXTT7E5jnSIiuVe/GSrnQfmEAZuDdgETpBDu7t4L3AJsBHYAa919m5ndYWZXR5ttBFrMbDvwBPC37t6SqaJFRHJimMFUIDArQkIKUyEB3H09sH7QttVxjx34QvRLRKTwHNkPR+qHDKZCJNynjStjbFlKkZoVukJVRCQVDYkHUyFYS/3GKNxFRFJRvxkwmH7ukF37AnSTjhiFu4hIKuo3Q9VZUD5+wOZw2AkdCtYFTKBwFxFJTf2WhP3tB4510d0brDnuoHAXERnZ4QY42jjsTBmFu4hIvhl2MDV4FzCBwl1EZGT1m8GKEg6mxi5gCtKKkKBwFxEZWf1mqDobyoaenYdaO5g6roxx5cGZ4w4KdxGR4bknHUyF2FK/wTprB4W7iMjwDtfDsaaE/e0QzAuYQOEuIjK8YQZT3SPruAetvx0U7iIiw6vfDFYM0xcN2dV8tIuu3nDgZsqAwl1EZHj1m+GUc6B06Nn5voDOcQeFu4hIcikMpkLw5riDwl1EJLm2ELQfgJnDh3vQFg0DhbuISHL9g6lLEu4OtbYzZWwp4wM2xx0U7iIiydVvhqISqF6YcHdkjnvwumRA4S4iklz/YGpFwt1BneMOCncRkcRGGEx198BenQoKdxGRxA69CR0Hk16ZeuBoN1294UBewAQKdxGRxPoHU5PNlImsBqk+dxGRfFK/GYpKoXrolakA+w5F57hP1Zm7iEj+qN8M1QugpDzh7v457uqWERHJEyMMpkKkW2by2FImVJRmsbDUKdxFRAZr3QOdh5IOpkJw13GPUbiLiAw2wmAqRMN9cjAHU0HhLiIyVP1mKC6DUxYk3B2Z4x7cC5hA4S4iMlT95siSA0kGU1uOddPZEw7kgmExCncRkXjuUP/iCIOpwV3qN0bhLiIS7+Bu6GobYTA1dgFTnp+5m9lyM3vZzHaZ2aph2n3EzNzMlqavRBGRLEphMHVfgNdxjxkx3M2sGFgDrAAWANeZ2ZBRBjObAHweeDbdRYqIZE39Ziguh6pzkjYJtXYwaUwpEwM6xx1SO3NfBuxy993u3g08AKxM0O7/AHcBnWmsT0Qku+q3RG6GXVKWtEnQZ8pAauE+C9gb9zwU3dbPzJYAs939d2msTUQku8JhaBh+MBWCfwETpGFA1cyKgG8CX0yh7c1mVmtmtc3NzSf71iIi6XVwN3QdHnYw9a113IM7UwZSC/d9wOy45zXRbTETgEXAk2a2B7gIWJdoUNXd73H3pe6+tKqq6sSrFhHJhBQGUw8e66ajp68gztw3AfPMbK6ZlQHXAutiO929zd0r3X2Ou88BngGudvfajFQsIpIp9ZuhpAKqzk7aJOirQcaMGO7u3gvcAmwEdgBr3X2bmd1hZldnukARkayp3xJZv704+SyYnY2HAZhTOS5bVZ2QklQauft6YP2gbauTtL305MsSEcmycDjSLXP+dcM2e7iukZopY5h3yvgsFXZidIWqiAhAyy7oPjrsYGpbRw9/3HWAFYumY2ZZLO74KdxFRCClwdTHd+6np89Zce6MLBV14hTuIiIQHUwdA5VnJW2yfmsj0ydWcEHN5CwWdmIU7iIiEL0y9VwoTjwUebSrl6deaWb5oukUFQW7SwYU7iIiEO6LXJk6TH/7Ezub6O4Ns2LR9CwWduIU7iIiB16FnmPDhvvDdY1Uji9n6ZypWSzsxCncRURGGEzt6O7j8Z1NXLmwmuI86JIBhbuISGQwtXQsVM5PuPupV5rp6OnjqjyYJROjcBcRqd8C08+DouKEuzfUNTBlbClvn5sfXTKgcBeR0a6vFxpfStrf3tXbx2M7mrhiwXRKivMnMvOnUhGRTDjwCvS0Jw33P756gKNdvSw/Nz9mycQo3EVkdBthMHVDXSMTKkq4+IzKLBZ18hTuIjK61W+GsvEw7cwhu3r6wjy6fT/vP6easpL8isv8qlZEJN2GGUz9y2sttHX05MVaMoMp3EVk9BphMHVDXQPjyoq5ZF5+dcmAwl1ERrPmndDbmbC/vbcvzCPb9vO+c6qpKE08RTLIFO4iMnr1D6YOPXN/bs9BWo51581aMoMp3EVk9KrfDGUTYOoZQ3Y9XNdIRWkRl55VlYPCTp7CXURGr/otMON8KBoYheGw83BdI5fOP4WxZSndjTRwFO4iMjr19UDj1oT97S+82UrTkS5W5NmFS/EU7iIyOjXtgL6uhP3tG+oaKSsu4n1nn5KDwtJD4S4io1OSwVT3SJfMJfMqmVBRmoPC0kPhLiKjU/1mKJ8EU+YO2PxSqI19hzry8sKleAp3ERmd6rfAjPOGDKaur2ugpMh4/znVOSosPRTuIjL69HbD/rqkXTLvPLOSSWPzt0sGFO4iMho1bYe+7iEzZbY3HOaNlva8vXApnsJdREafJIOpD9c1UmRwxYL87pIBhbuIjDbu8OIDMGHmkMHU9VsbePvcaUwbX56j4tJH4S4io8urj8Kbf4F3fxHM3tq8/wivNR/jqjy+cCmewl1ERo9wGB67I3LGvuSGAbvWb23EDK5cWBjhnp+LJoiInIhtv4L9W+HD34figbNhNtQ1sPS0KZwysSJHxaVXSmfuZrbczF42s11mtirB/i+Y2XYze8nMHjOz09JfqojISejrgSf+CaoXwaKPDNj1+oFj7Gw8wvJF+X3hUrwRw93MioE1wApgAXCdmS0Y1GwzsNTdzwMeBP453YWKiJyUzT+Fg7vhfV8ZcuHShroGAJYXwBTImFTO3JcBu9x9t7t3Aw8AK+MbuPsT7t4effoMUJPeMkVETkJPBzx1F8x+O8y/csjuDVsbOX/2ZGZNHpOD4jIjlXCfBeyNex6KbkvmRmBDoh1mdrOZ1ZpZbXNzc+pVioicjOfugSMNcNlXB8yQAdh7sJ2t+9q4qoDO2iHNs2XM7HpgKfD1RPvd/R53X+ruS6uq8vPuJiKSZzrb4I//D868HOZcPGT3w3WNAKwooP52SG22zD5gdtzzmui2AczscuDLwHvcvSs95YmInKQ/fxs6WuGy1Ql3b6hrYOHMiZw6bWyWC8usVM7cNwHzzGyumZUB1wLr4huY2WLg34Cr3b0p/WWKiJyAo03wl3+FhR+K3E5vkMa2Tl5481BBrCUz2Ijh7u69wC3ARmAHsNbdt5nZHWZ2dbTZ14HxwH+Y2RYzW5fk5UREsucP34DeTnjv3yfc/XB0lky+r92eSEoXMbn7emD9oG2r4x5fnua6REROzqE3ofY+WPwJqDwzYZMNdY3Mrx7PGVXjs1xc5mn5AREpTE9+DTB4z5DrLgFoPtLFc3sOFtxAaozCXUQKT9NOePHnsOwmmJR45vYj2xtxhxUFslDYYAp3ESk8T/wjlI6Dd30haZMNWxs5vXIcZ1VPyGJh2aNwF5HCsu952PFbeOfnYNy0hE1aj3Xzl90tLF80HRt0UVOhULiLSGF57A4YWwnv+EzSJr/b2kBf2LmqAGfJxCjcRaRw7H4y8nXJF6E8cXfL9vrD3Ll+B+fPnszCmROzWl42KdxFpDC4R87aJ9bA0k8nbLL/cCc33r+JiWNKuee/va1gu2RA4S4ihWLnQ5H+9ktXQenQG260d/fy1/fXcrijhx/ccCHVBXJTjmR0JyYRyX/hPnj8H6FyPpx/3ZDdfWHntge2sK2+je/fsJQFBdwdE6NwF5H899IvoHknfOx+KB4aa3c9vJNHtu/n9g8u4H1nV+egwOxTt4yI5LfeLnjiTphxASxYOWT3z559k3ue3s0N7ziNT148NwcF5obO3EUkvz3/I2h7E67+1pAbcfzh1Wa+8ps63ntWFV/5wOC7gxY2nbmLSP7qOgpPfx3mXAKnv3fArlf2H+EzP32BeaeM59sfX0JJ8eiKO525i0j+eva7cKwZrv35gLP25iNdfOqHm6goK+a+T17I+PLRF3Wj678yESkc7QfhT9+Gs/4LzL6wf3NnTx83/biWlmNd/OCGpcwsoJteH4/R99+ZiOS/cB/8/nboOgzve+tGHOGw88W1L/Ji6BDfu/5tnFczOXc15pjCXUTyy+F6+NXNsOcPcNFnoPqtgdJvPPoyv9vawJeuOpsrFxbmUr6pUriLSP7YuR5+8xno7YaVa+CCT/TvWlu7lzVPvMZ1y07lpktOz2GRwaBwF5Hg6+mAR74Cm+6F6efBR3844NZ5f37tAF/61VYumVfJHSsXFvSaMalSuItIsDXtgAc/DU3b4R23wGWroaS8f/drzUf5nz99gbmV4/jOx5dQOsqmPCajcBeRYHKP3OB645ciy/d+4pcw7/IBTQ4e6+bTP9pESZFx3ycvZNKY0hwVGzwKdxEJnvaDsO5zkZUez7gMPvQ9GH/KgCZdvX3c/ONaGto6eeDmi5g9dWyOig0mhbuIBMueP8Ivb4pcnHTFP0VmxBS91dXi7jyz+yDfeuwVat9o5TsfX8ySU6fksOBgUriLSDD09cJTX4On74ZpZ8B1v4eZF/Tv7ukL89BL9Xz/D6+zrf4wU8eV8X8/dC4fOG9mDosOLoW7iORe6xvwy7+G0HNwwfWw4i4oHw9AW3sPP3vuTe7/8x4aD3dyRtU47vzwuXxo8SwqSotzXHhwKdxFJLfqfgm/vS3y+CM/gHM/CsCbLe3c96fXWVu7l/buPi4+cxp3fvhc3jO/iqIiTXUcicJdRLKvsw1efzoS7Nt+DTXL4CP34pNP44U3DnLv06/zyPZGiouMD54/kxvfNZeFMyfluuq8onAXkcwL90H9FnjtMdj1GIQ2gfdB2Xh499/Re8nfsnFHC/f+7M9s2XuIiRUl/I/3nMEN75jD9EmFfa/TTFG4i0hmtO2D1x6PBPruJ6GjFbDIIOm7boMzLuNI1WLWbt7PD7/5R0KtHZw2bSx3rFzIR5bUMG4ULtObTjp6IpIePR3wxp9g1+ORUG/eEdk+fjrhectprn4XdeWL2XqolFcbj/Lyi0d4/cCT9IWdC+dM4SsfWMDl51RTrP70tEgp3M1sOfAtoBj4vrt/bdD+cuDHwNuAFuAad9+T3lJFJBD6euBIIxxpiKzQ2LoHXn8K3vgz9HbixeUcrHwbL5/+ef4QPp8nD1Xx2gvH6O4NA3swg1OnjmV+9QSuXFjN+xdM54LZo3dp3kwZMdzNrBhYA7wfCAGbzGydu2+Pa3Yj0OruZ5rZtcBdwDWZKFhEMsQ9MtAZC+0jDXC4AY7U44frCR+ObC9qP4DhA350f/lpPFu6nIe6z+Hpzvl0Hous/TJzUgXzp1dwyfwq5ldP4KzqCZx5ynjGlGkKY6alcua+DNjl7rsBzOwBYCUQH+4rgdujjx8EvmNm5u4DPwEiMrJwGMK9kS/viz6OfA/39hDu6yUc7sX7egj39RHu68H7egn3dNHb1U5f1zHC3ccId7UT7m6HnmPR7x1Ydzv0dlDU24H1tFPU20FRXwelPUcZ29VEWbhzSDmtTKAxPIVGn0KjL2I/U2j0qTT6FPb7VBp8KiUllcyvHM/8cybw1ekTmF89nnnVE5hYobVeciWVcJ8F7I17HgLenqyNu/eaWRswDTiQjiLjbfrVt6iquzfdL1vQ8rcH88TPDQafWdow5xlD2sY9jz1+a5v3H8/4bZHnkfeJbS8ijOEDvopwiH4fuD0MQDHhaJvEijixe2P2udFBeeTLyzhGOR2U0eEVdFDGMWbSzCLaSis5WlpFR0U13eOm0zeumnHjxjFpTCkTK0qZNKaUs8aUcuGYyONJY0qZNLZUIR5AWR1QNbObgZsBTj311BN6jZLx0zg4dm46yxoV/CQiPpf/OZxM3Qxa03vgaw23b+BrxMd5f3sbvM3eeh2LxHfkZ4uibQ23oki76GvGP+/fbiWErRiKivGiyGOzyGMvKgErhqKS/v1Y5DHFJZHtZeOhdCxF5WMpjn4vKh9HaVkFpSXFlBUXUVZSRGmxMa6kqP95eUkxFaVFWge9gKQS7vuA2XHPa6LbErUJmVkJMInIwOoA7n4PcA/A0qVLT+i0bPEV18MV15/Ij4qIjBqp/IW3CZhnZnPNrAy4Flg3qM064Ibo448Cj6u/XUQkd0Y8c4/2od8CbCQyFfI+d99mZncAte6+DvgB8BMz2wUcJPIfgIiI5EhKfe7uvh5YP2jb6rjHncDH0luaiIicKN1sUESkACncRUQKkMJdRKQAKdxFRAqQwl1EpABZrqajm1kz8MYJ/nglGVjaIA1U1/FRXccvqLWpruNzMnWd5u5VIzXKWbifDDOrdfelua5jMNV1fFTX8Qtqbarr+GSjLnXLiIgUIIW7iEgBytdwvyfXBSShuo6P6jp+Qa1NdR2fjNeVl33uIiIyvHw9cxcRkWEEOtzNbLmZvWxmu8xsVYL95Wb2i+j+Z81sThZqmm1mT5jZdjPbZmafT9DmUjNrM7Mt0a/ViV4rA7XtMbOt0fesTbDfzOxfosfrJTNbkoWazoo7DlvM7LCZ3TaoTdaOl5ndZ2ZNZlYXt22qmT1qZq9Gv09J8rM3RNu8amY3JGqTxpq+bmY7o7+nX5tZwjtIj/Q7z1Btt5vZvrjf11VJfnbYf78ZqOsXcTXtMbMtSX42I8csWTbk7PPl7oH8IrK88GvA6UAZ8CKwYFCbzwDfiz6+FvhFFuqaASyJPp4AvJKgrkuBh3JwzPYAlcPsvwrYQOS2QRcBz+bgd9pIZJ5uTo4X8G5gCVAXt+2fgVXRx6uAuxL83FRgd/T7lOjjKRms6QqgJPr4rkQ1pfI7z1BttwP/K4Xf9bD/ftNd16D93wBWZ/OYJcuGXH2+gnzm3n9jbnfvBmI35o63Erg/+vhB4DLL8H3C3L3B3V+IPj4C7CByD9l8sBL4sUc8A0w2sxlZfP/LgNfc/UQvXjtp7v40kXsOxIv/HN0P/NcEP3ol8Ki7H3T3VuBRYHmmanL3R9y9N/r0GSJ3QMu6JMcrFan8+81IXdEM+Cvg5+l6vxRrSpYNOfl8BTncE92Ye3CIDrgxNxC7MXdWRLuBFgPPJtj9DjN70cw2mNnCLJXkwCNm9rxF7lc7WCrHNJOuJfk/uFwcr5hqd2+IPm4EqhO0yeWx+zSRv7gSGel3nim3RLuM7kvSzZDL43UJsN/dX02yP+PHbFA25OTzFeRwDzQzGw/8ErjN3Q8P2v0Cka6H84FvA/+ZpbLe5e5LgBXAZ83s3Vl63xFZ5BaNVwP/kWB3ro7XEB75GzkwU8jM7MtAL/DvSZrk4nf+XeAM4AKggUgXSJBcx/Bn7Rk9ZsNlQzY/X0EO9+O5MTc2zI25083MSon88v7d3X81eL+7H3b3o9HH64FSM6vMdF3uvi/6vQn4NZE/jeOlckwzZQXwgrvvH7wjV8crzv5Y91T0e1OCNlk/dmb2SeADwCeioTBECr/ztHP3/e7e5+5h4N4k75mTz1o0Bz4M/CJZm0wesyTZkJPPV5DDPZA35o725/0A2OHu30zSZnqs79/MlhE5zhn9T8fMxpnZhNhjIgNydYOarQP+u0VcBLTF/bmYaUnPpnJxvAaJ/xzdAPwmQZuNwBVmNiXaDXFFdFtGmNly4O+Aq929PUmbVH7nmagtfpzmQ0neM5V/v5lwObDT3UOJdmbymA2TDbn5fKV7xDidX0Rmd7xCZNT9y9FtdxD5wANUEPkzfxfwHHB6Fmp6F5E/q14CtkS/rgL+BvibaJtbgG1EZgg8A7wzC3WdHn2/F6PvHTte8XUZsCZ6PLcCS7P0exxHJKyh+CkAAAAAn0lEQVQnxW3LyfEi8h9MA9BDpF/zRiLjNI8BrwK/B6ZG2y4Fvh/3s5+OftZ2AZ/KcE27iPTBxj5jsVlhM4H1w/3Os3C8fhL9/LxEJLhmDK4t+nzIv99M1hXd/qPY5yqubVaO2TDZkJPPl65QFREpQEHulhERkROkcBcRKUAKdxGRAqRwFxEpQAp3EZECpHAXESlACncRkQKkcBcRKUD/HzuubkXmBsFgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 20         # Количество шагов\n",
    "x_begin = 0    # Начальная скорость\n",
    "x_end = 1      # Конечная скорость\n",
    "k = 0.01       # Небольшой костыль\n",
    "\n",
    "# Тут многовато уравнений в системе, но это на будущее\n",
    "eq = LSQR_Helper(n+n+1+n+3, 2*n+1)\n",
    "\n",
    "# x[i+1] = x[i] + u[i]\n",
    "for i in range(0, n):\n",
    "    eq.set_a(i,    1)\n",
    "    eq.set_a(i+1, -1)\n",
    "    eq.set_a(i+n+1,  1)\n",
    "    eq.set_b(0)\n",
    "    eq.begin_new_row()\n",
    "\n",
    "\n",
    "# x[0] = x_begin\n",
    "eq.set_a(0,1)\n",
    "eq.set_b(x_begin)\n",
    "eq.begin_new_row()\n",
    "\n",
    "# x[n] = x_end\n",
    "eq.set_a(n, 1)\n",
    "eq.set_b(x_end)\n",
    "eq.begin_new_row()\n",
    "\n",
    "solve_and_plot(eq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Получилось как-то плохо. Нужно добавить больше дополнительных ограничений. Попробуем добавить условие, чтобы в конце ускорение становилось равно 0.\n",
    "$$\n",
    "u_{n-1} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/scipy/sparse/compressed.py:746: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.\n",
      "  SparseEfficiencyWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8XGW97/HPL7fm2iZN0vSSpi3QAqVQC7EgXjbKxbZHQT1yc6uoHDh6YCNHt+fF2e6DbNzH1xbdHi+bo4LwQjwqVrxVbSmIiIqALW2gNy6llnZyadM0TS+Z3J/zx5pJJ8lMMm3mtla/79crr5lZ65mZX9dMfn3yPL/1LHPOISIiwZKX7QBERCT1lNxFRAJIyV1EJICU3EVEAkjJXUQkgJTcRUQCSMldRCSAlNxFRAJIyV1EJIAKsvXGNTU1bv78+dl6exERX3rhhRcOOOdqJ2qXteQ+f/58Nm7cmK23FxHxJTN7I5l2GpYREQkgJXcRkQBSchcRCaCsjbnH09/fTygUoqenJ9uhJFRcXEx9fT2FhYXZDkVEJKGcSu6hUIiKigrmz5+PmWU7nDGcc3R0dBAKhViwYEG2wxERSWjCYRkze9DM9pvZ1gT7zcy+aWY7zewlMzv/ZIPp6emhuro6JxM7gJlRXV2d039ZiIhAcmPuDwErxtm/ElgY+bkZ+PZkAsrVxB6V6/GJiEASwzLOuT+a2fxxmlwFPOy86/U9Z2aVZjbLOdeaohhFRHLGsd4BHvrLbnr7B0/6NS49u46lcytTGNVYqRhznwPsjXkcimwbk9zN7Ga83j0NDQ0peGsRkcz63Y59fGX9KwCc7B/yM6YW+yK5J805dx9wH0BjY6OuzC0ivhPqDAOw4+4VlBTlZzmaxFJR594MzI15XB/Z5jt33HEH99577/Dju+66i69+9atZjEhEck2oM0xNeVFOJ3ZITc99DXCrmT0CXAh0pWK8/V9+vY3tLYcnHVysxbOn8oX3npNw/7XXXsvtt9/OLbfcAsDq1atZv359SmMQEX8LdXYzp7Ik22FMaMLkbmY/Bi4BaswsBHwBKARwzn0HWAusAnYC3cDH0xVsui1btoz9+/fT0tJCe3s7VVVVzJ07d+Inisgpo7kzzNmzpmY7jAklUy1z/QT7HXBLyiKKGK+HnU5XX301jz76KG1tbVx77bVZiUFEctPQkCN0KMzli+uyHcqEcuoM1Vxw7bXXctNNN3HgwAGefvrpbIcjIjnkwNFe+gaGqK/K/WEZLRw2yjnnnMORI0eYM2cOs2bNynY4IpJD9kYqZeqrSrMcycTUc49jy5Yt2Q5BRHJQqLMbQD13EZEgida4z1FyFxEJjlBnmOqyIkqLcn/QQ8ldRCRJzYfCvhiSASV3EZGkhTq7fTEkA0ruIiJJcc7R3Bn2RaUMKLmLiCSl/WgvvT6pcQcldxGRpISGa9yV3EVEAiPkoxOYQMl9jN27d7NkyZLhx1/96le56667sheQiOSE6AlMflgREnL5DNV1d0Bbis8UnXkurPy31L6miJwSQp1hppcVUTYld9NmLPXcRUSSEOr0T4075HLPPUs97IKCAoaGhoYf9/T0ZCUOEcktoc5uzppZke0wkqae+yh1dXXs37+fjo4Oent7+c1vfpPtkEQky6I17n4Zb4dc7rlnSWFhIXfeeSfLly9nzpw5nHXWWdkOSUSy7MDRvkiNuz8qZUDJPa7bbruN2267LdthiEiO8NNSv1EalhERmYDfatxByV1EZEJ+Wsc9KueSu3e97dyV6/GJSOqFOrupKi2k3Cc17pBjyb24uJiOjo6cTaDOOTo6OiguLs52KCKSQSEfrQYZlVP/DdXX1xMKhWhvb892KAkVFxdTX1+f7TBEJINCnd0sqvNPjTvkWHIvLCxkwYIF2Q5DRGSYc47mQ2HeddaMbIdyQnJqWEZEJNd0HOujp3/IVycwgZK7iMi4/FgGCUruIiLjGj6Babp67iIigTFc465hGRGR4Ah1dlNZWkhFcWG2QzkhSu4iIuPw2zruUUkldzNbYWavmNlOM7sjzv4GM3vKzDab2Utmtir1oYqIZF6oM0x9pb8mUyGJ5G5m+cC9wEpgMXC9mS0e1eyfgdXOuWXAdcD/TXWgIiKZ5pwj1Nkd2J77cmCnc26Xc64PeAS4alQbB0yN3J8GtKQuRBGR7IjWuPsxuSdzhuocYG/M4xBw4ag2dwGPm9k/AGXAZSmJTkQki5qHV4MM4LBMkq4HHnLO1QOrgB+Y2ZjXNrObzWyjmW3M5fVjREQg9gQm//Xck0nuzcDcmMf1kW2xbgRWAzjnngWKgZrRL+Scu8851+ica6ytrT25iEVEMiR6ApOf1nGPSia5bwAWmtkCMyvCmzBdM6rNHuBSADM7Gy+5q2suIr4W6gwzraSQqT6rcYckkrtzbgC4FVgP7MCritlmZneb2ZWRZp8FbjKzF4EfAx9zuboou4hIkvxaKQNJLvnrnFsLrB217c6Y+9uBt6Y2NBGR7Ap1hjmttizbYZwUnaEqIhKHV+PuvyswRSm5i4jEcfBYH+H+Qd8Oyyi5i4jE4dd13KOU3EVE4vDrUr9RSu4iInE0H/JvjTsouYuIxBXqDDO1uIBpJf6rcQcldxGRuPxcKQNK7iIicfn5BCZQchcRGcPvNe6g5C4iMkZndz/dff6tcQcldxGRMaKrQSq5i4gEiN9PYAIldxGRMfy8jnuUkruIyCjNnWEqfFzjDkruIiJj+L1SBpTcRUTG8JK7f4dkQMldRGQEr8bd3ycwgZK7iMgIh7r7OdY3qGEZEZEgOV4GqZ67iEhgBOEEJlByFxEZIQgnMIGSu4jICKHObiqm+LvGHZTcRURGCHWGfX1mapSSu4hIjOZD/j+BCZTcRUSGHV/HXT13EZHA6Ar3c7R3QMldRCRIglIpA0ruIiLDglLjDkruIiLDoj33ueq5i4gER6gzTMWUAqaWFGQ7lElTcheRU8+mh+H+d4FzIzaHOruZU1WCmWUpsNRJKrmb2Qoze8XMdprZHQnaXGNm281sm5n9KLVhioik0Mu/heYXoOfQiM1BKYMEmPBvDzPLB+4FLgdCwAYzW+Oc2x7TZiHwP4G3Ouc6zWxGugIWEZm0ls3e7aE9UFIFeDXuzZ1hLjqtOouBpU4yPfflwE7n3C7nXB/wCHDVqDY3Afc65zoBnHP7UxumiEiKHG6Fo/u8+4f2HN8cHuBIQGrcIbnkPgfYG/M4FNkWaxGwyMyeMbPnzGxFvBcys5vNbKOZbWxvbz+5iEVEJiPaa4cRyX1vgMogIXUTqgXAQuAS4HrgfjOrHN3IOXefc67ROddYW1uborcWETkBrU1geVBYOiK5B+kEJkhizB1oBubGPK6PbIsVAp53zvUDfzOzV/GS/YaURCkikiotm6H2LMBGJfdTr+e+AVhoZgvMrAi4Dlgzqs0v8XrtmFkN3jDNrhTGKSIyec55yX3Wm6CyYUzPvTwA67hHTZjcnXMDwK3AemAHsNo5t83M7jazKyPN1gMdZrYdeAr4nHOuI11Bi4iclMMtcKwdZi87ntwjte7RMsgg1LhDcsMyOOfWAmtHbbsz5r4DPhP5ERHJTdHJ1NlvgsE+6D3s1bqXVBHq7A7MkAzoDFUROZW0NoHlQ90Sr+cOw0MzzZ1h5lQquYuI+E/LZphxNhSVxiT3vXSF+yM17sGolAEldxE5VcROpsKInnvQKmUgyTF3ERHf6wpBd4c33g7esgNF5V5ynxqsGndQcheRU8XwZOoy79ZsuGImVBFN7uq5i4j4S2sT5BVA3TnHt0WTe3k3ZUX5VJYGo8YdNOYuIqeK6GRqYUzvPJrcO8PUV5UGpsYdlNxF5FQwejI1qrIBervo7GgP1JAMKLmLyKng0B4Idx4fb4+KVsx07VFyFxHxndgzU2NFkvv0vlbmKLmLiPhMaxPkFXpnpsaqnAfAHDsQqDJIUHIXkVNBy2aoWwwFU0ZuL6liIL+UejugYRkREV9xDlqaxk6mAphxpHgW9daunruIiK907vZWfhw9mRrRXlBHQ94BqgJU4w5K7iISdIkmUyOaXS311h6oGndQcheRoGttgvwimLE47u7X+6up4BiED2U4sPRScheRYGvZ7C05MHoyNWJHuNK707U3g0Gln5K7iASXc9DyYvzJVKAr3M+rvdO9BzHXUw0CJXcRCa6Du6C3K+FkanNnmJCr8R4ouYuI+MQEk6mhzm46qWCwoAQOaVhGRMQfWpsgfwrUnh13d/OhMGC4aQ1w6I3MxpZmSu4iElwtTTBzCRQUxd0d6gxTUphP/vR5GpYREfGFoaHEZ6ZGhDq7qa8qwSLrugeJkruIBNPBXdB3JOFkKhC5SEeJtzpkzyHo6cpggOml5C4iwTTBZCowfAWm4XXdAzSpquQuIsHU2gQFxVB7Vtzdh3v66Qr3H++5Q6CGZnSBbBEJppbNMPNcyI+/IFhzZxgg0nOf5m0MUHJXz11EgmdoCFoTn5kK3pAM4F2BqbQaCkuV3EVEclrHTug7OsFkajeANyxjBtPmBqrWXcldRIInycnU4sI8qssiNfCVDYFaPCyp5G5mK8zsFTPbaWZ3jNPuP5uZM7PG1IUoInKCWpugoARqzkzYpDlSKTO8jnvAat0nTO5mlg/cC6wEFgPXm9mYhZHNrAL4NPB8qoMUETkhLZth1nmQn7hmJHSoe+R1UysbINwJPYczEGD6JdNzXw7sdM7tcs71AY8AV8Vp90Xgy0BPCuMTETkxQ4PQ+tK4k6kQcwJTVLQcMiBDM8kk9zlA7L82FNk2zMzOB+Y6536bwthERE7cgdeg/9i4k6lHevo51N0/8qLYlfO824AMzUx6QtXM8oCvAZ9Nou3NZrbRzDa2t7dP9q1FRMZKYjLVWw2S+D33Uyi5NwNzYx7XR7ZFVQBLgD+Y2W7gImBNvElV59x9zrlG51xjbW3tyUctIpJIa5NXs16zKGGT0MGYE5iiymq8SdhTKLlvABaa2QIzKwKuA9ZEdzrnupxzNc65+c65+cBzwJXOuY1piVhEZDwtm2HWUsjLT9gkWuM+pzKm524GlcGpdZ8wuTvnBoBbgfXADmC1c26bmd1tZlemO0ARkaQNDkDblgknU19uO0L5lAJqyket8x6gcsik1pZxzq0F1o7admeCtpdMPiwRkZNw4FXo7x53MnVgcIgntu/jnWfNOF7jHlXZAM2b0hxkZugMVREJjiQmU/+6+yAdx/pYuWTm2J2VDRA+CL1H0hRg5ii5i0hwtDZBUTlUn5GwyWNb2yguzOOSM+MUdQRoXXcldxEJjgkmU4eGHI9tbeOSRTMoLYozKh2gWncldxEJhiQmUzft6WT/kV5WnhtnSAYCVeuu5C4iwdD+Mgz0jDuZunZLG0X5ebzrrBnxG5TVeldvCkA5pJK7iATDBJOpzjke29rKOxbVUFEc/+pMXq17MMohldxFJBham6CoAqafHnf3i6EuWrp6WLFk1vivM22ukruISM5o2ez12vPip7V1W1spyDMuP7tu/NdRz11EJEcM9kPbVq9SJg7nHOu2tHHxGTVMK00wJBM1XOt+NA2BZo6Su4j43/4dMNibcDJ1e+th9hzsZlW8E5dGC8i67kruIuJ/w5Op8ZP7ui1t5BlcvniCIRkITK27kruI+F9rE0yZBlULxuxyzrF2aysXnVZNdfmUiV8rILXuSu4i4n8tm2H20riTqa/tP8qu9mPx15KJp3xGIGrdldxFxN8G+mDftoRnpq7b0oYZvPucJJO7WSDKIZXcRcTf9m+Hwb7E4+1bW2mcV8WMqcXJv2alkruISHaNc2bqrvajvNx2hJUTnbg0WgBq3ZXcRcTfWpugOP5k6rqtbQCsSHa8PaqyAbo7fF3rruQuIv7Wstkbkhl9VSW8tduXzq1kduy1UpMRLYf0ca27kruI+NdAL+zbHncyde/BbrY0dyV34tJoAbhoh5K7iPjXvm0w1B93MvWxyJDMCY+3Q0xy9285pJK7iPjXOJOp67a2cs7sqTRUl57465bNgPwpvp5UVXIXEf9qbYKSquNj5NHNXWE27TmU/IlLo+Xl+b4cUsldRPwrwWTq+uiQzLknMSQT5fMTmZTcRcSf+nu81SDjTKau3drGorpyTq8tP/nX93mtu5K7iPjTvm0wNDBmMrX9SC8bdh88uYnUWJUN0H0A+o5N7nWyRMldRPypZZN3O2oydf22NpyDleee5Hh71PDSv/4sh1RyFxF/am2C0mpvbDzGY1vbOK2mjDPrKib3+j5f+lfJXUT8p6cLXv4tzLt4xGRq57E+nt3VwYolM7E4Z6yekOErMim5i4hkxl/+A8Kd8PZ/HLH5ie37GBxyrJpMlUxUeR3kF6nnLiKSEUfb4dl7YfH7xoy3r93aSn1VCefMnjr598nL83U5pJK7iPjLn/4dBnrgXf88YnNXuJ9ndh5g1bmzJj8kE+XjcsikkruZrTCzV8xsp5ndEWf/Z8xsu5m9ZGZPmtm8eK8jIjIph/bAxgdg2d9DzcIRu57csY/+QXfiy/uOx8dnqU6Y3M0sH7gXWAksBq43s8Wjmm0GGp1z5wGPAvekOlAREf7wZcDg78b0MVm3tY1Z04p5U31l6t6vsgGOtUNfd+peM0OS6bkvB3Y653Y55/qAR4CrYhs4555yzkX/9c8B9akNU0ROeftfhhd/BMtvgmlzRuw62jvA06+28+5zZpKXl6IhGfD1uu7JJPc5QOy/LBTZlsiNwLp4O8zsZjPbaGYb29vbk49SROSpf4XCMnjbZ8buenk/fQNDqamSieXjWveUTqia2YeBRuAr8fY75+5zzjU65xpra2tT+dYiEmTNL8COX8PFt0JZ9Zjd67a2UlM+hQvmVaX2fX28rnsyyb0ZiD0FrD6ybQQzuwz4PHClc643NeGJiABP3u2djfqWW8bsCvcN8tTL7axYUkd+KodkAMpnQl6hL5cgSCa5bwAWmtkCMysCrgPWxDYws2XAd/ES+/7Uhykip6xdT8OuP8DbPwtTxi4p8PSr7YT7Bye/UFg8Pl7XfcLk7pwbAG4F1gM7gNXOuW1mdreZXRlp9hWgHPipmTWZ2ZoELycikjzn4Ml/gan10Hhj3CbrtrZSVVrIhQumpycGn9a6FyTTyDm3Flg7atudMfcvS3FcIiLe+jHNL8CV34LC4jG7ewcGeXLHfv7TubMoyE/TOZnT5sKr69Pz2mmkM1RFJDcNDcLvvwjVC2Hph+I2+fNrBzjaOzD55X3HUzkPju2H/nD63iMNlNxFJDe9tBraX4Z3fR7y4w8yrNvaRkVxARefXpO+OIYrZvw1qarkLiK5Z6AP/vAlmLUUzr4qbpP+wSGe2L6PyxfXUVSQxlTm01r3pMbcRUQy6oWHvGT6nv/jVazE8ezrHXSF+9NTJRPLp7Xu6rmLSG7pPQp//ArMexucfmnCZr9qaqGsKJ+3L0zjkAxARbTW3V89dyV3Ecktz3/bm8C87AsjrrIU61dNzfxsU4irG+dSXJif3njy8mFavZK7iMhJ6z4Iz3wLFq2EucvjNtmw+yCf++lLXLhgOv+06uzMxFXZ4LvFw5TcRSR3PPN16D0Ml/6vuLvf6DjGzQ9vZE5VCd/9yAXpnUiN5cMTmZTcRSQ3HG6F578L514NdeeM2d3V3c/HH9qAAx782JupLC3KXGyVDXB0n69q3ZXcRSQ3/PEeGBqAd/7TmF19A0N88v+9wN6D3dz3kUYW1JRlNrZoxUxXKLPvOwlK7iKSfQd3waaH4YKPwfQFI3Y55/j8L7bw7K4O7vngeSxP1xoy4/FhOaSSu4hk31Nf8soN3/G5Mbu+/fTr/PSFELddupD3L8vSRd58eCKTkruIZFfbVtjyKFz0Sa+mPMZvX2rlnsde4cqls/nvly1M8AIZUDEL8gqU3EVEkvb7L0LxVHjrp0ds3rSnk8+sbqJxXhX3fPA8LEHNe0b4sNZdyV1EsmfPc/DqY15iLzl+iby9B7u5+eGN1E0t5rsfuSD9Jyolw2flkEruIpIdvUfg8X+Gshlw4SeHNx/u6ecTD22gb2CIBz/2ZqrLp2QxyBiVDb5aGVILh4lI5jW/AI/e6FWfvP+7UOSVNvYPDnHLDzfxtwPHePgTyzljRnmWA40xrQGOtkF/T9wLh+Qa9dxFJHOGhuDPX4cHrvBq2j+2Fs67BvBKHr+wZht/eu0AX3r/uVx8RpoXBDtRPqt1V89dRDLjSBv84r96F7tefBW89xsjxtkf+PPf+NHze/jUJadzzZvnZi/ORGJr3WvOyG4sSVByF5H0e3U9/PJT0NcN7/0mnP/RESs+rt/Wxv9eu4NV587kc1ecmcVAx+GzWncldxFJn/4e+N0X4PnvQN258MEHoHZk8t4S6uL2R5o4r76Sr13zJvLysljyOB6f1boruYtIerS/4k2a7tsCF34KLrtrzERky6EwN35/A9PLivjeRxtzo+QxkfwCmDpHyV1ETlHOwabvw7o7oKgUPrQaFr17RJP+wSHWbmnl6797jXDfID/41IXUVuRIyeN4fFTrruQuIqkT7oRffxq2/wpOu8Qrc4xZUqAr3M8jf93DQ3/ZTWtXD6fVlnHfRxs5c2ZF1kI+IZXzYOfvsh1FUpTcRSQ13ngWfvZfvFrwy++Gt/zD8MWt9x7s5sFn/sbqDXs51jfIW06r5l/ft4R3njkjd8fY46n0T627kruITM7ggHdB6z/e4/Vsb3wc5lwAeOvDPPCnv7Fuayt5Zrx36WxufNsClsyZluWgT1JlpETzcDNUn57dWCag5C4iJ+dwC7z+e3jh+xD6Kyy9HlZ9hcHCch7f0sr9f9rFpj2HmFpcwM3vOJ0bLp7HrGkl2Y56cmJr3ZXcRSQQ+sPwxl+8hL7zSWjf4W0vnwkfuJ9jZ36A1Rv38uAzG9h7MMzc6SXc9d7FXN04l7IpAUk1Pqp1D8gRF5GUcw7aXz6ezN94BgZ6IL8IGt4Cb7oeTr+U1uLTeOjZN/jRz5/kSM8AF8yr4vOrzubyxTPJ99N4ejIqZoPlK7mLiM90H/SWB3j9SXj9KW9sGaBmEQPLbiBUfTEv5S9m+4EhXn3tCK/8aT/Nh94gz2Dlklnc+PYFnN9QNe5b+Fp+AUzzR617UsndzFYA3wDyge855/5t1P4pwMPABUAHcK1zbndqQxWRlOkPe2PmR1rhcCsceMVL5i2bwA0xWDSNfbUXsa32Bn7ffy7PHyxl95+PMeQAXqMw3zi9tpzz51XxoQsbuHLpbOZOL832vyozKucFI7mbWT5wL3A5EAI2mNka59z2mGY3Ap3OuTPM7Drgy8C16QhYRMYxNATdHXCkxUvaMbfucCtDXV5Cz+89NPJp5LG7+GyeKbqWNcfOYlPPaQwezifPYH51KYvqKnjP0tmcWVfBorpy5teUUZh/ii4qW9ngDVXluGR67suBnc65XQBm9ghwFRCb3K8C7orcfxT4DzMz55xLYawi/uAcuCEYGvSWtR3+GcQN9eMGBxgcHMAN9uOGBhkciNwfHGCgv4fB3mMM9RxjsO8YQ31hXN8xXF839Hfj+rqx/jA2EMb6u8kbCJM3GCZvIExRXxelfQcocP0jwhnCOEAlrUNV7HNVtLk30+aq2Oem00YVbW46ra6a6SVVnDmzgsaZFXyorpxFdRWcXlue20sCZENlg/cXz0AvFOTuWbXJJPc5QOzlR0LAhYnaOOcGzKwLqAYOpCLIWBt+/g1qt96f6pfNef6clprc/+0W83yboJ9go95rxHMj949vc8PHM+425zAS/+ThIHI7cvsQAIUMjhOn93Myfd5+l0+YKYQpIuwit0yJ3C/hMDUcsDdzuKCWY1NqCZfU0V86E8pnUF5SzLSSQqaWFDKtpJClkdvotumlRZQUKYknZVqk1r0rlNPlkBmdUDWzm4GbARoaGk7qNQrKqzlYuiCVYfmGm0SKz9Z/DpOJGRixLOzY1xr5OO57RZ4fm86H29vobTbczttnOMuL3I58jEVSeszj6O2QFeLy8sHycHkFOCvAWT7kFUS2F0BePi6vwFtlMC//+G1BMRSWYEVl5BWVej9TyiiYUkZBYRGFBXkU5edRVJBHSX4eU6OP8/OYUpinXnYm1C2Gxe/LdhQTSia5NwOxK+fXR7bFaxMyswJgGt7E6gjOufuA+wAaGxtPqlu37IoPwxUfPpmniohM3uxlcM33sx3FhJL563ADsNDMFphZEXAdsGZUmzXADZH7HwR+r/F2EZHsmbDnHhlDvxVYj1cK+aBzbpuZ3Q1sdM6tAR4AfmBmO4GDeP8BiIhIliQ15u6cWwusHbXtzpj7PcDVqQ1NRERO1ilaqCoiEmxK7iIiAaTkLiISQEruIiIBpOQuIhJAlq1ydDNrB944yafXkIalDVJAcZ0YxXXicjU2xXViJhPXPOdc7USNspbcJ8PMNjrnGrMdx2iK68QorhOXq7EprhOTibg0LCMiEkBK7iIiAeTX5H5ftgNIQHGdGMV14nI1NsV1YtIely/H3EVEZHx+7bmLiMg4cjq5m9kKM3vFzHaa2R1x9k8xs59E9j9vZvMzENNcM3vKzLab2TYz+3ScNpeYWZeZNUV+7oz3WmmIbbeZbYm858Y4+83Mvhk5Xi+Z2fkZiOnMmOPQZGaHzez2UW0ydrzM7EEz229mW2O2TTezJ8zstchtVYLn3hBp85qZ3RCvTQpj+oqZvRz5nH5hZpUJnjvuZ56m2O4ys+aYz2tVgueO+/ubhrh+EhPTbjNrSvDctByzRLkha98v51xO/uAtL/w6cBpQBLwILB7V5r8B34ncvw74SQbimgWcH7lfAbwaJ65LgN9k4ZjtBmrG2b8KWId36aGLgOez8Jm24dXpZuV4Ae8Azge2xmy7B7gjcv8O4Mtxnjcd2BW5rYrcr0pjTFcABZH7X44XUzKfeZpiuwv4xyQ+63F/f1Md16j9/w7cmcljlig3ZOv7lcs99+ELczvn+oDohbljXQVEL4nyKHCpmaX1inLOuVbn3KbI/SPADrxryPrBVcDDzvMcUGlmszL4/pcCrzvnTvbktUlzzv0R75oDsWK/R98H4l1D7d3AE865g865TuAJYEW6YnLOPe6cG4g8fA4RXBWqAAAC+UlEQVTvCmgZl+B4JSOZ39+0xBXJAdcAP07V+yUZU6LckJXvVy4n93gX5h6dREdcmBuIXpg7IyLDQMuA5+PsfouZvWhm68zsnAyF5IDHzewF865XO1oyxzSdriPxL1w2jldUnXOuNXK/DaiL0yabx+4TeH9xxTPRZ54ut0aGjB5MMMyQzeP1dmCfc+61BPvTfsxG5YasfL9yObnnNDMrB34G3O6cOzxq9ya8oYelwLeAX2YorLc5584HVgK3mNk7MvS+EzLvEo1XAj+Nsztbx2sM5/2NnDMlZGb2eWAA+GGCJtn4zL8NnA68CWjFGwLJJdczfq89rcdsvNyQye9XLif3E7kwNzbOhblTzcwK8T68Hzrnfj56v3PusHPuaOT+WqDQzGrSHZdzrjlyux/4Bd6fxrGSOabpshLY5JzbN3pHto5XjH3R4anI7f44bTJ+7MzsY8B7gL+PJIUxkvjMU845t885N+icGwLuT/CeWfmuRfLAB4CfJGqTzmOWIDdk5fuVy8k9Jy/MHRnPewDY4Zz7WoI2M6Nj/2a2HO84p/U/HTMrM7OK6H28Cbmto5qtAT5qnouArpg/F9MtYW8qG8drlNjv0Q3Ar+K0WQ9cYWZVkWGIKyLb0sLMVgD/A7jSOdedoE0yn3k6Youdp3l/gvdM5vc3HS4DXnbOheLtTOcxGyc3ZOf7leoZ41T+4FV3vIo36/75yLa78b7wAMV4f+bvBP4KnJaBmN6G92fVS0BT5GcV8Engk5E2twLb8CoEngMuzkBcp0Xe78XIe0ePV2xcBtwbOZ5bgMYMfY5leMl6Wsy2rBwvvP9gWoF+vHHNG/HmaZ4EXgN+B0yPtG0Evhfz3E9Evms7gY+nOaadeGOw0e9YtCpsNrB2vM88A8frB5Hvz0t4iWvW6Ngij8f8/qYzrsj2h6Lfq5i2GTlm4+SGrHy/dIaqiEgA5fKwjIiInCQldxGRAFJyFxEJICV3EZEAUnIXEQkgJXcRkQBSchcRCSAldxGRAPr/KZNRd+cyBigAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#u[n-1] = 0\n",
    "eq.set_a(2*n, 1)\n",
    "eq.set_b(0)\n",
    "\n",
    "solve_and_plot(eq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Все равно неудовлетворительно. Больше ограничений богу ограничений! Это ограничение я подсмотрел в исходной статье. Пусть скорость будет стремиться к целевой как можно быстрее. Для этого мы будем требовать, чтобы все $x_i = x_{end}$: \n",
    "$$\n",
    "\\begin{cases}\n",
    "...\\\\\n",
    "x_0 = x_{end}\\\\\n",
    "x_1 = x_{end}\\\\\n",
    "... \\\\\n",
    "x_n = x_{end} \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Но чтобы это ограничение не мешало задаче (ведь в начале скорость равна $x_{start}$), каждое уравнение для этого ограничения умножим на некоторый небольшой коэффициент, например 0.01, чтобы остальные уравнения имели больший вес. Как-то так."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/scipy/sparse/compressed.py:746: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.\n",
      "  SparseEfficiencyWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHFlJREFUeJzt3X9wHPd53/H3g98ACZIASNMSQYp0zNShpNTSYGjHjly1dh1K04pNG1lSk6kTe8TJVKqsJk2HHXcYVv2jdaIkM2nZpGzjOvGklhS1SVmXKqMkrtWmliraliVTsixaokxQP0jhjj90B/AOwNM/dg88gnfAHniHw+7385rB4O527+7h4vDhF9/dfdbcHRERyZaOdhcgIiLNp3AXEckghbuISAYp3EVEMkjhLiKSQQp3EZEMUriLiGSQwl1EJIMU7iIiGdTVrjdev369b926tV1vLyKSSt/85jffcfcNi63XtnDfunUrR48ebdfbi4ikkpm9nmQ9TcuIiGSQwl1EJIMU7iIiGaRwFxHJIIW7iEgGLRruZvZFMzttZt+ts9zM7LfN7LiZPW9mNze/TBERaUSSkfuXgF0LLL8N2B5/7QF+5+rLEhGRq7Hoce7u/pSZbV1gld3AH3h0vb6nzWydmV3j7m82qcbLvf4NeOUI/I190NHYrFJpepb/9JevUbg43ZLSRCQ9ers7GR3qZ3RogM3D/WxY3YuZNf19Zmadt89PcTJX5GR+kvF8kY9/YCM3jq5t+ntVa8ZJTJuAk1X3x+PHrgh3M9tDNLpny5YtS3u3N74F/+e34KOfg/6hhp569ESOf/XE9+Jalvb2kpy7trOsXPMvH93X3REF/VA/m4cH2ByHfvTYAGsHuuu8jpMrlDiZn4wDvMjJXBTiJ3NFTp2dpDxz6c3MYP3q3lSEe2LufhA4CDA2Nra0K3MPjETfi7mGw/2dQgmAJ//xx9i+cXBJby8i2TBZmokCOA7j6mA++nqeC1OX/4U/2Nc1F/jrV/fGo/FJTuaLFEszl607vKqHzUP9XL9pLT91w3vj50X/cWwa6qe3q7Pl/75mhPspYHPV/dH4sdaYC/cJGPmRhp6aj8N93UBPs6sSkZTp7+lk+8bBugO9c5NlTuaK8Qh8Mg7+Ij84U+CZ13K8d00fm4cH+Mj7Ry6FdzzSX93bts4uc5pRwSHgfjN7BPgQcK5l8+0AA8PR92Ku4afmi5Vwr/3nlYhIxdr+btZuWssNm1o7fdIqi4a7mX0FuBVYb2bjwK8C3QDu/rvAYeB24DhQBH6hVcUCl4/cG5QvlFjT10V3pw7vF5FsS3K0zD2LLHfgvqZVtJirCPdcsczwKk3JiEj2pW8I27MaOnuWPHIfUriLSADSF+5m0eh9KSP3Qolh7UwVkQCkL9wB+oeXvENVI3cRCUE6w31guOGRe+VEA825i0gIUhrujU/LTJZnuDg9y5CmZUQkAMGEey4+gWl4lY5xF5HsS2+4T+ZhdmbxdWP5QhlAI3cRCUJ6wx2HqXOJn1I5O1U7VEUkBCkOdxqampkLd43cRSQAKQ33Sn+Z5OF+ac5d4S4i2ZfScF/CyL1QwixqBiQiknXBhHuuWGJdfzedHbp6hIhkX0rDvfFpmXyhrJ2pIhKMdIZ79wB09TU8566+MiISinSG+1zzsOT9ZdRXRkRCks5wh4b7y2jkLiIhSXG4J29B4O4auYtIUIII93cvTlOecYZ07VQRCUTKwz3ZnPvZYtxXRiN3EQlEusN96izMTC+66tzZqZpzF5FApDvcIeoOuYicmoaJSGBSHO7JT2TKq6+MiAQmxeGevAWBpmVEJDTpDff+BkbuxRKdHcZgX1eLixIRWRnSG+4NjdzLDA1006GmYSISiBSHe2Nz7rpIh4iEJL3h3t0P3asSHeue09mpIhKY9IY7JD5LNa++MiISmJSH+zBMLj5yzxfLDK1S6wERCUfKw33xkftc0zCN3EUkIInC3cx2mdnLZnbczPbWWL7FzL5mZt82s+fN7Pbml1pDgnA/PzXNzKzrBCYRCcqi4W5mncAB4DZgB3CPme2Yt9o/Bx5z95uAu4F/1+xCa0rQPKxydqpG7iISkiQj953AcXd/1d1LwCPA7nnrOLAmvr0WeKN5JS5gYAQunofpUt1VKn1lNHIXkZAkCfdNwMmq++PxY9X2Az9nZuPAYeAf1XohM9tjZkfN7OiZM2eWUO48A0PR9wV2qs6N3BXuIhKQZu1QvQf4kruPArcDXzazK17b3Q+6+5i7j23YsOHq3zXBWarqKyMiIUoS7qeAzVX3R+PHqn0WeAzA3b8B9AHrm1HgghKEe36u3a8OhRSRcCQJ92eB7Wa2zcx6iHaYHpq3zg+BjwOY2Y8RhXsT5l0WkWjkXqa701jdq6ZhIhKORcPd3aeB+4EjwEtER8UcM7OHzOyOeLVfBu41s+8AXwF+3t29VUXPSTJyj/vKmKlpmIiEI9Fw1t0PE+0orX5sX9XtF4GPNre0BOba/i6wQ7VY0pEyIhKcdJ+h2tUDvWsWDfd1A5pvF5GwpDvcIeovs8jRMhq5i0hoMhDuC7cgyBfLOjtVRIKT6XCfmXXOas5dRAKU/nDvH647535+ssysq6+MiIQn/eG+wMhdfWVEJFQZCPdhKBegPHnFIvWVEZFQZSDcKycyXTk1o74yIhKqDIX7lVMz6isjIqHKeLiXAc25i0h4sh3uhRI9XR30d3cuc1EiIu2VnXCfzF+xKFcoMaymYSISoPSHe398NaY6c+46UkZEQpT+cO/sgr51NcM96iujnakiEp70hzvUbR6mvjIiEqqMhHvts1TVEVJEQpXZcJ+emeXcpEbuIhKmDIX75Weonp3UMe4iEq6MhHs851512Vb1lRGRkGUk3EdgegrKxbmH1FdGREKWnXCHy+bdK60H1FdGREKUsXC/NO8+1zRMI3cRCVDGwv3SyL0yLaNwF5EQZSzcq0buhRL93Z3096hpmIiEJxvh3j8cfa8euevC2CISsIyE+zrALt+hWihpZ6qIBCsb4d7RGXWHvGzkrrNTRSRc2Qh3uKIFQV59ZUQkYJkOd43cRSRUGQv36GiZ0vQsFy5Oa+QuIsFKFO5mtsvMXjaz42a2t846nzKzF83smJn95+aWmUBVT/ezk+orIyJh61psBTPrBA4AfxMYB541s0Pu/mLVOtuBfwZ81N3zZvaeVhVcV2Vaxp18IW49MKCjZUQkTElG7juB4+7+qruXgEeA3fPWuRc44O55AHc/3dwyExgYgdkylN5V0zARCV6ScN8EnKy6Px4/Vu1HgR81s780s6fNbFezCkysqgXBXF8ZTcuISKAWnZZp4HW2A7cCo8BTZnaju5+tXsnM9gB7ALZs2dKkt45VhXuusB7QhTpEJFxJRu6ngM1V90fjx6qNA4fcvezurwHfJwr7y7j7QXcfc/exDRs2LLXm2gYqLQhycxfqWKc5dxEJVJJwfxbYbmbbzKwHuBs4NG+dPyEatWNm64mmaV5tYp2Lqx65F0us7u2it0tNw0QkTIuGu7tPA/cDR4CXgMfc/ZiZPWRmd8SrHQEmzOxF4GvAr7j7RO1XbJGBS83D1FdGREKXaM7d3Q8Dh+c9tq/qtgO/FH+1R+9asM545F7WkTIiUle5XGZ8fJypqal2l1JXX18fo6OjdHcvbaDarB2q7dfRMXciU75QYmS1wl1EahsfH2dwcJCtW7diZu0u5wruzsTEBOPj42zbtm1Jr5Gd9gMwdyJTvljSyF1E6pqammJkZGRFBjuAmTEyMnJVf1lkMNxz8Zy7wl1E6lupwV5xtfVlLNyHmS1OUCjNqPWAiAQtY+E+gheig3Q0cheRkGUu3Dsmc4Brzl1EVqy9e/dy4MCBufv79+/n4Ycfbup7ZOdoGYD+YcxnWENRI3cRSeRf/PdjvPjG+aa+5o5r1/Crf/v6usvvuusuHnzwQe677z4AHnvsMY4cOdLUGrIV7vFZqkN2QX1lRGTFuummmzh9+jRvvPEGZ86cYWhoiM2bNy/+xAZkMtyHuaBL7IlIIguNsFvpzjvv5PHHH+ett97irrvuavrrZzLch+yCmoaJyIp21113ce+99/LOO+/w9a9/vemvn7EdqlF/mWt7inR3ZuufJiLZcv3113PhwgU2bdrENddc0/TXz+TIfVNPsc2FiIgs7oUXXmjZa2dreNs7yDRdbOwqtLsSEZG2yla4m3HO1rC+4912VyIi0lbZCncgzyBDpnAXkbBlLtzfmV3NWm/uCQkiImmTqXCfLM0wMbuKwdlz7S5FRKStMhXuuWKJvA/SP61wF5GwZSrc84USOQbpLZ2D2Zl2lyMi0jaZCvdcIRq5G7MwpdG7iIQrU+GeL5bI+WB0pzjR3mJEROo4ceIEN9xww9z9hx9+mP379zf1PTJ1hmq+UCJPdbhvb2s9IpICT+yFt5p8puh7b4Tb/nVzX7NBmRq554rleeEuIhKmzI3cZ3qHwVG4i0gybRhhd3V1MTs7O3d/amqq6e+RsZF7aa55mMJdRFaqjRs3cvr0aSYmJrh48SJf/epXm/4emRu5968ahIt9UMy1uxwRkZq6u7vZt28fO3fuZNOmTXzgAx9o+ntkKtxzhRKjQwMwOaxwF5EV7YEHHuCBBx5o2etnalomXywxvKo7mprRtIyIBCwz4e7u5Atlhlb1RFdkUriLSMAyE+6F0gylmVmGB3o0cheR4CUKdzPbZWYvm9lxM9u7wHp/z8zczMaaV2Iy+UIJIB65K9xFZGHu3u4SFnS19S0a7mbWCRwAbgN2APeY2Y4a6w0CnwOeuaqKligXh/vcyH3qLMxMt6MUEVnh+vr6mJiYWLEB7+5MTEzQ19e35NdIcrTMTuC4u78KYGaPALuBF+et9y+BLwC/suRqrkK+OG/kDjCZh9Ub2lGOiKxgo6OjjI+Pc+bMmXaXUldfXx+jo6NLfn6ScN8EnKy6Pw58qHoFM7sZ2Ozu/8PM2hruw5UdqhBNzSjcRWSe7u5utm3b1u4yWuqqd6iaWQfwm8AvJ1h3j5kdNbOjzf4fM1coA1XTMqB5dxEJVpJwPwVsrro/Gj9WMQjcAPwvMzsBfBg4VGunqrsfdPcxdx/bsKG5I+p8oUSHwWBfV9W0jE5kEpEwJQn3Z4HtZrbNzHqAu4FDlYXufs7d17v7VnffCjwN3OHuR1tScR25YomhgR46OuzyaRkRkQAtGu7uPg3cDxwBXgIec/djZvaQmd3R6gKTyhdK0c5UgH6Fu4iELVFvGXc/DBye99i+OuveevVlNS5XKEXz7QA9A9A9oP4yIhKszJyhmi+WGFrVfekBncgkIgHLTLjnCuXoMMgK9ZcRkYBlItzdnbPxDtU5GrmLSMAyEe4XLk4zPevzRu4KdxEJVybCfa5p2BUjd+1QFZEwZSLc55qGzR+5XzwP06U2VSUi0j6ZCPdKX5l1A9VHy8THuussVREJUCbCfa6vTPXIfe5EJoW7iIQnE+F+2YU6KtQ8TEQClolwzxVLdHUYg71VJ9wq3EUkYJkI90pfGTO79KDCXUQClolwv6yvTMWA5txFJFyZCPcr+soAdPVCz6BG7iISpIyE+7y+MhXqLyMigcpGuBfm9ZWpUAsCEQlU6sN9dtbJF0t1Ru4KdxEJU+rD/fxUmVlngZG7dqiKSHhSH+65uROYuq9cODCi9gMiEqTUh3ulr0ztkfsQlN6F8tQyVyUi0l6pD/eafWUqKicyafQuIoFJfbjX7OVeobNURSRQqQ/3XLFGL/cKhbuIBCr14Z4vlOjp6mCgp/PKhQp3EQlU+sO9GPWVuaxpWMVcuGvOXUTCkvpwzxXKl/dxr9Y/FH3XyF1EApP6cI/OTq1xjDtAZzf0rVW4i0hw0h/u9frKVKgFgYgEKPXhnqvXV6ZC4S4iAUp1uE/PzHJussy6hUbu/cPaoSoiwUl1uJ+bLOMOwwN15txBzcNEJEiJwt3MdpnZy2Z23Mz21lj+S2b2opk9b2Z/bmbXNb/UK831lVlwWkYX7BCR8Cwa7mbWCRwAbgN2APeY2Y55q30bGHP3HwceB36t2YXWsmBfmYqBEZiehFJxOUoSEVkRkozcdwLH3f1Vdy8BjwC7q1dw96+5eyU9nwZGm1tmbbmF+spU6CxVEQlQknDfBJysuj8eP1bPZ4EnrqaopPIL9ZWpULiLSIC6mvliZvZzwBjw1+os3wPsAdiyZctVv9+CvdwrFO4iEqAkI/dTwOaq+6PxY5cxs08AnwfucPeLtV7I3Q+6+5i7j23YsGEp9V4mXyjR391Jf62mYRXqLyMiAUoS7s8C281sm5n1AHcDh6pXMLObgH9PFOynm19mbblCeeEpGdDIXUSCtGi4u/s0cD9wBHgJeMzdj5nZQ2Z2R7zarwOrgT8ys+fM7FCdl2uqfLFU+9qp1frXAaZwF5GgJJpzd/fDwOF5j+2ruv2JJteVSG6xvjIAHZ1RwCvcRSQgqT5DNV9MEO4QTc3oOqoiEpBUh3uusEjTsAo1DxORwKQ23Mszs1yYmk4+ctfRMiISkNSG+6UTmBbZoQrqLyMiwUlvuMd9ZRZsGlZRmZZxb3FVIiIrQ3rDvTJyTzotM1OC0rstrkpEZGVIb7gXErT7rdCJTCISmNSGey5J07AKhbuIBCa14V4Zua9b6CpMFeovIyKBSW245wplVvd20du1QNOwiv7h6LtG7iISiNSGe75YSjZqh+hQSNDIXUSCkdpwT3x2KkDfOrAOjdxFJBipDffEfWUAOjqiqRmFu4gEIrXh3tDIHdRfRkSCktpwzydp91tN/WVEJCCpDPeL0zMUSjPJ+spUqL+MiAQkleF+tthAX5kKTcuISEBSGe65QgN9ZSrUPExEApLKcG+or0zFwAj4DEyda1FVIiIrRyrDvaG+MhUDOktVRMKRynBvqK9MhfrLiEhAUhnuucqFOhqdcwddKFtEgpDKcM8XSwz2ddHd2UD5mpYRkYCkMtwbPjsV1NNdRIKSynBvqK9MRe8a6OhSuItIEFIb7g2P3M10IpOIBCOd4V4oNz5yB/WXEZFgpDLcozn3Bg6DrNDIXUQCkbpwnyzNMFmeaezs1Ao1DxORQKQu3PPFJfSVqdDIXUQCkSjczWyXmb1sZsfNbG+N5b1m9mi8/Bkz29rsQityS+krU9E/DJN5mJ1tclUiIivLouFuZp3AAeA2YAdwj5ntmLfaZ4G8u78f+C3gC80utKIycl/yDlWfhamzTa5KRDKvVIDjfwZP7oODfx3ePd3uihbUlWCdncBxd38VwMweAXYDL1atsxvYH99+HPi3Zmbuze+vO9fud6k7VCE6YqZyxqqISC3lKRh/Fk78b3jtKRg/CrNl6OiG0TEovAOr39PuKutKEu6bgJNV98eBD9Vbx92nzewcMAK804wiq821+13qyB3g0Z+F7oEmVpVl7e5/b21+f1nRegdh3RYYug7WbY2/XxeFrjX42ZmZhje+Da99PQrzk8/A9BRYB1zzQfiJ+2Dbx2DLh6FnVUv+Oc2UJNybxsz2AHsAtmzZsqTXeM+aPm7Zvp61/UsYuY+OwfU/DRffXdJ7B6vRX5Jm0YVVZEEOk2fh+/8TCmcuX9TVXxX6111+e+g66B+K9r29/UIU5K89Ba//XyjF2bDxBhj7DGy9Ba77CPSvW/5/3lVKEu6ngM1V90fjx2qtM25mXcBa4IrDUtz9IHAQYGxsbEm/ubffeA2333jNUp4a/YDu/NLSnisiK1epAGd/CPnXo+9nX4f8iej7D5+Bi/Mu0tO7NvqjsHLxnpH3w49/KhqZb70FVq1f7n9B0yUJ92eB7Wa2jSjE7wb+/rx1DgGfBr4B/AzwF62YbxcRqalnFbznx6KvWibzcfC/fuk/gJkSXPdR2HYLrLl2eetdBouGezyHfj9wBOgEvujux8zsIeCoux8Cfg/4spkdB3JE/wGIiKwM/UPR17UfbHclyybRnLu7HwYOz3tsX9XtKeDO5pYmIiJLlbozVEVEZHEKdxGRDFK4i4hkkMJdRCSDFO4iIhmkcBcRySCFu4hIBlm7TiQ1szPA60t8+npa0JSsCVRXY1RX41ZqbaqrMVdT13XuvmGxldoW7lfDzI66+1i765hPdTVGdTVupdamuhqzHHVpWkZEJIMU7iIiGZTWcD/Y7gLqUF2NUV2NW6m1qa7GtLyuVM65i4jIwtI6chcRkQWs6HA3s11m9rKZHTezvTWW95rZo/HyZ8xs6zLUtNnMvmZmL5rZMTP7XI11bjWzc2b2XPy1r9ZrtaC2E2b2QvyeR2ssNzP77Xh7PW9mNy9DTX+lajs8Z2bnzezBeess2/Yysy+a2Wkz+27VY8Nm9qSZvRJ/H6rz3E/H67xiZp9ucU2/bmbfi39Of2xmNa/zttjPvEW17TezU1U/r9vrPHfB398W1PVoVU0nzOy5Os9tyTarlw1t+3y5+4r8IrowyA+A9wE9wHeAHfPW+YfA78a37wYeXYa6rgFujm8PAt+vUdetwFfbsM1OAOsXWH478ATRBcY+DDzThp/pW0TH6bZlewEfA24Gvlv12K8Be+Pbe4Ev1HjeMPBq/H0ovj3Uwpo+CXTFt79Qq6YkP/MW1bYf+CcJftYL/v42u655y38D2Lec26xeNrTr87WSR+47gePu/qq7l4BHgN3z1tkN/H58+3Hg42atvZqzu7/p7t+Kb18AXgI2tfI9m2g38AceeRpYZ2ZLvCDtknwc+IG7L/Xktavm7k8RXS2sWvXn6PeBv1PjqT8FPOnuOXfPA08Cu1pVk7v/qbtPx3efJrp28bKrs72SSPL725K64gz4FPCVZr1fwprqZUNbPl8rOdw3ASer7o9zZYjOrRP/IpwDRpalOiCeBroJeKbG4p8ws++Y2RNmdv0yleTAn5rZN81sT43lSbZpK91N/V+4dmyvio3u/mZ8+y1gY4112rntPkP0F1cti/3MW+X+eMroi3WmGdq5vW4B3nb3V+osb/k2m5cNbfl8reRwX9HMbDXwX4AH3f38vMXfIpp6+KvAvwH+ZJnK+kl3vxm4DbjPzD62TO+7KDPrAe4A/qjG4nZtryt49DfyijmEzMw+D0wDf1hnlXb8zH8H+BHgg8CbRFMgK8k9LDxqb+k2WygblvPztZLD/RSwuer+aPxYzXXMrAtYC0y0ujAz6yb64f2hu//X+cvd/by7vxvfPgx0m9n6Vtfl7qfi76eBPyb607hakm3aKrcB33L3t+cvaNf2qvJ2ZXoq/n66xjrLvu3M7OeBvwX8bBwKV0jwM286d3/b3WfcfRb4D3Xesy2ftTgH/i7waL11WrnN6mRDWz5fKzncnwW2m9m2eNR3N3Bo3jqHgMpe5Z8B/qLeL0GzxPN5vwe85O6/WWed91bm/s1sJ9F2bul/Oma2yswGK7eJdsh9d95qh4B/YJEPA+eq/lxstbqjqXZsr3mqP0efBv5bjXWOAJ80s6F4GuKT8WMtYWa7gH8K3OHuxTrrJPmZt6K26v00P13nPZP8/rbCJ4Dvuft4rYWt3GYLZEN7Pl/N3mPczC+iozu+T7TX/fPxYw8RfeAB+oj+zD8O/D/gfctQ008S/Vn1PPBc/HU78IvAL8br3A8cIzpC4GngI8tQ1/vi9/tO/N6V7VVdlwEH4u35AjC2TD/HVURhvbbqsbZsL6L/YN4EykTzmp8l2k/z58ArwJ8Bw/G6Y8B/rHruZ+LP2nHgF1pc03GiOdjKZ6xyVNi1wOGFfubLsL2+HH9+nicKrmvm1xbfv+L3t5V1xY9/qfK5qlp3WbbZAtnQls+XzlAVEcmglTwtIyIiS6RwFxHJIIW7iEgGKdxFRDJI4S4ikkEKdxGRDFK4i4hkkMJdRCSD/j8Il0pG5n5ZOAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# x[i] = x_end   \n",
    "for i in range(0, n+1):\n",
    "    eq.set_a(i, k)\n",
    "    eq.set_b(x_end*k)\n",
    "    eq.begin_new_row()    \n",
    "    \n",
    "\n",
    "solve_and_plot(eq)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Уже лучше. Мы потребовали сойтись к целевому состоянию $x_{end}$ как можно быстрее - мы это получили. Но в этом случае ускорение слишком большое (реальная механическая система не может же любое ускорение выдать). Мы хотим, чтобы система как можно быстее сошлась, но при этом с как можно меньшим ускорением.\n",
    "\n",
    "Поэтому потребуем, чтобы все время ускорение было как можно меньше.\n",
    "$$\n",
    "\\begin{cases}\n",
    "...\\\\\n",
    "u_0 = 0 \\\\\n",
    "u_1 = 0 \\\\\n",
    "...\\\\\n",
    "u_{n-1} = 0 \\\\\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Аналогично, домножим на небольшой коэффициент, чтобы ничего не ломалось."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/scipy/sparse/compressed.py:746: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.\n",
      "  SparseEfficiencyWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8XOV97/HPb0arLdmWLXnRAjZg4tgYsCNMSpYStmLfBm5aqOE2zUahaaEJN91IkxJK76stCTf33rS+oWQjpL0BSm5aN7FZSghZCsEGjNdghDFobNmWJdmWtS9P/zhH8mg8I43kmTmaM9/36zWvOcszMz8fjb9z5jnPnGPOOUREJFwiQRcgIiKZp3AXEQkhhbuISAgp3EVEQkjhLiISQgp3EZEQUriLiISQwl1EJIQU7iIiIVQU1AtXV1e7xYsXB/XyIiJ56aWXXjrqnKuZqF1g4b548WK2bt0a1MuLiOQlM3srnXbqlhERCSGFu4hICCncRURCSOEuIhJCCncRkRCaMNzN7JtmdsTMdqZYb2b2FTNrMrPtZrY682WKiMhkpLPn/hBw7Tjr1wJL/dttwFfPvCwRETkTE45zd879xMwWj9PkeuBh512v7wUzm2Nmi5xzLRmqUaYR5xyDw46BoWH6B4fp9+8HhsYuG/Dvh4Ydw84xPAxDzuGcY2gYb5lLvc45x7DzXs8BzuHfO78OcLi45X59/rLEmk9Nxy0f0yZ++fiXnsy3K1OaJVnG6QtH2lmyZcmeZMxjDDPvsWYQ8VeY2egyG5m3scsjZkTNiESMaMSfj8Qt8+cjo8sYXVYcjVBaHKG0KEppUYSSogilRd58cdRS1p1Lw8OOw529NLf30NzeTXNHN1cuW8DK+tlZfd1M/IipDmiOm4/5y04LdzO7DW/vnrPOOisDLy3pcs7R1T/EiZ4Bjvu3kekTvYOj86eWjbQZpG9wiIEh54X20HDehZsUJjO8wI9GKC2O+qHvBX9JUYSZpVGqZpQwd2YJVTNKmFdRMjo/cquaUUJJ0fgdHM45OroHRoO7ub3Hv+8m1tHDgY4e+oeGx9RVXVGaF+GeNufcg8CDAI2NjYqIDHPO0Xqyj72HTrL3cCevH+nktUOdvNXWzbGeAYaGU29yM6gsLWJWeTGz/ds51RXMKi+irDhKcTRCcdTbMyqJ2uh0cdT7zzMyXRw1v02E4qKIt8fl75mZcWo+4u29JVsX8ffeRvf0GNnz83b/xuwFxq2L30kbbR83Pzo95t9tKZaPv62nwx5hOlyST+JkH85udJ1Lsiy+XYpvQad9k/K+eZFk+eg3MX962DnvG57/DW7k297Q8NjpYf+b3cj04LBj0P+22Dc4TN/gEH2DcfMDQ/7ysev7BrxvlV19g+w+eIL27n6OdQ+k3IaVpUVUzRwb+jNKohw81kvMD/Gu/qExj6maUUzD3BksXzSLa1YsoKFqBg1zZ9BQVU5dVTmlRdGUr5cpmQj3A0BD3Hy9v0yyqL2r3wvww528driTvYe9QI9/k1bNKOb8BZVcs2Ihc2eeCu1ZZf59+an7ytIiIpH8CCxJX7IPofE/lwrzPTA4NMyxngE6uvpp6+qno6uf9u5+2k/6913e7UhnL79sOcHJvkEWzS6nYW457z5n3mhwN8z1QryiNLAzu4zKRAUbgTvM7BHgUuC4+tsz68CxHp57rZW9hzv920mOnuwbXV9ZVsT5CypZe8Eizl9QwTsWVLJ0QSXVFSV5s4cpEqSiaITqilKqK0pZGnQxGTJhuJvZd4HLgWoziwFfAIoBnHMPAJuAdUAT0A18PFvFFppX3u7g6z97kyd2HmJo2DGjJMrS+RV84B01nL+gkvMXVnL+ggoWzipTiIvIGOmMlrl5gvUOuD1jFRW4waFhntp9mG/87E1eequDyrIifve9S1h/SQOL581U14mIpCX4jiEBoLN3gMe2xvjWz98k1tHDWXNncM8Hl3NDY8O06L8Tkfyi1AhYrKObb//Hfh55sZnOvkEuWVzF5//Lcq5evoCo9tJFZIoU7gGJ708HWLdyEbe8dwkXN8wJuDIRCQOFew4NDTue2nWIryf0p3/kssXUzSkPujwRCRGFe47867YD3P/UazS399Awt5wvfHA5N6o/XUSyRMmSA0/sbOHOR7exsm42n1v3Tq5evlD96SKSVQr3LNvWfIw7H93GxQ1z+O6t76asOPs/OxYR0cU6sqi5vZvf/fYWaipL+dpHGhXsIpIzCvcsOd4zwCce2kL/4DDf+tglVFeUBl2SiBQQdctkwcDQMLf/08u8ebSLh29Zw3nzK4MuSUQKjMI9w5xzfP77O/lZ01G+dMOFXHZuddAliUgBUrdMhj3w3D4e3drMH15xHjc2Nkz8ABGRLFC4Z9APt7dw3xO/5LqLavnM1ecHXY6IFDCFe4a89FYH//2xbTSeXcUXb7hQp+AVkUAp3DPg7bZubnt4K4tml/GghjyKyDSgcD9Dx7sH+PhDLzI47PjWxy5h7sySoEsSEVG4n4n+wWE++Y8v8XZ7Nw/+zrs4p6Yi6JJERAANhZwy5xx//v0dPL+vjf+1/iIuPWde0CWJiIzSnvsUbXi2icdfinHnVUv50Kr6oMsRERlD4T4FG189yP1P7eVDq+r49JVhuVa6iISJwn2Stu5v54//+VXWLJnL3/7mSg15FJFpSeE+CfuPdnHrw1upm1POP3z4XZQWacijiExPCvc0DQ4Nc8u3twDwrY9dQpWGPIrINKbRMml68c123mjt4is3r2Jx9cygyxERGZf23NO0aWcL5cVRrn7ngqBLERGZkMI9DUPDjid2HuaKZfMpL1E/u4hMfwr3NGzd387Rk32sXbkw6FJERNKicE/D5p2HKC2K8IF3zA+6FBGRtCjcJzA87Ni8s4XL31HDzFIdfxaR/KBwn8ArzR0cPtHHupWLgi5FRCRtaYW7mV1rZq+ZWZOZ3ZVk/Vlm9qyZvWJm281sXeZLDcamHYcoiUa4Ypm6ZEQkf0wY7mYWBTYAa4HlwM1mtjyh2eeBx5xzq4CbgP+b6UKD4Jxj844W3n9+NZVlxUGXIyKStnT23NcATc65fc65fuAR4PqENg6Y5U/PBg5mrsTgbGs+xsHjvay9QF0yIpJf0jlCWAc0x83HgEsT2twDPGVmfwjMBK7KSHUB27zzEMVR4yr9cElE8kymDqjeDDzknKsH1gHfMbPTntvMbjOzrWa2tbW1NUMvnR3OOTbtaOE951Uze4a6ZEQkv6QT7geAhrj5en9ZvFuAxwCcc88DZUB14hM55x50zjU65xpramqmVnGO7DxwglhHD+vUJSMieSidcN8CLDWzJWZWgnfAdGNCm7eBKwHM7J144T69d80nsGlnC9GIcfVydcmISP6ZMNydc4PAHcCTwB68UTG7zOxeM7vOb/ZHwK1m9irwXeBjzjmXraKzbWSUzGXnztOpfUUkL6X1k0vn3CZgU8Kyu+OmdwPvyWxpwdnT0sn+tm5ue/+5QZciIjIl+oVqEpt3thAxuGaFumREJD8p3BM45/jhjhYuXTKP6orSoMsREZkShXuC14+cZF9rF+t0el8RyWMK9wSbdrRgBr+2QuEuIvlL4Z5g845DXHL2XObPKgu6FBGRKVO4x2k6cpLXDnfqiksikvcU7nGe2NkCwLUXKNxFJL8p3ONs2nGI1WfNYdHs8qBLERE5Iwp33/6jXexuOaErLolIKCjcfZt3HgLUJSMi4aBw923e2cJF9bOpr5oRdCkiImdM4Q40t3ezPXacteqSEZGQULgDT/hdMmvVJSMiIaFwxzt3+4raWZw9b2bQpYiIZETBh/vBYz288vYxjZIRkVAp+HBXl4yIhFHBh/vmnS0sW1jJOTUVQZciIpIxBR3uh0/0svWtDtbqItgiEjIFHe5P7jqEc+jc7SISOgUd7pt2tHDe/AqWLqgMuhQRkYwq2HBv7ezjxTfbWacDqSISQgUb7k/tPsSwQ79KFZFQKthw37zjEEuqZ7JsobpkRCR8CjLc27v6eX5fG2svWIiZBV2OiEjGFWS4P737EEPDTr9KFZHQKshw37TjEA1zy1lROyvoUkREsqLgwv149wA/bzrKugsWqUtGREKr4ML96T2HGRx2GiUjIqFWcOG+eUcLtbPLuKh+dtCliIhkTf6F+8FX4KdfntJDT/QO8NPXj7J2pbpkRCTc0gp3M7vWzF4zsyYzuytFm98ys91mtsvM/l9my4zz9gvwzF/CiZZJP/SVt4/RPzTMFcvmZ6EwEZHpo2iiBmYWBTYAVwMxYIuZbXTO7Y5rsxT4LPAe51yHmWUvPWtXefct22DW5PrNYx3dACyp1hWXRCTc0tlzXwM0Oef2Oef6gUeA6xPa3ApscM51ADjnjmS2zDgLV4JFvO6ZSYp19FAcNRbMKstCYSIi00c64V4HNMfNx/xl8c4Hzjezn5vZC2Z2baYKPE3JTKhZNuVwr51TTjSi/nYRCbdMHVAtApYClwM3A18zszmJjczsNjPbamZbW1tbp/5qtavg4DZwblIPi3V0U19VPvXXFRHJE+mE+wGgIW6+3l8WLwZsdM4NOOfeBPbihf0YzrkHnXONzrnGmpqaqdbshXvXEThxcFIPi3X0UD9nxtRfV0QkT6QT7luApWa2xMxKgJuAjQlt/gVvrx0zq8brptmXwTrHGjmoOomumd6BIVo7+7TnLiIFYcJwd84NAncATwJ7gMecc7vM7F4zu85v9iTQZma7gWeBP3HOtWWraBasgEjRpML9wLEeAOrnKtxFJPwmHAoJ4JzbBGxKWHZ33LQDPuPfsq+4HOa/c1LhHuvww71K3TIihW5gYIBYLEZvb2/QpaRUVlZGfX09xcXFU3p8WuE+LdWugj0/8A6qpvFr05Ex7uqWEZFYLEZlZSWLFy+elr9Wd87R1tZGLBZjyZIlU3qO/Dv9wIhFF0NPOxx7O63mI2Pc51dqjLtIoevt7WXevHnTMtgBzIx58+ad0TeL/A33SR5UbW7v1hh3ERk1XYN9xJnWl7/hvmAFRIrTDvdYR4+6ZESkYORvuBeVegE/mXDXGHcRKRD5G+6Q9i9VeweGOHqyjwYNgxSRaeCuu+5iw4YNo/P33HMP999/f0ZfI39Hy4AX7i99C9r3wbxzUzbTMEgRSeUv/20Xuw+eyOhzLq+dxRc+uCLl+vXr13PnnXdy++23A/DYY4/x5JNPZrSG/A938Lpmxg13DYMUkelj1apVHDlyhIMHD9La2kpVVRUNDQ0TP3AS8jvc578ToqVeuK+8IWUz7bmLSCrj7WFn04033sjjjz/OoUOHWL9+fcafP7/DPVrsnd+95dVxm50a416ao8JERMa3fv16br31Vo4ePcpzzz2X8efP7wOqcOqg6vBwyiaxjm7q5pQT0Rh3EZkmVqxYQWdnJ3V1dSxaNLmryqUjv/fcwQv3LV+D9jeg+rSzDAMjY9zVJSMi08uOHTuy9tzh2HOHcce76wdMIlJo8j/cq8+HovKU4T4yxl3hLiKFJP/DPVoEiy5MGe6nhkGqW0ZECkf+hzt4XTMtr8Lw0GmrmkeHQWrPXUQKR3jCfaAbju49bZXGuItIIQpPuEPSrplYRzcl0YjGuItIQQlHuM87D0oqUoR7D3VVGuMuIoUlHOEeicKii1KGu/rbRaTQhCPcweuaObQDhgbGLD7Q0a1wF5FpZf/+/VxwwQWj8/fffz/33HNPRl8j/3+hOqJ2FQz2QusvvfPNAD39Qxw92a+DqSKS2ua7vB3DTFq4Etb+bWafc5LCtecOY7pmDhzTqX5FpDCFZ8+9agmUzvZOIrb6I4DGuItIGgLYwy4qKmI47mSHvb29GX+N8Oy5RyJQO/agqsa4i8h0tGDBAo4cOUJbWxt9fX384Ac/yPhrhGfPHbyumRe+CoP9UFQyOsa9pkJj3EVk+iguLubuu+9mzZo11NXVsWzZsoy/RrjCfdHFMNQPR3ZD7cXE2jXGXUSmp0996lN86lOfytrzh6dbBk47qBrTMEgRKVDhCveqxVA2Jy7c9QMmESlM4Qp3M/+ye6/Q3T9IW5fGuItIYUor3M3sWjN7zcyazOyucdr9ppk5M2vMXImTVLsKjuzmYGsHoGGQIpKccy7oEsZ1pvVNGO5mFgU2AGuB5cDNZrY8SbtK4NPAL86oojNVuwqGBzm2fxugYZAicrqysjLa2tqmbcA752hra6OsrGzKz5HOaJk1QJNzbh+AmT0CXA/sTmj3V8B9wJ9MuZpM8A+qDsVeBi6gQXvuIpKgvr6eWCxGa2tr0KWkVFZWRn19/ZQfn0641wHNcfMx4NL4Bma2Gmhwzv3QzIIN99n1MKOastbtlBRdSLXGuItIguLiYpYsWRJ0GVl1xgdUzSwCfBn4ozTa3mZmW81sa9Y+Mf2DqtUndlM/R2PcRaQwpRPuB4CGuPl6f9mISuAC4Mdmth94N7Ax2UFV59yDzrlG51xjTU3N1KueSO0qFvbvZ8lsBbuIFKZ0wn0LsNTMlphZCXATsHFkpXPuuHOu2jm32Dm3GHgBuM45tzUrFaejdhVRhmksOzBxWxGREJow3J1zg8AdwJPAHuAx59wuM7vXzK7LdoFT0VPjnc99he0LuBIRkWCkdW4Z59wmYFPCsrtTtL38zMs6M7HBOcx2c1jctzfoUkREAhGuX6j6mju62T68hJrOxNGaIiKFIZThHuvoYcfwOZQda4K+k0GXIyKSc6EN9z2RczEcHNoedDkiIjkX0nDvpm2Wf4aEuCsziYgUinBdrMMX6+hhzrw6iNYp3EWkIIV0z90/j7t/+l8RkUITunDv6hukvavfD/eLoa0Jeo8HXZaISE6FLtwPHOsB/FP9jlx2r+XVACsSEcm90IV7rKMb8C/SsWjsNVVFRApF6A6oxjpG9tzLYWYZzDlL4S4iBSeEe+49lBZFqBk5j7sOqopIAQphuHdTV1WOmX+639pV0LEfutsDrUtEJJdCF+7N7T1jr5u66GLvXgdVRaSAhC7cYx3dXn/7iFo/3NU1IyIFJFThfrJvkI7uARri99zLq6BqicJdRApKqML9QPxImXi1q+DgtgAqEhEJRqjCfcwY93i1q+D429B1NICqRERyL2ThHvfr1Hgjv1TV3ruIFIiQhXs3pUURqitKxq5YdJF3r353ESkQIQt372yQo2PcR5TNgnlLFe4iUjBCGO4zkq/UL1VFpICELNy7Tz+YOqJ2FXQehM5DuS1KRCQAoQn3kTHu4+65gw6qikhBCE24pxzjPmLhSrCIumZEpCCEJtyb21OMcR9RWgHV71C4i0hBCE24n/oBU4puGfDOM3PwFXAuR1WJiAQjROHeQ1lxkjHu8WpXQdcROHEwd4WJiAQgVOFeXzXj9DHu8Rou9e53fi83RYmIBCQ84X5snGGQI2ovhnOvgJ99GXqP56YwEZEAhCfc/V+nTujKu6GnA/7j77JflIhIQEIR7p29Axwbb4x7vNpVsOJD8PwGOHkk+8WJiAQgrXA3s2vN7DUzazKzu5Ks/4yZ7Taz7Wb2jJmdnflSUztwbIIx7ok+8HkY7IOffCmLVYmIBGfCcDezKLABWAssB242s+UJzV4BGp1zFwKPA1/MdKHjibWnONVvKtXnwerfga3f8i6eLSISMunsua8Bmpxz+5xz/cAjwPXxDZxzzzrnuv3ZF4D6zJY5vpQX6RjPr/4ZRKLw7N9kqSoRkeCkE+51QHPcfMxflsotwOZkK8zsNjPbamZbW1tb069yAiNj3OfNHGeMe6JZtXDp78H2R+HwrozVIiIyHWT0gKqZfRhoBJJ2ZjvnHnTONTrnGmtqajL2ummNcU/mPXdC6Sx45q8yVouIyHSQTrgfABri5uv9ZWOY2VXA54DrnHN9mSkvPWmNcU9mxlx476dh72Z4+4XMFyYiEpB0wn0LsNTMlphZCXATsDG+gZmtAv4BL9hzPr6wuT3NMe7JXPpJqFgA/36PzjkjIqExYbg75waBO4AngT3AY865XWZ2r5ld5zf7ElAB/LOZbTOzjSmeLuNO9A5wvGeAhnRHyiQqmQm/+qfw9vPw+tOZLU5EJCBF6TRyzm0CNiUsuztu+qoM15W2U+dxn2K4A6z+KPzH38MzfwnnXQWRUPy2S0QKWN6nWGyii3SkI1oMV3weDu/UScVEJBRCEO5TGOOezIrfgAUr4dn/AYP9GahMRCQ4IQj3HsqLo8ydzBj3ZCIRuOoL3i9WX/52RmoTEQlKCMLdGwY56THuyZx3FZz9Hnjui9DfdebPJyISkBCE+xkMg0xkBld+wbta0wtfzcxziogEICThfgYjZRKddSmcvxZ+/hXobs/c84qI5FBeh/vIGPeM7bmPuPIvoO8E/Px/Z/Z5RURyJK/DPSNj3JNZsAIuXA+/+AddTFtE8lJeh3tGxrin8oHPwvAQPHdf5p9bRCTL8jrcm9szNMY9marF0PgJePk7cLQp888vIpJFeR3usY4eZpRkYIx7Ku//Yygq837YJCKSR/I83DM4xj2ZivnwK7fDru/DwVey8xoiIlmQ5+Ge4WGQyVx2B5TPhWfuze7riIhkUJ6H+xQv0jEZZbPhfX8Eb/wI9j2X3dcSEcmQvA334z0DnOgdzH64A1zyuzCrzjslsC7oISJ5IG/DPWtj3JMpLoPLPwsHXoJfPKCAF5FpL2/DPWOn+k3XRTfDOR+AJ+6Cxz8OPcdy87oiIlOQx+Gewz13gGgRfPh73onF9vwbPPBeeOv53Ly2iMgk5XW4zyiJUjWjOHcvGonC+z4Dn3gKIkXw0Dp49q9haDB3NYiIpCGPwz3LY9zHU/8u+ORPvfPPPHefF/Idb+W+DhGRFPI43HMwxn08pZXwoQfgN74OR/Z43TQ7Hg+uHhGROHkb7s25GOOejgtv9Pbia5bB926B7/8+9HUGXZWIFLi8DPfjPQN09g7SEOSee7yqxfDxzfCrfwbbH4EH3gexl4KuSkQKWF6Ge86HQaYjWgQf+HP42A9haAC+eQ389MveaYNFRHIsT8M9x8MgJ+Psy+D3fwbLft37RevD1+uCHyKSc3ke7tNozz1eeRXc+BBcvwEOvAxfvQz2/CDoqkSkgORpuHczsyTKnFyOcZ8sM1j1Yfi9n8Ccs+HR34Z/vAFe/Bq0vaFTGIhIVhUFXcBUjAyDDGSM+2RVnwe3PA0/vR9efQSanvaWzzkbzrsSzr0ClrzfO/ukiEiG5HG4T9MumWSKSryDrZd/Ftr3eacPfuNHsP0x2PpNsCjUN8K5ftjXrfZ+DSsiMkV5Gu7drFlcFXQZk2cG8871bmtu9UbVNL94Kux//Dfw47/29uLPudwL+nOvhDkNQVcuInkmrXA3s2uB/wNEga875/42YX0p8DDwLqANWO+c25/ZUj0jY9yn5UiZyYoWw+L3eLcr/wK62uDNH3tB3/Qj2P2vXrt5S+GsS2FWPcxaBJW1p+5nzPU+NEQkd/o6oXgmRKbvYcsJw93MosAG4GogBmwxs43Oud1xzW4BOpxz55nZTcB9wPpsFDwtx7hnysx5cMFvejfnoPU1f6/+Gdj7FHS1AgkHYqOlULkQZtV696PBv8hftggqFkBxuT4ERKZisA8O7/RGvh142buuw9G9cMcWqF4adHUppbPnvgZocs7tAzCzR4Drgfhwvx64x59+HPh7MzPnMj8kZFqPcc8kM5i/zLv9yh94y4YGoPMQdLZ4Y+dH7/1lLdth75Mw0J3k+SJQPMML+eJyf3rG2GUlMxPW+dPRIu8smJEi7/hApMg7JhCJWz5mPq6NRbwb5k9bknlLsd7i7vGmR7bNmGkS2pGwbszC5Ns6nXYTPibNx075edN6ghTPle7yFNu4UHYMhofg6Otw0A/xAy97wT7U762fWQN17/J2wEpmBlvrBNIJ9zqgOW4+Blyaqo1zbtDMjgPzgKOZKDLetB/jnk3RYq//fbw+eOeg9/jYD4CTR7zAH+g5dd/f5c/3QHebPx23rL+L074liIxK+NAd+RAfvSV8WMffEtuP7AxEi+Pui0/NJ11X5N0XlUFphRe0JRXeCf2STZfM9OajccOnnYPjMS/ED/p75Qe3Qb9/bqiSSqi9GN79+1C72gv12fV580GX0wOqZnYbcBvAWWedNaXnWLN4LnetXTa9x7gHyQzK53i3+e+c+vM4530dHej29maGB+Nu48y7oVPLhgbADQPOu3cubt4lzCeu9+dHPmBGvwS6038jkKpdsn/T6QvTbDfBY9J+7BSfN62Hxz/eTX65i59Otr3d6dt3ZNuP/O1Gb4l/y8TbyGOG/PfPgPd+GR707we86yQM9vrLBuPaDJyaH+zzT9SX5raLlvpBXwH93dDt739GimHhSrhovRfitau9Lpc8HrWWTrgfAOJ3Fev9ZcnaxMysCJiNd2B1DOfcg8CDAI2NjVN6J6+sn83Keo0Jzzoz79qxxWVBVyIyPuf8b5snvaDv7/Km+7vGzved9Jf76yJRWHSxN/R4wQVQVBr0vySj0gn3LcBSM1uCF+I3Af8toc1G4KPA88ANwI+y0d8uInIaMyiZ4d0q5gddzbQxYbj7feh3AE/iDYX8pnNul5ndC2x1zm0EvgF8x8yagHa8DwAREQlIWn3uzrlNwKaEZXfHTfcCN2a2NBERmarpOwJfRESmTOEuIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhZEH91sjMWoG3pvjwarJw3poMUF2To7omb7rWprom50zqOts5VzNRo8DC/UyY2VbnXGPQdSRSXZOjuiZvutamuiYnF3WpW0ZEJIQU7iIiIZSv4f5g0AWkoLomR3VN3nStTXVNTtbryss+dxERGV++7rmLiMg4pnW4m9m1ZvaamTWZ2V1J1pea2aP++l+Y2eIc1NRgZs+a2W4z22Vmn07S5nIzO25m2/zb3cmeKwu17TezHf5rbk2y3szsK/722m5mq3NQ0zvitsM2MzthZncmtMnZ9jKzb5rZETPbGbdsrpk9bWav+/dVKR77Ub/N62b20SzX9CUz+6X/d/q+mc1J8dhx/+ZZqu0eMzsQ9/dal+Kx4/7/zUJdj8bVtN/MtqV4bFa2WapsCOz95Zyblje8c8e/AZwDlACvAssT2vwB8IBFOGyAAAAEBElEQVQ/fRPwaA7qWgSs9qcrgb1J6roc+EEA22w/UD3O+nXAZrwrH78b+EUAf9NDeON0A9lewPuB1cDOuGVfBO7yp+8C7kvyuLnAPv++yp+uymJN1wBF/vR9yWpK52+epdruAf44jb/1uP9/M11Xwvr/Cdydy22WKhuCen9N5z33NUCTc26fc64feAS4PqHN9cC3/enHgSvNsnv1Wudci3PuZX+6E9iDd4HwfHA98LDzvADMMbNFOXz9K4E3nHNT/fHaGXPO/QTvgjLx4t9H3wb+a5KH/hrwtHOu3TnXATwNXJutmpxzTznnBv3ZF/Aub5lzKbZXOtL5/5uVuvwM+C3gu5l6vTRrSpUNgby/pnO41wHNcfMxTg/R0Tb+f4TjwLycVAf43UCrgF8kWf0rZvaqmW02sxU5KskBT5nZS+ZdjDxROts0m24i9X+4ILbXiAXOuRZ/+hCwIEmbILfdJ/C+cSUz0d88W+7wu4y+maKbIcjt9T7gsHPu9RTrs77NErIhkPfXdA73ac3MKoDvAXc6504krH4Zr+vhIuDvgH/JUVnvdc6tBtYCt5vZ+3P0uhMysxLgOuCfk6wOanudxnnfkafNEDIz+xwwCPxTiiZB/M2/CpwLXAy04HWBTCc3M/5ee1a32XjZkMv313QO9wNAQ9x8vb8saRszKwJmA23ZLszMivH+eP/knPv/ieudcyeccyf96U1AsZlVZ7su59wB//4I8H28r8bx0tmm2bIWeNk5dzhxRVDbK87hke4p//5IkjY533Zm9jHg14Hf9kPhNGn8zTPOOXfYOTfknBsGvpbiNQN5r/k58BvAo6naZHObpciGQN5f0znctwBLzWyJv9d3E7Axoc1GYOSo8g3Aj1L9J8gUvz/vG8Ae59yXU7RZONL3b2Zr8LZzVj90zGymmVWOTOMdkNuZ0Gwj8BHzvBs4Hvd1MdtS7k0Fsb0SxL+PPgr8a5I2TwLXmFmV3w1xjb8sK8zsWuBPgeucc90p2qTzN89GbfHHaT6U4jXT+f+bDVcBv3TOxZKtzOY2Gycbgnl/ZfqIcSZveKM79uIddf+cv+xevDc8QBne1/wm4EXgnBzU9F68r1XbgW3+bR3wSeCTfps7gF14IwReAC7LQV3n+K/3qv/aI9srvi4DNvjbcwfQmKO/40y8sJ4dtyyQ7YX3AdMCDOD1a96Cd5zmGeB14N+BuX7bRuDrcY/9hP9eawI+nuWamvD6YEfeYyOjwmqBTeP9zXOwvb7jv3+24wXXosTa/PnT/v9msy5/+UMj76u4tjnZZuNkQyDvL/1CVUQkhKZzt4yIiEyRwl1EJIQU7iIiIaRwFxEJIYW7iEgIKdxFREJI4S4iEkIKdxGREPpPQwZN2JiizrgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# u[i] = 0   \n",
    "for i in range(0, n):\n",
    "    eq.set_a(i+n+1, k)\n",
    "    eq.set_b(0)\n",
    "    eq.begin_new_row()\n",
    "\n",
    "solve_and_plot(eq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вот теперь хорошо.\n",
    "\n",
    "### 2. Ищем коэффициент\n",
    "Ищем коэффициент $F$ в уравнении, как описывалось ранее для линейной регресии:\n",
    "$$\n",
    "u = -Fx\n",
    "$$\n",
    "\n",
    "**КАЗАЛОСЬ БЫ, ВСЕ ОК**\n",
    "\n",
    "Но когда я решил в таком виде, получилась какая-то хрень. Можно посмотерть на графиики и убедиться, что они не совсем на коэффициент отличаются. Поэтмому я сделал так и все заработало:\n",
    "$$\n",
    "u = kx + b\n",
    "$$\n",
    "\n",
    "Хз почему."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/scipy/sparse/compressed.py:746: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.\n",
      "  SparseEfficiencyWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8XOV97/HPbxZZkiXLsiXbsiVjg83mBTDGbAmFCyRmM0mbxBDShCaB5jaUZm3JTS4h3DY3JNDepNC0hFKWLODQJLWJKVACSUMwYBZbXgDLYEuyZEuWZFv2aJnluX/MSB7LkjWWZuZoRt/36zUvneU5c346Gn3nzJlznmPOOUREJL/4vC5ARETST+EuIpKHFO4iInlI4S4ikocU7iIieUjhLiKShxTuIiJ5SOEuIpKHFO4iInko4NWKKyoq3Jw5c7xavYhITnrttdf2Oucqh2vnWbjPmTOH9evXe7V6EZGcZGY7U2mnwzIiInlI4S4ikocU7iIieUjhLiKShxTuIiJ5aNhwN7MHzazFzDYNMd/M7AdmVmdmG81sSfrLFBGR45HKnvtDwPJjzL8CmJ943Az8cPRliYjIaAx7nrtz7ndmNucYTa4FHnHx+/WtM7PJZlblnGtOU40yxjnnCEcd4WiMcDRGbzQWH48MGI/GCEfi485BzDmiMUfMxZ8j6pKGE9NjzhFLGu6b5wDnwCXVQP801z+v7y6SfdOS25I0P95m4O+VPG/wZQbdHsNvsOFaZJ7ZsWenuLgltTw87ehV2ID1mcWXjf88cryvva9/3uF2WHyNfp/hT0z3+wy/z/DZwJ/gS7RLnu73QYHfT0HAR0HAx4TEz/5hv++oenNROi5imgU0JI03JqYdFe5mdjPxvXtmz56dhlXLSDnnONQbpbM7zIGuCAe6w0cMH+gK09ndN5z42R2hsytMZ0+kP6jDUUdvNOb1r5NzvMyOsfDeMtYV+AcJ/YCPCQE/EwI+JhUFKS8uoLw4SPnEgqOHJwaZXFRAQcC7rzWzeoWqc+5+4H6ApUuX6iWWYdGYo7EjRF3LQba1HKQu8djZdogD3RGisWP/CQqDPkoLg0wqDDCpKMjkoiCzpxRTMsHPhICfoN8I+n0E/fEX/hHjfh/BwIBxf7xNwO87vGdl8T0qny9pODHd74vvmQ0c9iX23vr29ujf20va0+PovUGGmkby/AF7mEfMG7rdeJL8Kal/2sB5R7Tva3Pkckd/yjr8iYy+T2oDpjvijfs+5cVi8U9y0cQnvGjfp75Y4pNh3/T+NvFlw5EYPZEYvdEoPeH4p8nexLSeSN9wtH9ab9K0nkiMPQe6eav5AB2hMF3h6JDbqmRCgPKJfW8Eh98Arl5cxdknTBnlX+LY0hHuu4CapPHqxDTJkp5IlB17Q/3hXdd6kG17Onlv7yF6Iof3qqeVTmDetBKuWFTFlOICShOhPakwyKSiAJMKg/3TSgsDTAj4PfytZKzqe2Mb/P1t/L3pdYejdIR66TgUjv8M9dJxqJeOUPio4Xf3HmTfoTCnzZiUE+G+GrjFzB4DzgX263h75uwL9fL82y1s2xPfG9/ecpCd7aH+vXAzqC4vYl5lCRedXMm8yhJOmlbCvGkllBUFPa5eJP8UBv1UlRVRVVaU8jIuC8fGhg13M/sZcDFQYWaNwDeBIIBz7p+BtcCVQB0QAv4sU8WOZ1uaDvDISzv41Zu76A7HCPiMuRUTOWVGKVctrmJeIsBPrCihqEB73CJjWTYO66Vytsz1w8x3wOfTVpH0C0dj/Oem3Tzy0g5e3dFBUdDPh8+q5vplNZxWNYmgX9egicjgPOvyV4bW0tnNz15u4Ccv76Sls4cTphbzjatO46Nn11BWrEMrIjI8hfsY4Zzj9foOHv7DTp7a1Ew46rj4lEruOn8Of3RyJT7f+PuiSkRGTuHuse5wlNVvNvHwSzvY3HSA0sIAnzx/Dn963gnMqZjodXkikqMU7h5paA/x43U7eXx9A/tCYU6ZXsrffXghHzpzFhMn6M8iIqOjFMmyhvYQ31qzmefeasFnxgcXTOeT58/h3LlTxvWFMSKSXgr3LGrt7OGGB16mI9TL5y+exw3nzT6uc2NFRFKlcM+Szu4wN/7bK7R0dvPTm85jyexyr0sSkTymcM+CnkiUz/34Nd7a3ckDn1yqYBeRjNNVMBkWizm+vGoDL9a1cdefLOaSU6d5XZKIjAMK9wxyznHnk1t4cmMzt11xKh85u9rrkkRknFC4Z9A/vbCdh/6wg09fOJc/v+hEr8sRkXFE4Z4hq9Y38L2n32bFGTP5xlWn6TRHEckqhXsGPLd1D1/7RS3vn1/B3R89Q10HiEjWKdzT7LWd7Xz+p69zetUkfviJsz29zZaIjF9KnjTatqeTTz+0nhmTCvm3PzuHEnUjICIeUbinSdO+Lj754CsE/T4e+fS5VJRM8LokERnHFO5psC/Uy6cefIXO7ggPf/ocZk8t9rokERnnFO6j1NUb5bMPr2dnW4j7P3k2C2aWeV2SiIi6HxiNSDTGX/7sdV6r7+De65dwwUkVXpckIgJoz33EnHN8/Zeb+K+tLXxrxQKuWlzldUkiIv0U7iN0zzPv8Pj6Bv7yf8zjk+fP8bocEZEjKNxH4KEX3+Pe5+u47pwavnT5yV6XIyJyFIX7cXqxbi/fenILl502nb/90EJ1KyAiY5LC/Tg9+tJOpk4s4N6Pn0XAr80nImOT0uk4dHaH+c3bLVy1qIrCoN/rckREhqRwPw7PbN5DbyTGijNnel2KiMgxKdyPw5qNTcyaXMRZNbpNnoiMbQr3FLUf6uX32/Zy9RlV6sJXRMY8hXuKntrUTCTmuGaxDsmIyNincE/Rmg1NnFg5kQUzJ3ldiojIsFIKdzNbbmZvm1mdmd02yPzZZva8mb1hZhvN7Mr0l+qd3fu7efm9dq5ZPFPntYtIThg23M3MD9wHXAGcDlxvZqcPaPYNYJVz7izgOuCf0l2ol35d24xzcM0ZOiQjIrkhlT33ZUCdc+5d51wv8Bhw7YA2Dug7XlEGNKWvRO+t3tDE6VWTmDetxOtSRERSkkq4zwIaksYbE9OS3QF8wswagbXAX6alujGgvi3EhoZ92msXkZySri9Urwcecs5VA1cCj5rZUc9tZjeb2XozW9/a2pqmVWfWmo3xDyFXq0tfEckhqYT7LqAmabw6MS3ZZ4BVAM65l4BC4Kg7Vzjn7nfOLXXOLa2srBxZxVm2ZkMTS2ZPpmaKbp0nIrkjlXB/FZhvZnPNrID4F6arB7SpBy4FMLPTiId7buyaH8M7ezp5a3cnK3RIRkRyzLDh7pyLALcATwNbiZ8Vs9nM7jSzFYlmXwZuMrMNwM+AG51zLlNFZ8uaDU34DK7UIRkRyTEp3UPVObeW+BelydNuTxreAlyY3tK85ZxjzYYmzj9pKtNKC70uR0TkuOgK1SHU7trPjraQuhsQkZykcB/Cmg1NBP3G8oUzvC5FROS4KdwHEYs5ntzYzEXzK5lcXOB1OSIix03hPoj1Ozto3t+tC5dEJGcp3AexZkMThUEfl58+3etSRERGROE+QCQaY21tM5eeOp2JE1I6mUhEZMxRuA/wh+1ttB3q1SEZEclpCvcBVm9oomRCgItPyY3uEUREBqNwT9ITifL0pt18YMF0CoN+r8sRERkxhXuS377dSmdPRH3JiEjOU7gnWb2hifLiIBfOO6pDSxGRnKJwTwj1RnhuawtXLqoi6NdmEZHcphRLeHbLHrrCUZ0lIyJ5QeGesGZDM9MnTWDZnClelyIiMmoKd2B/KMxv32nh6sUz8fnM63JEREZN4Q48vXk34ajTWTIikjcU7sRvgj17SjGLq8u8LkVEJC3Gfbi3dvbwYt1erjmjCjMdkhGR/DDuw/2pTc3EHKw4Y5bXpYiIpM24D/fVbzZx8vQSTplR6nUpIiJpM67Dfde+Ltbv7NAXqSKSd8Z1uD+5oQmAq3UTbBHJM+M63NdsbOKM6jLmVEz0uhQRkbQat+H+butBNu06oO4GRCQvjdtwX7OhGTO4anGV16WIiKTduAx35xyrN+zinDlTqCor8rocEZG0G5fhvrW5k+2th3RIRkTy1rgM9zUbm/D7jCsXzvC6FBGRjBh34e6cY82GJi6cV8HUkglelyMikhE5F+5vblnFD/7jBpxzI1r+jYZ9NHZ06cIlEclrgVQamdly4PuAH3jAOfedQdp8DLgDcMAG59zH01hnv807n+dH+zbysT0bmDHjzONe/qXtbQBcftr0dJcmImNYOBymsbGR7u5ur0tJSWFhIdXV1QSDwREtP2y4m5kfuA+4HGgEXjWz1c65LUlt5gNfAy50znWY2bQRVZOCxbP/CFp+T+32p0YU7g3tISpKJlBWPLINJiK5qbGxkdLSUubMmTPme4B1ztHW1kZjYyNz584d0XOkclhmGVDnnHvXOdcLPAZcO6DNTcB9zrmORGEtI6omBaecdCVB56htfnVEy9e3h5g9Rac/iow33d3dTJ06dcwHO4CZMXXq1FF9ykgl3GcBDUnjjYlpyU4GTjazF81sXeIwTkYUFE7iVBdk48H6ES0fD/fiNFclIrkgF4K9z2hrTdcXqgFgPnAxcD3wIzObPLCRmd1sZuvNbH1ra+uIV7aoeCZbXDeRSM9xLReOxmja10WNwl1E8lwq4b4LqEkar05MS9YIrHbOhZ1z7wHvEA/7Izjn7nfOLXXOLa2srBxpzSycdiZdPmP7jt8c13LN+7qJORTuIpL3Ugn3V4H5ZjbXzAqA64DVA9r8ivheO2ZWQfwwzbtprPMIi0/8AACbdjx3XMvVt4cAdFhGRLLutttu47777usfv+OOO7j77rsztr5hw905FwFuAZ4GtgKrnHObzexOM1uRaPY00GZmW4Dnga8659oyVfTs6vcxKRajdu/G41pO4S4iXlm5ciWrVq3qH1+1ahUrV67M2PpSOs/dObcWWDtg2u1Jww74UuKRceb3s8g3kY1de45rufr2EEG/MX1SYYYqE5Fc8K01m9nSdCCtz3n6zEl885oFQ84/66yzaGlpoampidbWVsrLy6mpqRmy/WilFO5j0cJJc/nR/s2EDrVSPDG14/cN7SGqy4vx+3LnG3MRyR8f/ehHeeKJJ9i9e3dG99ohh8N9cdW5xA5sYXPdrznnjBtTWqahI6QvU0XkmHvYmbRy5Upuuukm9u7dy29/+9uMrivn+pbps3D+NQBsavjvlJfRBUwi4qUFCxbQ2dnJrFmzqKrK7I2CcnbPfcrU+cyKQm3HOym1398VZl8orC9TRcRTtbW1WVlPzu65AywuKKc23JFS24bEmTI15Qp3Ecl/OR3uC8tPZbffaG3ZMmzbxo5EuGvPXUTGgZwO98U1FwGwse7Xw7btP8d9qsJdRPJfTof7qfOuJOAcm3a/Mmzb+vYQk4uDTCpUV78ikv9yOtwLi6cw3wWo7dwxbNv69i4dbxeRcSOnwx1gcdEMNsW6iEUjx2zXoK5+RWQcyflwX1h5Bod8xns7XxiyTTTm2NWhrn5FZPzI+XBfPPdyADa+9+yQbfYc6KY3GtOeu4iMGzkf7nNO+CNKYo5NrUP3EKneIEVkLNixYwcLFy7sH7/77ru54447MrKunA93nz/IAl8RtV3NQ7bpC/cadT0gIuNEznY/kGxxyRwe7NxKd1cHhUXlR81vbA/hM5g5WeEuIsBTt8HuNHcDMGMRXPGd9D7nKOT8njvAwqpziJqxddvgFzPVt4eYObmIoD8vfl0RyVGBQIBYLNY/3t3dnbl1ZeyZs2jRvKtg26PUNv6OsxZ/4qj59ToNUkSSebSHPX36dFpaWmhra6OkpIQnn3yS5cuXZ2RdeRHuldMWMCPqqG1/a9D59e1dXHbatCxXJSJypGAwyO23386yZcuYNWsWp556asbWlRfhDrAoOJna3qN7iOzqjbL3YI/OcReRMeHWW2/l1ltvzfh68uYg9OLyU9jlh7a2bUdMb1BvkCIyDuVNuC+sfj8Am7atOWJ6fZvOcReR8Sdvwv30+Vfhc47a5iN7iNQFTCIyHuVNuBdPrGSe81N74N0jpte3h5hY4Ke8WF39isj4kTfhDrCocDq1sRCxWLR/WmNHiJopxZiZh5WJiGRXXoX74orFdPqM+vrf90/TOe4iMh7lVbgvnHspALXvPQOAc07hLiJj3o033sgLL7yQ1ufMq3A/ac6lFMUctS1vAtB6sIfucEz3TRWRcSdvLmIC8AcKWGCF1IaaAGho7wLQ7fVEZEzYsWMHV199NZs2bQLiXf4ePHiQsrIyCgoK0rquvAp3gEWls3m08x16uw/Q0K4LmERk7Pv+97+f9ufMv3CfvpTIwW28tf0p6tuXAFBdrq5+ReSwu165i7eG6ItqpE6dcip/s+xv0vqco5HSMXczW25mb5tZnZnddox2f2JmzsyWpq/E47No3pUA1Na/QH17iBmTCikM+r0qR0Sk35jq8tfM/MB9wOVAI/Cqma12zm0Z0K4U+Cvg5UwUmqoZM85kWtRR276V+oMh3X1JRI7i1R52Nrv8TWXPfRlQ55x71znXCzwGXDtIu/8D3AVk7q0oRQuDZdT2tNHYHtLxdhEZM5K7/L388ss97/J3FtCQNN4InJvcwMyWADXOuV+b2VfTWN+ILJo8n9+0v4b/YD2zp9R4XY6ISL+c6fLXzHzA3wNfTqHtzWa23szWt7a2jnbVQ1o060IA5ha9pguYRGRcSiXcdwHJu7/ViWl9SoGFwAtmtgM4D1g92Jeqzrn7nXNLnXNLKysrR171MBbMvxpzjtKibQp3ERmXUgn3V4H5ZjbXzAqA64DVfTOdc/udcxXOuTnOuTnAOmCFc259RipOQUlpFXOiRk9Ri465i8i4NGy4O+ciwC3A08BWYJVzbrOZ3WlmKzJd4EjVuHJ2FfZSoa5+RSTBOed1CSkbba0pXcTknFsLrB0w7fYh2l48qorSpCg6h/3BDpqaX6am5gKvyxERjxUWFtLW1sbUqVPHfBfgzjna2tooLCwc8XPk3RWqfVq7FkLhG9Ruf1rhLiJUV1fT2NhIJk/mSKfCwkKqq6tHvHxehrtzjjc7TmNCmaO25XWu9LogEfFcMBhk7ty5XpeRNXnV5W+f/V1h9vf4mBcNUnto1/ALiIjkmbwM976bYs+fMJOt1ku455DHFYmIZFdeh/vplUvoNeOdd//T44pERLIrr8P9nFOvAqB25wseViMikn15Ge4N7V1MnVjASbOXMSXmqG3b7HVJIiJZlafhHu8N0nw+FvknUduz1+uSRESyKi/Dvb491N+nzKKyk3jP7ziwv97jqkREsifvwj0SjbFrX9fhcJ8Vv4Bp87YnvSxLRCSr8i7cm/d3E425/jswLTz5GgBqd73kZVkiIlmVd+HekDhTpq83yEmTqpkTNWr313lZlohIVuVduPedBpncj/uiwgpqI524pBvTiojks7wM94DPqCo7fGPsRVMX0OY3mptf87AyEZHsyctwry4vwu873KXnotkXA1D77lMeVSUikl15F+4NHV1H3X3plJOuoMA5andrz11Exof8C/fEBUzJggXFnEoBtYcaPapKRCS78ircO7vDtB/qHfSm2IuKZ7HV9RAJd3lQmYhIduVVuDe0x4N70HCfvoQun7H93WeyXZaISNblVbj3nQZZUz5IuJ/4QQA27nw+qzWJiHghr8K9sePoc9z71Mw6j8kxR+3eTdkuS0Qk6/Iq3OvbQ0wqDFBWHDxqnvl8LPSXUNvd4kFlIiLZlXfhPnvq0XvtfRZNOontvhiHOpuzWJWISPblX7gPckimz6JZ5+HM2KIeIkUkz+VNuMdijsaOrkG/TO2zcP7VAGxs/H22yhIR8UTehHtLZw+9kdhRFzAlK588l5oY1O5TD5Eikt/yJtwH6w1yMAsLplIb2Q/OZaMsERFPjLtwP7vyLFr8xvo3HshGWSIinsibcG9oD2EGMycXHbPdivd/k2kxuPvNe4lFw1mqTkQku/Iq3GeWFVEQOPavVFQ0mVtP/BCb/TGe+u3tWapORCS78ibc69tD/fdNHc4177+D01yA7+9YQ/ehvRmuTEQk+1IKdzNbbmZvm1mdmd02yPwvmdkWM9toZs+Z2QnpL/XYhjvHPZnP5+crZ32BZr/x42e/kOHKRESyb9hwNzM/cB9wBXA6cL2ZnT6g2RvAUufcYuAJ4LvpLvRYusNRWjp7Ug53gGVnfIqLfWU80PEGba1bMlidiEj2pbLnvgyoc86965zrBR4Drk1u4Jx73jkXSoyuA6rTW+ax9XUYdqxz3AfzxYv+L91m/PC/tPcuIvkllXCfBTQkjTcmpg3lM8CgNys1s5vNbL2ZrW9tbU29ymH0d/V7nOF+4gnv56MT5/JETxPvbtP9VUUkf6T1C1Uz+wSwFPjeYPOdc/c755Y655ZWVlambb31bamd4z6Y/3nZDyhycM8fvqkLm0Qkb6QS7ruAmqTx6sS0I5jZZcDXgRXOuZ70lJea+vYuigv8TJ1YcNzLTimfy03TLuB3dLHu1XszUJ2ISPalEu6vAvPNbK6ZFQDXAauTG5jZWcC/EA/2rHeY3tARoqa8GDMb0fI3XHYPM2Nw96YfEQ13p7k6EZHsGzbcnXMR4BbgaWArsMo5t9nM7jSzFYlm3wNKgJ+b2ZtmtnqIp8uIhvbQcR9vTzZhQil/Nf863vY7Vj//tTRWJiLijUAqjZxza4G1A6bdnjR8WZrrSplzjvr2EBecVDGq57nigq/x47pf8I+Nz/DBg7spLpmRpgpFRLIv569QbTvUS6g3yuwUr04divl8fHXZX9Pq9/HwM7emqToREW/kfLj39wZ5jNvrpeqs01dyeWAK/7Z/M63Nb4z6+UREvJLz4d7Qd477Me7AdDy+eMn3CJtx72++nJbnExHxQt6Ee3Wawr1m5jKuLz2FX4ZbeHvrr9LynCIi2Zbz4V7fHmJa6QSKCvxpe84/v/wHlDq4Z93f4mKxtD2viEi25EW4j+TK1GMpmzSLz828hJd8Pfx+3T1pfW4RkWzI+XBvaO9Ke7gDXHfJd5kdM+556xEivaHhFxARGUNyOtx7IzGa93dRnYFwDxYU8cXTbmS7H37x3FfT/vwiIpmU0+HetK+LmBtZh2GpuPTcL7KEQu5rfoGD+xuGX0BEZIzI6XDvP8c9Q+FuZnzlvG/Q7vfx4LO6sElEcofCfRiLTrmWK4LTeKRzG7sb12VsPSIi6ZTT4d7QEaIg4GNa6YSMrucLl/4DzuAHz+vYu4jkhtwO9/YQ1eVF+Hwj6+o3VTOnL+YTZQtZE9vH5tqfZnRdIiLpkNPhnolz3Ify2ct/QHnMcff67+Gi0aysU0RkpHI73NuyF+6lJdP4i5rlrPdFeP7Fb2dlnSIiI5Wz4b4/FOZAdyRr4Q7wkYu/zdyYj3vqHqelQV+uisjYlbPh3tCR3g7DUhEIFPD1s79Ciw8+/OxneOqZL4H6nhGRMShnwz0bp0EO5tzFf8rPP/AwcwIT+evmZ/nKI+fTsXtjVmsQERlOzod7zSjvwDQSc2aezcM3vMit09/Pcxziw2uv57cvfBOcy3otIiKDyelwnzKxgNLCoCfrD/iD3LT8n3jsknuZ4pvALTt/wTcfvYiD7ds9qUdEJFnOhntDe4ia8uzvtQ90ygkX89gNf+AzU5bwq1gHf/KrFbz6kroJFhFv5Xa4Z/l4+1AKgoV84ZqHefiCbxPwBfn0Ow9x108uo7uz2evSRGScyslwj8YcjR2Z6cd9NM48eQU/v/6/ua70FH4c2cPHVl3Optcf8LosERmHcjLcm/d3EYm5MRfuAMUTSvn6Hz/Bv5z9Nxzy+fjExv/HvY9fTbirw+vSRGQcyclwb2jvArJ/GuTxuGDhJ/jlyt9wZfFs/qV7Jzf89CLqNq/yuiwRGSdyNNz7ToMcu+EOMKm4gm9/bC3/sPB/stsHK1+5k4f+fSXRnoNelyYieS4nw72+PYTfZ1SVFXpdSkouO/sv+MUf/5oLJ0zjnoNb+NMfn8+jv7yedzY8glPQi0gGBLwuYCTq20PMmlxEwJ87700VZbP5/vXPseal73D/O6v47oFN8OYmprz2HZZZCedWLObc+SuombccAgVelysiOc6cR1dVLl261K1fv35Ey37ovhcpmRDgx589N81VZU9zx3bWbXmcl3f9npdDjey1+N9hViTKucFyzp12NstO/QgVs98Hvtx5ExORzDKz15xzS4drl5N77o0dIS4/fbrXZYxKVflJfPjC/8WHAecc7+15g3Vbf87Lu1/h2Z4WfrHnedjzPPPCUc4tnMG5M89j6WkrKZ2xGCyzNycRkdyXUrib2XLg+4AfeMA5950B8ycAjwBnA23ASufcjvSWGneoJ8Leg71j/svU42FmnDhjCSfOWMLHgWgsytaG37Hu7V/ycusbPBFu4ScNa/DXr2ZBBM4unklV8XQqS6qonDSbaZNPpGLqyQTLaiCQ2VsOisgoRXrBxSCY2e8Mhw13M/MD9wGXA43Aq2a22jm3JanZZ4AO59w8M7sOuAtYmYmC+7r6HcunQY6W3+dn4QmXsPCES/gs0BPpZsP2/2Rd3Rpead/Mo71NRMLNsP9N2HV4ufJolIoYTLMCKgPFVBaUUVlUEX8TKK1hWvlJVEydT7B0FgSL9AlAJNMiPbBnM7GmN+hqeo3Q7lpC7duYfNnfUXbOZzK66lT23JcBdc65dwHM7DHgWiA53K8F7kgMPwHca2bmMnBAv74t/8N9oAmBQpad8iGWnfIhIL5n39HVTmvHNlo76mjdv5OWzl3s7WqlpaeDveGDbIsdoq3nANHexvibQJJJ0SjFzlHooBCjEB+F5j/88AUp9AUp8hcwwVdAYaCQQn8hRYEiJgSKKAwW4fcF8PuC8Z8WwO8PEvAF8Pn6hgvw9U33B/H7CvD7g/j9QXy+AD4LYGaYz4/hw8yHz+fHMMwXn+czf7yNBTCfJdr5MZ8PSLwx9b9B2ZHDQ85jwLwjJg7+Bzietsdc5jiWH86o35iTlj/iuY61nQZbxgYfztCOg3MOh8M5R4xY/3jMxYdjLnbk+IA2fdMjLkI0FiUSi/QPh2Nhoi4+LRqLEnGCnKOGAAAIIklEQVSR+PzEo29eOBamO9JNKBKKP8IhuiJdhHo7CR1qJdTVRqjnAKFIFyEXpsuMrr7vzQqBmZX8b7efj2VkCx2WSrjPAhqSxhuBgd9k9rdxzkXMbD8wFdibjiKTNXSM/QuYMs3v81MxsZKKiZWcVn3BkO0OvwnU9b8JtB7cRVvXXrqjPfFHLNz/2BeL0k2UbtdLd9TFHwb0ZO93S5Ul9huS43xgnPTPc2Ac3T7ldaU47Xh4+ZlptHtcA5cf7vmObh//7Z3F57lB2iZPd2P4E2ahBSjCKI5GKIr0UhyLUexiTLEgxYXlFE+soLhkJsWTT6C4pIqiYDHFwWIWVyzOeG1Z/ULVzG4GbgaYPXv2iJ5j2ZwpfO2KUykr8qar31xy5JvA+SN6jpiL0RPtoSfSQ3c4RFfPfnp6DxCN9hKNholGe4lEw8RiYaLRMBEXnxbr29uJhonGwkRjkfjDRYhGIzhiuAF7W67/Z3waA/bISNpjg8Pd57ukKOj7sOiSxvvDY8C8ZINNG7LtKD+Qjv4DrWN0bw8uxaXdoIPx5e2IGUc8n+ubf3ghS0wf+ITm3OE3Zuf634T7hg+/iR897ItF48u7GD4XxeccFosmxmOHh2Px+ZaY5otF8Ud7CUR68ANB5wg4hx8IOEeg76cDP/GfAfMRCBYTCBbj9wUo3reLolgEP0BxBcw8E2afCVVnxIfLajw/7JlKuO8CapLGqzniSO8RbRrNLACUEf9i9QjOufuB+yF+KuRICl5UXcai6rKRLCoj4DMfRYEiigJFUDgZSmd6XZJIesSiEA5B76HE4+Aww6H4cLgLFn4Mqs6MB/mkWZ4H+WBSCfdXgflmNpd4iF8HfHxAm9XAp4CXgI8Av8nE8XYRkbTx+WFCafyRh4YN98Qx9FuAp4mfCvmgc26zmd0JrHfOrQb+FXjUzOqAduJvACIi4pGUjrk759YCawdMuz1puBv4aHpLExGRkdJ17SIieUjhLiKShxTuIiJ5SOEuIpKHFO4iInlI4S4ikoc8u1mHmbUCO0e4eAUZ6LcmjVTf6Ki+0RvrNaq+kTvBOVc5XCPPwn00zGx9Knci8YrqGx3VN3pjvUbVl3k6LCMikocU7iIieShXw/1+rwsYhuobHdU3emO9RtWXYTl5zF1ERI4tV/fcRUTkGMZ0uJvZcjN728zqzOy2QeZPMLPHE/NfNrM5WaytxsyeN7MtZrbZzP5qkDYXm9l+M3sz8bh9sOfKYI07zKw2se71g8w3M/tBYvttNLMlWaztlKTt8qaZHTCzLwxok/XtZ2YPmlmLmW1KmjbFzJ41s22Jn+VDLPupRJttZvapLNX2PTN7K/H3+6WZTR5i2WO+FjJc4x1mtivp73jlEMse8/89g/U9nlTbDjN7c4hls7IN08Y5NyYfxPuO3w6cCBQAG4DTB7T5C+CfE8PXAY9nsb4qYEliuBR4Z5D6Lgae9HAb7gAqjjH/SuAp4nc5Ow942cO/9W7i5+96uv2Ai4AlwKakad8FbksM3wbcNchyU4B3Ez/LE8PlWajtA0AgMXzXYLWl8lrIcI13AF9J4TVwzP/3TNU3YP49wO1ebsN0PcbynvsyoM45965zrhd4DLh2QJtrgYcTw08Al5pl535Xzrlm59zrieFOYCvxG4XnkmuBR1zcOmCymVV5UMelwHbn3Egvaksb59zviN9wJlny6+xh4EODLPpB4FnnXLtzrgN4Flie6dqcc8845yKJ0XXEb4PpmSG2XypS+X8ftWPVl8iOjwE/S/d6vTCWw30W0JA03sjR4dnfJvEC3w9MzUp1SRKHg84CXh5k9vlmtsHMnjKzBVktLH4X4mfM7LXEzckHSmUbZ8N1DP0P5eX26zPdOdecGN4NTB+kzVjYlp8m/klsMMO9FjLtlsShoweHOKw1Frbf+4E9zrltQ8z3ehsel7Ec7jnBzEqAfwe+4Jw7MGD268QPNZwB/CPwqyyX9z7n3BLgCuDzZnZRltc/LDMrAFYAPx9kttfb7ygu/vl8zJ1iZmZfByLAT4Zo4uVr4YfAScCZQDPxQx9j0fUce699zP8/JRvL4b4LqEkar05MG7SNmQWAMqAtK9XF1xkkHuw/cc79YuB859wB59zBxPBaIGhmFdmqzzm3K/GzBfgl8Y++yVLZxpl2BfC6c27PwBleb78ke/oOVyV+tgzSxrNtaWY3AlcDNyTefI6SwmshY5xze5xzUedcDPjREOv29LWYyI8/Bh4fqo2X23AkxnK4vwrMN7O5ib2764DVA9qsBvrOSvgI8JuhXtzpljg+96/AVufc3w/RZkbfdwBmtoz49s7Km4+ZTTSz0r5h4l+8bRrQbDXwycRZM+cB+5MOP2TLkHtLXm6/AZJfZ58C/mOQNk8DHzCz8sRhhw8kpmWUmS0H/hpY4ZwLDdEmlddCJmtM/h7nw0OsO5X/90y6DHjLOdc42Eyvt+GIeP2N7rEexM/meIf4t+hfT0y7k/gLGaCQ+Mf5OuAV4MQs1vY+4h/PNwJvJh5XAp8DPpdocwuwmfg3/+uAC7JY34mJ9W5I1NC3/ZLrM+C+xPatBZZm+e87kXhYlyVN83T7EX+jaQbCxI/7fob49zjPAduA/wKmJNouBR5IWvbTiddiHfBnWaqtjvix6r7XYN/ZYzOBtcd6LWRx+z2aeH1tJB7YVQNrTIwf9f+ejfoS0x/qe90ltfVkG6broStURUTy0Fg+LCMiIiOkcBcRyUMKdxGRPKRwFxHJQwp3EZE8pHAXEclDCncRkTykcBcRyUP/HwaEXbL2fweqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k = -0.6187760590682561, b = 0.6182123620096673\n"
     ]
    }
   ],
   "source": [
    "sol = lsqr(eq.a, eq.b)[0]\n",
    "v = sol[0:n]\n",
    "u = sol[n+1:]\n",
    "\n",
    "\n",
    "eq1 = LSQR_Helper(len(v), 2)\n",
    "for i in range(0, len(v)):\n",
    "    eq1.set_a(0, v[i])\n",
    "    eq1.set_a(1, 1)\n",
    "    eq1.set_b(u[i])\n",
    "    eq1.begin_new_row()\n",
    "    \n",
    "k,b = lsqr(eq1.a, eq1.b)[0]\n",
    "\n",
    "u_ = []\n",
    "for vi in v:\n",
    "    u_.append(k*vi+b)\n",
    "\n",
    "plt.plot(v, label='v')\n",
    "plt.plot(u, label='u')\n",
    "plt.plot(u_, label=\"u'\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "print(\"k = {}, b = {}\".format(k, b))"
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
No results found