limit_cycles.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook provides examples to go along with the [textbook](https://underactuated.csail.mit.edu/limit_cycles.html).  I recommend having both windows open, side-by-side!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import mpld3\n",
    "import numpy as np\n",
    "from IPython.display import display\n",
    "from pydrake.all import DirectCollocation, PiecewisePolynomial, Solve, eq\n",
    "from pydrake.examples import VanDerPolOscillator\n",
    "\n",
    "from underactuated import running_as_notebook\n",
    "\n",
    "if running_as_notebook:\n",
    "    mpld3.enable_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finding the limit cycle of the Van der Pol oscillator\n",
    "\n",
    "by setting up a simple trajectory optimization, where the timestep, $h$, is a decision variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO(russt): visualize this in meshcat so we can watch the convergence again.\n",
    "\n",
    "\n",
    "def vdp_limit_cycle():\n",
    "    plant = VanDerPolOscillator()\n",
    "    context = plant.CreateDefaultContext()\n",
    "\n",
    "    dircol = DirectCollocation(\n",
    "        plant,\n",
    "        context,\n",
    "        num_time_samples=61,\n",
    "        minimum_timestep=0.01,\n",
    "        maximum_timestep=0.5,\n",
    "    )\n",
    "    prog = dircol.prog()\n",
    "\n",
    "    # Constrain all timesteps, $h[k]$, to be equal, so the trajectory breaks are evenly distributed.\n",
    "    dircol.AddEqualTimeIntervalsConstraints()\n",
    "\n",
    "    # Initial state on the surface of section (and velocity > .1).\n",
    "    prog.AddBoundingBoxConstraint(\n",
    "        [0.0, 0.1], [0.0, 10.0], dircol.initial_state()\n",
    "    )\n",
    "\n",
    "    # Periodicity constraint.\n",
    "    prog.AddLinearConstraint(eq(dircol.final_state(), dircol.initial_state()))\n",
    "\n",
    "    # Help the solver with an initial guess (circular trajectory).\n",
    "    samples = np.linspace(0, 2 * np.pi, 10)\n",
    "    x_guess = np.vstack(\n",
    "        ([2 * np.sin(t) for t in samples], [2 * np.cos(t) for t in samples])\n",
    "    )\n",
    "    initial_x_trajectory = PiecewisePolynomial.FirstOrderHold(samples, x_guess)\n",
    "    dircol.SetInitialTrajectory(PiecewisePolynomial(), initial_x_trajectory)\n",
    "\n",
    "    fig = plt.figure(figsize=(6, 6))\n",
    "    (h,) = plt.plot([], [], \".-\")\n",
    "    plt.xlim((-2.5, 2.5))\n",
    "    plt.ylim((-3.0, 3.0))\n",
    "    plt.axis(\"equal\")\n",
    "\n",
    "    def draw_trajectory(t, x):\n",
    "        h.set_xdata(x[0, :])\n",
    "        h.set_ydata(x[1, :])\n",
    "        fig.canvas.draw()\n",
    "        if plt.get_backend() == \"MacOSX\":\n",
    "            plt.pause(1e-10)\n",
    "\n",
    "    if False:\n",
    "        dircol.AddStateTrajectoryCallback(draw_trajectory)\n",
    "\n",
    "    result = Solve(prog)\n",
    "    assert result.is_success()\n",
    "\n",
    "    x_trajectory = dircol.ReconstructStateTrajectory(result)\n",
    "\n",
    "    x_knots = np.hstack(\n",
    "        [\n",
    "            x_trajectory.value(t)\n",
    "            for t in np.linspace(\n",
    "                x_trajectory.start_time(), x_trajectory.end_time(), 100\n",
    "            )\n",
    "        ]\n",
    "    )\n",
    "    plt.plot(x_knots[0, :], x_knots[1, :], \"b-\", marker=\".\")\n",
    "    display(mpld3.display())\n",
    "\n",
    "\n",
    "vdp_limit_cycle()"
   ]
  },
  {
   "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.10.9 (main, Dec 15 2022, 10:44:50) [Clang 14.0.0 (clang-1400.0.29.202)]"
  },
  "vscode": {
   "interpreter": {
    "hash": "b0fa6594d8f4cbf19f97940f81e996739fb7646882a419484c72d19e05852a7e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}