lqr.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "from cvxpylayers.torch import CvxpyLayer\n",
    "from scipy.linalg import sqrtm\n",
    "from scipy.linalg import solve_discrete_are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "\n",
    "# generate problem data\n",
    "n, m = 4, 2\n",
    "noise = np.sqrt(0.25)\n",
    "Q0 = np.eye(n)\n",
    "R0 = np.eye(m)\n",
    "A = np.random.randn(n, n)\n",
    "A /= np.max(np.abs(np.linalg.eig(A)[0]))\n",
    "B = np.random.randn(n, m)\n",
    "W = noise**2 * np.eye(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8031165780081877\n"
     ]
    }
   ],
   "source": [
    "# compute lqr solution\n",
    "P = cp.Variable((n, n), PSD=True)\n",
    "R0cvxpy = cp.Parameter((m, m), PSD=True)\n",
    "\n",
    "objective = cp.trace(P@W)\n",
    "constraints = [cp.bmat([\n",
    "    [R0cvxpy + B.T@P@B, B.T@P@A],\n",
    "    [A.T@P@B, Q0+A.T@P@A-P]\n",
    "]) >> 0, P >> 0]\n",
    "R0cvxpy.value = R0\n",
    "result = cp.Problem(cp.Maximize(objective), constraints).solve()\n",
    "P_lqr = P.value\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.85426155, 0.12235637, 0.51550081, 1.05324043],\n",
       "       [0.12235637, 1.18209044, 0.09639788, 0.44752023],\n",
       "       [0.51550081, 0.09639788, 1.31473547, 0.6759572 ],\n",
       "       [1.05324043, 0.44752023, 0.6759572 , 2.86137885]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_lqr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up policy\n",
    "x = cp.Parameter((n, 1))\n",
    "P_sqrt = cp.Parameter((n, n))\n",
    "\n",
    "u = cp.Variable((m, 1))\n",
    "xnext = cp.Variable((n, 1))\n",
    "\n",
    "objective = cp.quad_form(u, R0) + cp.sum_squares(P_sqrt @ xnext)\n",
    "constraints = [xnext == A @ x + B @ u]\n",
    "prob = cp.Problem(cp.Minimize(objective), constraints)\n",
    "policy = CvxpyLayer(prob, [x, P_sqrt], [u])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize parameters\n",
    "Qt, Rt, At, Bt = map(torch.from_numpy, [Q0, R0, A, B])\n",
    "P_sqrt = torch.eye(n).double(); P_sqrt.requires_grad_(True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss(time_horizon, batch_size, P_sqrt, seed=None):\n",
    "    if seed is not None:\n",
    "        torch.manual_seed(seed)\n",
    "    x_batch = torch.randn(batch_size, n, 1).double()\n",
    "    P_sqrt_batch = P_sqrt.repeat(batch_size, 1, 1)\n",
    "    Qt_batch = Qt.repeat(batch_size, 1, 1)\n",
    "    Rt_batch = Rt.repeat(batch_size, 1, 1)\n",
    "    At_batch = At.repeat(batch_size, 1, 1)\n",
    "    Bt_batch = Bt.repeat(batch_size, 1, 1)\n",
    "    loss = 0.0\n",
    "    for _ in range(time_horizon):\n",
    "        u_batch, = policy(x_batch, P_sqrt_batch, solver_args={\"acceleration_lookback\": 0})\n",
    "        state_cost = torch.bmm(torch.bmm(Qt_batch, x_batch).transpose(2, 1), x_batch)\n",
    "        control_cost = torch.bmm(torch.bmm(Rt_batch, u_batch).transpose(2, 1), u_batch)\n",
    "        cost_batch = (state_cost.squeeze() + control_cost.squeeze())\n",
    "        loss += cost_batch.sum() / (time_horizon * batch_size)\n",
    "        x_batch = torch.bmm(At_batch, x_batch) + \\\n",
    "            torch.bmm(Bt_batch, u_batch) + \\\n",
    "            noise * torch.randn(batch_size, n, 1).double()\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "Kt = np.linalg.solve(R0 + B.T @ P_lqr @ B, -B.T @ P_lqr @ A)\n",
    "loss_lqr = loss(100, 6, torch.from_numpy(sqrtm(P_lqr)), seed=0).item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "it: 001, loss: 0.389, dist: 0.291, dist_P: 2.905\n",
      "it: 002, loss: 0.173, dist: 0.182, dist_P: 2.783\n",
      "it: 003, loss: 0.066, dist: 0.118, dist_P: 2.739\n",
      "it: 004, loss: 0.032, dist: 0.088, dist_P: 2.737\n",
      "it: 005, loss: 0.017, dist: 0.062, dist_P: 2.742\n",
      "it: 006, loss: 0.010, dist: 0.048, dist_P: 2.736\n",
      "it: 007, loss: 0.007, dist: 0.040, dist_P: 2.743\n",
      "it: 008, loss: 0.006, dist: 0.036, dist_P: 2.751\n",
      "it: 009, loss: 0.007, dist: 0.041, dist_P: 2.776\n",
      "it: 010, loss: 0.006, dist: 0.038, dist_P: 2.760\n",
      "it: 011, loss: 0.007, dist: 0.037, dist_P: 2.751\n",
      "it: 012, loss: 0.005, dist: 0.029, dist_P: 2.710\n",
      "it: 013, loss: 0.005, dist: 0.027, dist_P: 2.702\n",
      "it: 014, loss: 0.005, dist: 0.025, dist_P: 2.715\n",
      "it: 015, loss: 0.005, dist: 0.022, dist_P: 2.682\n",
      "it: 016, loss: 0.004, dist: 0.023, dist_P: 2.680\n",
      "it: 017, loss: 0.003, dist: 0.026, dist_P: 2.686\n",
      "it: 018, loss: 0.003, dist: 0.034, dist_P: 2.663\n",
      "it: 019, loss: 0.003, dist: 0.020, dist_P: 2.688\n",
      "it: 020, loss: 0.002, dist: 0.028, dist_P: 2.650\n",
      "it: 021, loss: 0.005, dist: 0.030, dist_P: 2.705\n",
      "it: 022, loss: 0.004, dist: 0.029, dist_P: 2.689\n",
      "it: 023, loss: 0.007, dist: 0.043, dist_P: 2.713\n",
      "it: 024, loss: 0.001, dist: 0.019, dist_P: 2.658\n",
      "it: 025, loss: 0.001, dist: 0.034, dist_P: 2.633\n",
      "it: 026, loss: 0.001, dist: 0.019, dist_P: 2.647\n",
      "it: 027, loss: 0.002, dist: 0.016, dist_P: 2.669\n",
      "it: 028, loss: 0.002, dist: 0.016, dist_P: 2.666\n",
      "it: 029, loss: 0.002, dist: 0.017, dist_P: 2.673\n",
      "it: 030, loss: 0.002, dist: 0.019, dist_P: 2.681\n",
      "it: 031, loss: 0.002, dist: 0.017, dist_P: 2.675\n",
      "it: 032, loss: 0.002, dist: 0.017, dist_P: 2.671\n",
      "it: 033, loss: 0.002, dist: 0.018, dist_P: 2.675\n",
      "it: 034, loss: 0.002, dist: 0.016, dist_P: 2.669\n",
      "it: 035, loss: 0.002, dist: 0.016, dist_P: 2.669\n",
      "it: 036, loss: 0.001, dist: 0.015, dist_P: 2.661\n",
      "it: 037, loss: 0.001, dist: 0.015, dist_P: 2.662\n",
      "it: 038, loss: 0.001, dist: 0.015, dist_P: 2.657\n",
      "it: 039, loss: 0.001, dist: 0.015, dist_P: 2.653\n",
      "it: 040, loss: 0.001, dist: 0.016, dist_P: 2.652\n",
      "it: 041, loss: 0.001, dist: 0.016, dist_P: 2.650\n",
      "it: 042, loss: 0.001, dist: 0.018, dist_P: 2.646\n",
      "it: 043, loss: 0.001, dist: 0.016, dist_P: 2.645\n",
      "it: 044, loss: 0.001, dist: 0.015, dist_P: 2.652\n",
      "it: 045, loss: 0.001, dist: 0.017, dist_P: 2.646\n",
      "it: 046, loss: 0.001, dist: 0.016, dist_P: 2.653\n",
      "it: 047, loss: 0.001, dist: 0.017, dist_P: 2.648\n",
      "it: 048, loss: 0.001, dist: 0.018, dist_P: 2.651\n",
      "it: 049, loss: 0.001, dist: 0.017, dist_P: 2.661\n",
      "it: 050, loss: 0.001, dist: 0.017, dist_P: 2.670\n"
     ]
    }
   ],
   "source": [
    "losses = []\n",
    "opt = torch.optim.SGD([P_sqrt], lr=.5)\n",
    "test_loss = None\n",
    "for k in range(50):\n",
    "    with torch.no_grad():\n",
    "        test_loss = loss(100, 6, P_sqrt.detach(), seed=0).item()\n",
    "        K_np = (torch.solve(-Bt.t() @ P_sqrt.t() @ P_sqrt @ At, Rt + Bt.t() @ P_sqrt.t() @ P_sqrt @ Bt).solution).detach().numpy()\n",
    "        dist = np.linalg.norm(K_np - Kt)\n",
    "        P = (P_sqrt.t() @ P_sqrt).detach().numpy()\n",
    "        dist_P = np.linalg.norm(P_lqr - P)\n",
    "        losses.append(test_loss)\n",
    "        print(\"it: %03d, loss: %3.3f, dist: %3.3f, dist_P: %3.3f\" % (k+1, test_loss - loss_lqr, dist, dist_P))\n",
    "    opt.zero_grad()\n",
    "    l = loss(100, 6, P_sqrt, seed=k+1)\n",
    "    l.backward()\n",
    "    opt.step()\n",
    "    if k == 25:\n",
    "        opt = torch.optim.SGD([P_sqrt], lr=.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAACnCAYAAAAPF9kJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXxUlEQVR4nO3db0wb9/0H8PcXm0D4axva/IPMdjptWacto0Tr+qhSiKZVyh60FNpNkZj0A9ZuDyathdLH1SpQNkV7sJ9Cqj1Zpy2Jf6qmrtpWvHR70Goi4GRtpKpa4qyBrckI4GgEHIz5/h74e5cz2IDPd9wZv18SSnx3Pr6c7Y+/d/f9fD9CSgkiokJVON0AIipNDB5EZAqDBxGZwuBBRKYweBCRKQweRGSK1+kGFKu5uVkGg0Gnm0G0I01OTt6RUj6Ua13JB49gMIiJiQmnm0G0IwkhPs23jqctRGQKgwcRmWJr8BBCdKqfYSFEZyHbqMc+9dNXTDu6urrw05/+tJhdENEatl3zEEK0AUhIKaMAIkKIeSFEVEqZ2OI2bQAmAcQA9BbTlomJCezatauYXRDRGnZeMA0DOAogqh7PqWWxLW5zRkoZsaIhfr8fc3NzVuyKSkAqlcL09DSSyaTTTSkJ1dXVaGlpQWVlZUHPsy14qA9+BACEED61LFbANmEhRAcAn2FbnRDiBIATjzzyyKZtCQQCmJ+fL+bPoRIyPT2N+vp6BINBCCGcbo6rSSkxOzuL6elphEKhgp67XRdMhwEcL2QbKeWIlDKqgsaQFlwM69+WUvY1NjZu+svZ8ygvyWQSTU1NDBxbIIRAU1OTqV6a7cFDXQQdllLGt7qNuoA6YNhEO50xhT2P8sPAsXVmj5Wtg8TUaUdMShlXPYfA2iCSaxsAcfWjCaw95SmE1vOQUvJNRWQRu++2nAGQUB/YsJTSr9ZNAjiGTG8i3zadQoiw2mawmLYEAgGk02ksLCygvr6+mF0RkWLnBdMYgEN51j2m/rvRNpbcaQEyPQ8AmJ+fZ/Agx0WjUcRiMbS1tQHIfLmdO3cOw8PD+vpEIoFwOIy5uTkkEgl0dnZuaR/9/f0YHBxEPB7H0NAQAODSpUs4fvw4Ojo6LP07Sj63ZSsCgQAAYG5uDgcPHnS4NVTORkdHcf36dT1QAEA8Hkc0GtXXA0Bf34NxkZFIBIODg/pzNtrH8PAwuru7MTY2pgeczs5O+P1+3LhxAz5f1n2HopRF8DD2PIic1N/fj7WTjofDYXR3dwMABgcH171POzs7MTg4iKGhIfh8vk33kUsgEEA8Htd7KlYoi+Bh7HlQ+fnRj36EK1eu2Po7jhw5gtOnT2+4TTQaRTic+6bhwMDAhuvD4TCi0Sh8Pt+G+8glHs/ce7AycABlEjzY86BSkUgkci4v9ItPO41JJBKYm5vD9evXrWheFrtv1WpXeY4CuJTrImi+bdTyBDIjTOPF3Kplz6O8bdYj2C4dHR16L2CtaDSK9vb2vO9R7ZQjHA5vuA/tomg4HLb8Aulatg0SMyS9RaSUgwDOrh0lmm8btd1x4wjTYtpSW1sLr9fLngc57sKFC+jv789alkgkEAgE4PP5MDw8jJGRkaz1o6Oj6Ovr009XNtrHdnJrYlw7Mr0OnRCizWzvQwiBQCDAngc5rrOzEz6fDyMjIwiHw/rdD62X0NfXh0gkgkgkAp/Pp5/GGO+sbLSPeDyOc+fOIR6PZ93KtYPYjnKTqicxKaXMOaZj7TZqaHpCSjmq1p0BMGY87TEkxvX+4x//2LQNX/ziF/GVr3wF58+fL/bPIZf7+OOPcfjwYaebUVLyHTMhxKSUsj3Xc1ybGLeZQhLjALDnQWQxVybG4cGFUo2W72Ka3+/nNQ8iC9k9DWFW0pvKVdnKNueRPWzdV8zdFoCZtURWc3Ni3AUVWIDMKU1ROKcHkbXcnBgXzbXcrEAggLt37yKdTsPj8Vi5a6KyVDalF7RRpvlG8BFRYcouePC6B5E1yiZ4cIg6kbUcL/qktmtbW9jJyqJPAHseRFZztOiT2q4DwLPIFHgysqzoE8CeB5HV7Ox5hJE9YjTnDOgquIzleP4ZKeUhKeWzawOOGex5EFnL0aJPm9iw6FOhtODBngc5JR6PZ80vunZeUiAz5SCAnPOXas9PJBJ6Vm08HkdHR4etCXD5bNdkQAXlrQCZok/a/4UQk7lOeQpRVVWFmpoa9jzIMdpUgcb5RY1GR0eRSCSyZgSLRCIYGRnBwMBAzucnEgmEQiFH3te2B4+t5LbkeU7YEEDWpfMXUm5Sw+S48vXkk0+uW9bV1YUXX3wRi4uLeOqpp9at7+npQU9PD+7cuZPzw/7CCy+gu7sbU1NTaG1tLbqN+eYvPXToEPr6+nJOXqz1TpzgeG5LHnE8mOMDyFH0qdCsWoDJceRe0Wg072Q+2vylmrm5OcRiMYyOjmJwcBCTk2vvNWwPR3NbpJQJFWC6AfiEEHE1e1jMyqJPGibHla+//OUvedfV1NRsuL65uXnD9Vb0OgBsWBbB2GMOBAJoa2vTr3NYWU6hEE7ntmh3W9blsVhZ9Enj9/tx7do1q3dLVBRt/tJ8px8TExM4c+ZM3ucPDw9vuN4uZTPCFGDPg9wpFovB5/NhcHBQPz2JxWJIJBKIRqPo6urKW25Bq8cCbH/eVlkFD6blk5O0+UUnJiYQiUQwOjqK48ePY3Z2FkBm/tJAIIBoNIq5uTmMjo4iFovpvYpEIoGxsTFMTEzoAUObzzQajeadVd0uZVG3RRMIBLC0tIRkMonq6mqnm0NlJhwO48KFC1nLjGUlgezCTOFwGJFIRJ/MuLOzM+fpydp9bpey63kAHGVKpSEcDmNsbMyWItVWKLueB5AJHvv27XO4NUSbGxvLlbnhDo5XjFPbtQFo10otGJ5rScU4DYeoE1nHlVm1hopx/erxBbVNUYw9D9rZpJRQ44toE2ZrN7k1q7YLOSrGFdsg9jzKQ3V1NWZnZ01/KMqJlBKzs7OmbiC4NavWB8BY1jtXqcqCsedRHlpaWjA9PY2ZmRmnm1ISqqur0dLSUvDzXJtVuxkziXGNjY0QQrDnscNVVlYiFAo53YwdzxUV43LYtGKcmcS4iooK+Hw+9jyILODWrFrLK8ZpmFlLZA23ZtUmrK4Yp+GcHkTWKCh4CCGOSCmvqP8fAyCllBdzbWtBVq2lFeM07HkQWaPQ0xb9tENK+WfkuPXqdkyOI7LGlnoeQoheAP3InFoMqcVzAJzJyCkC0/KJrLGl4CGlPCuEOI/MdYvLNrfJVlrPgyMQiYqz5dMWKeVdY+AQQgTtaJDdAoEA0uk0FhYWnG4KUUkr6JqHEOJ1IcTXhBD/C+D7Qoj/saldtuEQdSJrFHrBdFT1Po5KKV8BcMOGNtmKQ9SJrFHoOI82IcTXAJxTjzcc3rmVlPx8qfdCiGEAr6vNuozp+sVgz4PIGoUGjzgyH+QhdQcmb/DYSkr+Jqn3lha61rDnQWSNgoKHlPKyEALqWsclbcBYHmFkehzaYK9cmbE5U+9V7+OMXeUXAPY8iIpV6AXTZ5D5wAtkLpg+nW9bKWVESjmonpcvJT9f6j2gCl2r4k/rav0JIU4IIUbv3r1byJ/AngeRRQq9YJqQUg5JKc9KKb+PTBDZClOFrlWeSwTAkBaADOsLzqoFMtXBKisrGTyIilRo8Jhd83jTT+AmKfk5U+9Vb2PAsDznLGRmCCEQCAT0WhlEZE6hF0yPCiHakblwGkbmg58zMQ5Yn5KPTMFqYxA5j+yMWZ+qUwtkz9+xrtB1MVpaWjA1NWXV7ojKUqEXTM8arntMSClP5du2gJT8dan3dhW61gSDQVy9etXKXRKVnUJT8r8GoENK+X0hRKMxRX+tAlPyc21j+Z0WTSgUwjvvvMP8FqIimJlJ7BUgk+uCzDWKkhMKhZBMJnHr1i2nm0JUsgoNHn4VNDS+vFu6WDAYBADcuFFyo+uJXKPQ4HFIJcc9LYR4HSU4GRAAfWbtf/7zn842hKiEmblgegyZMRuXpJT/Z0+z7MWeB1HxCp4AWU0/+Gcb2rJtdu/ejT179jB4EBXB9rotbhUKhXjaQlQEWyvGFZmSn3O5VYLBIMbHx63cJVFZsa3nYUjJ1xLkzq7NTzGk5Os5LBstt1IoFMLNmzeRTqet3jVRWbDztCWM7GS4XPkpOVPyN1humVAohJWVFfzrX/+ycrdEZcO24FFkSv5GqfpQ+zSVkq/hHRei4mzXBdOCU/I3YzYlX6ON9WDwIDLH9uBhJiV/g+WWOXjwIIQQvONCZJKtwWNtSr7KkjU6j+zkOZ86tcm33DK7du3CgQMH2PMgMsm2W7VFpuTnXG61UCjE4EFkkm3Bw4KU/JzLrRQKhfDee+/Z/WuIdqSyHWEKZO64TE9PY3l52emmEJWcsg4eoVAIUkrcvHnT6aYQlZyyDx4AU/OJzCjr4MGBYkTm2X2rtk0I0bfB+k7102e4swIhxLC6tevb6PnFamlpgdfrZfAgMsHOW7UdyNSdncyzvg3ZdWrH8KA0pS11atfyeDw4ePAgT1uITLDzVm1Uq9WSZ5MOZOevQAgRViNRbalTm0swGGTPg8gEJ695JAA0GR4HsMU6tVbiQDEic2ydDGgT5wGcNTzWc1mklCPa/4UQk0KIqJQyK0XfKsFgELdv38bS0hJ2795tx68g2pEc63moYDCoehhhZBLfJrZap7bYlHwNb9cSmeP0rdqwGoY+h8xUgwlkgohxaHrOOrXFpuRrmJpPZI7dd1u6AfiEEHEtV8WYFAegzZA01w/YX6d2LW2sB3seRIWx9W4LsnsQ2nJjUtzI2vVq+bbcaQGAvXv3oqqqij0PogI5fdriuIqKCnz1q1/Fu+++Cyml080hKhllHzwA4Hvf+x4+/PBDXLp0yemmEJUMBg8A3/nOd1BTU4M33njD6aYQlQwGDwANDQ3o7u7Gb37zGywsLDjdHKKSwOCh9Pb2YmFhAb/97W+dbgpRSWDwUB5//HF86UtfwtmzZzffmIgYPDRCCPT29mJ8fBwffvih080hcj23zufRaUiMs7TM5EZOnjyJXbt2sfdBtAV2FrruANC/wXptPo+IlHIUaiTpdhS5zqepqQnPPPMM3nzzTSwtLW3XryUqSXbWqo0CGNtgk5zzeWAbilxvpLe3F4lEApHItg1yJSpJbpzPY9Mi13Z68skncfjwYQwMDGBqamq7fi1RyXEyeJxHdlDw5dswF6tS8nPsFxcuXMDi4iJOnDjBcR9EebhuPg9ssci1VSn5uTz66KM4f/48PvroIzz//PNIp9OW/w6iUuf0rdpc83nYXuR6K775zW/i5z//OX7/+99jYGBg8ycQlRk3zuexLUWut+IHP/gBPvnkE/zsZz9DVVUVhoaGUF9f71RziFxFlHoaent7u5yYmLBt/ysrK+jp6cGvf/1rBAIBvPTSS/jhD3/IIEJlQQgxKaVsz7XO6dMW1/N6vXjzzTcxPj6Ob3zjG3j11VcRCoXw8ssv4+LFiyySTWWLPY8CjY+P47XXXsMf//hHpFIp1NXV4dixY/jyl7+M+vp61NXVoa6uDslkErdv38Z//vMf3L59G+l0GnV1daitrUVtbS08Hg/S6TRWVlaQTqdRW1uLvXv3Ys+ePdi7dy8CgQDq6ur0fXq9XqRSKaysrCCVSkEIAa/XC6/XC4/Hg6qqKqhTwG2RTqcRj8dx9epVXL16FX6/H9/97nfh9/u3rQ1kv416HgweJi0sLODixYv4wx/+gD/96U/49NNPsbq6um67QCCAhx9+GF6vFwsLC7h37x7u3buHdDqtf/A9Hg8WFhaQSqVMt6eqqgr79u3Dvn37sH//fjQ3N6O+vh4NDQ1oaGjA7t279WDj9WYudaXTaayuriKdTqOyshKNjY1oaGhAY2MjZmZmcPnyZVy5cgWXL1/G7Oys3laPx4OZmRkkk8msNlRXV6Orqwt9fX144oknthzMVlZWsLy8nLX/1dVV3L9/X/9JJpNYWlrC4uIilpaWkE6nsWvXLlRWVqKysjIrkHo8HgDA8vIyUqmU3jusrq7Wf7xer94+IYQegKuqqvTn0w4PHvX19fKxxx7LWtbV1YUXX3wRi4uLeOqpp9Y9p6enBz09Pbhz5w46O9fXlHrhhRfQ3d2NqakpnDx5ct36H//4xzhx4gQ++eQT9Pc/GIGvfRD7+/tx7Ngx/Pvf/8bLL7+87vk/+clP8MQTT+CDDz7Aq6++qi/XPkS9vb1oamrC3/72N7z11ltIp9OQUqKiogJCCDz33HPYs2cP/v73v+O9996DlFJ/7v3799HU1IT//ve/mJ+fLyogAZmgVFdXh6qqKkgpIaXEt771Lfj9fkxPT+PatWuoqanB0tISPvvsM72XVVtbq7dLazeQKfG5f/9+LC8v49atW7h//77rpn/UAtHy8rLeNiEEhBDYvXs3fD4fKisr8dlnn+nLtb/P7/fjoYceAgBcu3ZNf6727/79+3Hw4EGsrKwgFsu+ibi6uoq9e/fC5/MhmUzixo0bWfuXUuLAgQNobGzE4uLiuvUA8Mgjj6C5uRnz8/P4+OOPsbq6mvVTW1ur93bXfm4A4PTp0zhy5Aii0Shee+01/PWvf80bPJws+rTjVFRUoKKiAvv378eBAwcwMzNT0PO1N+3Ro0dx5MgRNDQ0YHx8fN12J0+exBe+8AW8/fbb+hvU6Fe/+hVaW1tx7tw5/OIXv9B7GFJKrK6u4vTp0/D7/fjd736Ht956CwD0D3o6ncYrr7yChx9+GJOTk3jnnXfW7f+Xv/wlAODUqVO4ffs2AKCurg6f//zn8eijj+K5557DRx99hPfffx9TU1P67wYywaitrQ2VlZW4evUq7ty5A4/Ho7/56+vr8fTTT0MIgXfffRe3bt3Sj2tFRQVaW1vx0ksvwePx4NSpU5iamtL/LgBobW3F888/DwB44403MD8/j4qKCj3whUIhfPvb30YqlcLZs2dx7949/e9fXV1FMBjE448/jpWVFUQiEaysrOjrpZRobW3F4cOHkUql9HlvtR8gEzwOHDiAVCqlz8ivrZNSIpVKIZFIIJlM6vvWVFRUwOv1orGxEVVVVXq7teOnBYmamhp9n9o67XEikcDq6iqWlpb0Lxytl+XxePD1r38doVAI8/PzuHnz5ibvyI2VfM/DqdMWonLAuy1EZDkGDyIyhcGDiExh8CAiU0r+gqkQYgbAp1vYtBnAHZubUwy2rzhsX3Hyte9zUsqHcj2h5IPHVgkhJvJdNXYDtq84bF9xzLSPpy1EZAqDBxGZUk7BY9TpBmyC7SsO21ecgttXNtc8iMha5dTzICILMTHOAaoOTbsqdqUt68SDyZ/jTszbStbSXlPDFJw76jXe8cHDbS+Ympv1WQCThmValbx+9fiC2sYR6pgBwFEAl1TlPtccS3UMtXo+AS0Iu6V9qi0+ZObwPWN47IrXWAgxDOB19bDL7PHb0cHDTS+YRkoZVe0KGBbnrJLnxJtf9Yq0b8uIEGJeCBFVqx0/lurYDUopjwOICSEkgFEXvtbtAC4ZHrvmNQbQhsyXVwxAr2pLwcdvp1/zcLR0ZQEcrZK3RhjA8RxtccWxlFImVODQfr926ueK9qnf26Gdqhi46TU+I6U8JKV8VlUxAEwcv50ePNz0gpUEVXjcWHQc6tvRVcdSvbG7tW9KuKR96pitK1LmMmFVbK3TcIpa8PHb0actJWRLVfIcMIzsXohrSCljQoi4EOK6lPLQ5s/YNu2AXrT9KIAmIUQcLnqNpZQj2v+FEJOG09KC7PSeh2tesE24okqekfpGGpZSasfLFcdSCNGmFQTTutzqsSvaJ6WMaj/q94+pY+iK11j1NowlELUeRsHHb6f3PM4ju+KcGz6U6yrpualKHqC3MSaljBsu7rrlWLYj84Y30uocu6F9APTTqg5kThHi6li64TWOIzsoBLReHAo8fjt+hKnhxQKQ+WZwqi2lQL3pL+DBxbOwlNKv1rniWAoh+vDgGzNuuJXsiva5neE6RxiZLwltHEpBx2/HBw8issdOv+ZBRDZh8CAiUxg8iMgUBg8iMoXBg4hMYfCgbSGE8KnRjB2bb02lgMGDtoUaDTrHsRc7B4MHEZnC4EHbQhvyrk5f2tR8EVTCGDxouxwHMAuVSwFg0OH2UJEYPGi7dCCTkKVlxLoxu5kKwOBB22VOJbAdV6cunJSpxDF4kO3U9Q7tGscY1IQ5VNqYVUtEprDnQUSmMHgQkSkMHkRkCoMHEZnC4EFEpjB4EJEpDB5EZAqDBxGZ8v/YeOtE2jGGpgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x177.994 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from utils import latexify\n",
    "import matplotlib\n",
    "\n",
    "latexify(fig_width=4)\n",
    "fig, ax = plt.subplots()\n",
    "plt.xticks([0, 10, 20, 30, 40, 50])\n",
    "plt.semilogy(losses, color='k', label='COCP')\n",
    "plt.gca().yaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())\n",
    "plt.axhline(loss_lqr, linestyle='--', color='k', label='LQR')\n",
    "plt.legend()\n",
    "plt.ylabel(\"cost\")\n",
    "plt.xlabel(\"$k$\")\n",
    "plt.subplots_adjust(left=.15, bottom=.2)\n",
    "plt.savefig(\"lqr.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}