Birkhoff Docs (#1088)

* Added documentation for birkhoff method * Simplified initialization of Birkhoff transcription and updated examples to use new API * added an example using Birkhoff transcription to the docs * fixed references on moon landing problem * added moon landing problem to toc tree * Fixed a bug with min time climb where number of nodes was specified incorrectly and updated new load case test --------- Co-authored-by: Kaushik Ponnapalli <[email protected]>
OpenMDAO · Aug 6, 2024 · ea5faad · ea5faad
1 parent a7b7750
commit ea5faad
Show file tree

Hide file tree

Showing 34 changed files with 421 additions and 64 deletions.
diff --git a/docs/dymos_book/_toc.yml b/docs/dymos_book/_toc.yml
@@ -30,6 +30,7 @@ parts:
     - file: examples/length_constrained_brachistochrone/length_constrained_brachistochrone
     - file: examples/min_time_climb/min_time_climb
     - file: examples/mountain_car/mountain_car
+    - file: examples/moon_landing/moon_landing
     - file: examples/multi_phase_cannonball/multi_phase_cannonball
     - file: examples/multibranch_trajectory/multibranch_trajectory
     - file: examples/racecar/racecar

diff --git a/docs/dymos_book/api/transcriptions_api.ipynb b/docs/dymos_book/api/transcriptions_api.ipynb
@@ -102,6 +102,23 @@
     "om.show_options_table(tx)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Birkhoff Options"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tx = dm.transcriptions.Birkhoff(num_nodes=5, grid_type='cgl')\n",
+    "om.show_options_table(tx)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -122,6 +139,13 @@
     "tx = dm.transcriptions.ExplicitShooting(num_segments=1, order=5)\n",
     "om.show_options_table(tx)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {

diff --git a/docs/dymos_book/examples/moon_landing/moon_landing.ipynb b/docs/dymos_book/examples/moon_landing/moon_landing.ipynb
@@ -0,0 +1,287 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8f7334f0",
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# This cell is mandatory in all Dymos documentation notebooks.\n",
+    "missing_packages = []\n",
+    "try:\n",
+    "    import openmdao.api as om\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !python -m pip install openmdao[notebooks]\n",
+    "    else:\n",
+    "        missing_packages.append('openmdao')\n",
+    "try:\n",
+    "    import dymos as dm\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !python -m pip install dymos\n",
+    "    else:\n",
+    "        missing_packages.append('dymos')\n",
+    "try:\n",
+    "    import pyoptsparse\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !pip install -q condacolab\n",
+    "        import condacolab\n",
+    "        condacolab.install_miniconda()\n",
+    "        !conda install -c conda-forge pyoptsparse\n",
+    "    else:\n",
+    "        missing_packages.append('pyoptsparse')\n",
+    "if missing_packages:\n",
+    "    raise EnvironmentError('This notebook requires the following packages '\n",
+    "                           'please install them and restart this notebook\\'s runtime: {\",\".join(missing_packages)}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c9df3b73",
+   "metadata": {},
+   "source": [
+    "# Moon Landing Problem\n",
+    "\n",
+    "The Moon landing problem is a version of the soft landing problem presented in {cite}`Meditch1964`. The problem is simplified to have one degree-of-freedom and normalized such that the Moon's gravity is unity. The goal is to minimize the amount of fuel consumed or, stated differently, maximize the final mass, while bringing the lander down to the surface for a soft landing."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a8b1357a",
+   "metadata": {},
+   "source": [
+    "## State and control variables\n",
+    "\n",
+    "This system has three state variables, the altitude ($h$), velocity ($v$), and mass ($m$) of the lander.\n",
+    "\n",
+    "This system has one control variable, ($T$), the thrust applied to the vehicle.\n",
+    "\n",
+    "The dynamics of the system are given by\n",
+    "\n",
+    "\\begin{align}\n",
+    "  \\dot{h} &= v \\\\\n",
+    "  \\dot{v} &= -1 + \\frac{T}{m} \\\\\n",
+    "  \\dot{m} &= -\\frac{T}{2.349}\n",
+    "\\end{align}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "71a99c51",
+   "metadata": {},
+   "source": [
+    "## Problem Definition\n",
+    "\n",
+    "We seek to maximize the final mass of the vehicle while bringing it to a soft landing.\n",
+    "\n",
+    "\\begin{align}\n",
+    "  \\mathrm{Minimize} \\, J &= m_f\n",
+    "\\end{align}\n",
+    "\n",
+    "The initial conditions are\n",
+    "\\begin{align}\n",
+    "  h_0 &= 1 \\\\\n",
+    "  v_0 &= -0.783 \\\\\n",
+    "  m_0 &= 1\n",
+    "\\end{align}\n",
+    "and the terminal constraints are\n",
+    "\\begin{align}\n",
+    "  h_f &= 0 \\\\\n",
+    "  v_f &= 0\n",
+    "\\end{align}\n",
+    "\n",
+    "Additionally, the thrust is constrained to be positive but remain under 1.227.\n",
+    "\n",
+    "\\begin{align}\n",
+    "  0 \\le T \\le 1.227 \n",
+    "\\end{align}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "48fa15f1",
+   "metadata": {},
+   "source": [
+    "## Defining the ODE\n",
+    "\n",
+    "The following implements the dynamics of the Moon landing problem described above."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d1fa5d91",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import openmdao.api as om\n",
+    "\n",
+    "\n",
+    "class MoonLandingProblemODE(om.ExplicitComponent):\n",
+    "    def initialize(self):\n",
+    "        self.options.declare('num_nodes', types=int)\n",
+    "\n",
+    "    def setup(self):\n",
+    "        nn = self.options['num_nodes']\n",
+    "\n",
+    "        # inputs\n",
+    "        self.add_input('h', val=np.ones(nn), units=None, desc='Altitude')\n",
+    "        self.add_input('v', val=np.ones(nn), units='1/s', desc='Velocity')\n",
+    "        self.add_input('m', val=np.ones(nn), units=None, desc='Mass')\n",
+    "        self.add_input('T', val=np.ones(nn), units=None, desc='Thrust')\n",
+    "\n",
+    "        # outputs\n",
+    "        self.add_output('h_dot', val=np.ones(nn), units='1/s', desc='Rate of change of Altitude')\n",
+    "        self.add_output('v_dot', val=np.ones(nn), units='1/s**2', desc='Rate of change of Velocity')\n",
+    "        self.add_output('m_dot', val=np.ones(nn), units='1/s', desc='Rate of change of Mass')\n",
+    "\n",
+    "        # partials\n",
+    "        ar = np.arange(nn)\n",
+    "        self.declare_partials(of='h_dot', wrt='v', rows=ar, cols=ar, val=1.0)\n",
+    "        self.declare_partials(of='v_dot', wrt='m', rows=ar, cols=ar)\n",
+    "        self.declare_partials(of='v_dot', wrt='T', rows=ar, cols=ar)\n",
+    "        self.declare_partials(of='m_dot', wrt='T', rows=ar, cols=ar, val=-1/2.349)\n",
+    "        self.declare_partials(of='m_dot', wrt='T', rows=ar, cols=ar, val=-1/2.349)\n",
+    "\n",
+    "    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):\n",
+    "        v = inputs['v']\n",
+    "        m = inputs['m']\n",
+    "        T = inputs['T']\n",
+    "\n",
+    "        outputs['h_dot'] = v\n",
+    "        outputs['v_dot'] = -1 + T/m\n",
+    "        outputs['m_dot'] = -T/2.349\n",
+    "\n",
+    "    def compute_partials(self, inputs, partials, discrete_inputs=None):\n",
+    "        m = inputs['m']\n",
+    "        T = inputs['T']\n",
+    "\n",
+    "        partials['v_dot', 'T'] = 1/m\n",
+    "        partials['v_dot', 'm'] = -T/m**2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "18e69427",
+   "metadata": {},
+   "source": [
+    "## Solving the Moon landing problem with Dymos\n",
+    "\n",
+    "The optimal solution to this problem is known to have _bang-bang_ control. That is, the control has a \"jump\" that render it discontinuous in time. Capturing this behavior accurately requires the use of grid refinement for the Gauss-Lobatto and Radau pseudospectral transcriptions but the Birkhoff pseudospectral transcription can be used to handle this behavior without the use of any grid refinement. The following code shows the use of the Birkhoff pseudospectral transcription to solve the problem."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "65ab9d17",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import dymos as dm\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "p = om.Problem(model=om.Group())\n",
+    "p.driver = om.pyOptSparseDriver()\n",
+    "p.driver.declare_coloring()\n",
+    "p.driver.options['optimizer'] = 'IPOPT'\n",
+    "p.driver.opt_settings['hessian_approximation'] = 'limited-memory'\n",
+    "p.driver.opt_settings['print_level'] = 0\n",
+    "p.driver.opt_settings['linear_solver'] = 'mumps'\n",
+    "p.driver.declare_coloring()\n",
+    "\n",
+    "t = dm.Birkhoff(num_nodes=20)\n",
+    "\n",
+    "traj = p.model.add_subsystem('traj', dm.Trajectory())\n",
+    "phase = dm.Phase(ode_class=MoonLandingProblemODE, transcription=t)\n",
+    "\n",
+    "phase.set_time_options(fix_initial=True, fix_duration=False)\n",
+    "phase.add_state('h', fix_initial=True, rate_source='h_dot')\n",
+    "phase.add_state('v', fix_initial=True, rate_source='v_dot')\n",
+    "phase.add_state('m', fix_initial=True, lower=1e-3, rate_source='m_dot')\n",
+    "phase.add_control('T', lower=0.0, upper=1.227)\n",
+    "\n",
+    "phase.add_boundary_constraint('h', loc='final', equals=0.0)\n",
+    "phase.add_boundary_constraint('v', loc='final', equals=0.0)\n",
+    "\n",
+    "phase.add_objective('m', scaler=-1)\n",
+    "phase.set_simulate_options(atol=1.0E-1, rtol=1.0E-2)\n",
+    "\n",
+    "traj.add_phase('phase', phase)\n",
+    "\n",
+    "p.setup(check=True, force_alloc_complex=True)\n",
+    "\n",
+    "phase.set_time_val(initial=0.0, duration=1.0)\n",
+    "phase.set_state_val('h', [1.0, 0.0])\n",
+    "phase.set_state_val('v', [-0.783, 0.0])\n",
+    "phase.set_state_val('m', [1.0, 0.2])\n",
+    "phase.set_control_val('T', [0.0, 1.227])\n",
+    "dm.run_problem(p, simulate=False, simulate_kwargs={'times_per_seg': 100}, make_plots=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9b266ef1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from IPython.display import IFrame\n",
+    "IFrame(src=str(p.get_reports_dir() / 'traj_results_report.html'), width=800, height=1200)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "59334e2c",
+   "metadata": {},
+   "source": [
+    "### Notes on the solution\n",
+    "\n",
+    "We can see that the collocation solution accurately captures the jump in the thrust. The oscillatory behavior observed is a result of interpolation performed post solution rather than a property of the solution itself."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2e64faca",
+   "metadata": {},
+   "source": [
+    "## References\n",
+    "\n",
+    "```{bibliography}\n",
+    ":filter: docname in docnames\n",
+    "```"
+   ]
+  }
+ ],
+ "metadata": {
+  "celltoolbar": "Tags",
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}