Refactored the timeseries handling to be more general and less transcription-specific.

.github/workflows/dymos_docs_workflow.yml

@@ -215,14 +215,14 @@ jobs:
           mamba install jupyter-book jupyter -q -y
 
 
-      - name: Install optional dependencies
-        if: (github.event_name != 'workflow_dispatch' || matrix.NAME == 'latest') && matrix.OPTIONAL == '[all]'
+      - name: Install bokeh
+        if: (github.event_name != 'workflow_dispatch' || matrix.NAME == 'latest')
         shell: bash -l {0}
         run: |
           echo "============================================================="
           echo "Install additional packages for testing/coverage"
           echo "============================================================="
-          pip install bokeh
+          mamba install bokeh
 
       - name: Install Dymos
         if: (github.event_name != 'workflow_dispatch' || matrix.NAME == 'latest')

docs/dymos_book/_toc.yml

@@ -39,6 +39,7 @@ parts:
     - file: examples/robertson_problem/robertson_problem
     - file: examples/water_rocket/water_rocket
     - file: examples/cart_pole/cart_pole
+    - file: examples/brachistochrone/brachistochrone_tandem_phases
 - caption: Feature Reference
   chapters:
   - file: features/phases/phases_list

docs/dymos_book/examples/brachistochrone/brachistochrone_tandem_phases.ipynb

@@ -0,0 +1,365 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "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",
+   "metadata": {},
+   "source": [
+    "(examples:brachistochrone_tandem_phases)=\n",
+    "# Brachistochrone with tandem phases\n",
+    "\n",
+    "```{admonition} Things you'll learn through this example\n",
+    "- How to run two phases with different ODE's and different grids simultaneously in time.\n",
+    "```\n",
+    "\n",
+    "This is a contrived example but it demonstrates a useful feature of Dymos we call **tandem phases**.\n",
+    "Tandem phases are two phases that occur simultaneously in time (having the same start time and duration) but with different dynamics.\n",
+    "In practice, this can be useful when some of your dynamics are quite expensive and you can tolerate evaluating them on fewer nodes.\n",
+    "Or perhaps one phase has relatively rapid dynamics compared to the other one.  For instance, thermal responses tend to happen vary rapidly in an electric aircraft compared to changes in the flight dynamics state of the vehicle.\n",
+    "\n",
+    "In this example we'll evaulate the standard brachistochrone problem, but limit the arclength of the wire along which the bead travels.\n",
+    "The arclength is integrated as a state variable, and can be done so along with the typical states _x_, _y_, and _v_, but for the purposes of this contrive example we'll perform this integration in a separate phase that occurs at the same time."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The first phase to integrate the standard brachistochrone ODE\n",
+    "\n",
+    "```{admonition} Things to note about this phase\n",
+    "- The transcriptions for the two phases are delcared up front so that `tx1` may be used as both the transcription of the second phase, and for outputting the states of the first phase to the control input nodes of the second phase.\n",
+    "```\n",
+    "\n",
+    "This secondary timeseries is the key to making this sort of formulation work."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE\n",
+    "\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "plt.switch_backend('Agg')\n",
+    "import openmdao.api as om\n",
+    "import dymos as dm\n",
+    "\n",
+    "\n",
+    "p = om.Problem(model=om.Group())\n",
+    "\n",
+    "p.driver = om.pyOptSparseDriver()\n",
+    "p.driver.options['optimizer'] = 'SLSQP'\n",
+    "p.driver.options['print_results'] = False\n",
+    "p.driver.declare_coloring()\n",
+    "\n",
+    "# The transcription of the first phase\n",
+    "tx0 = dm.GaussLobatto(num_segments=10, order=3, compressed=False)\n",
+    "\n",
+    "# The transcription for the second phase (and the secondary timeseries outputs from the first phase)\n",
+    "tx1 = dm.Radau(num_segments=20, order=9, compressed=False)\n",
+    "\n",
+    "#\n",
+    "# First Phase: Integrate the standard brachistochrone ODE\n",
+    "#\n",
+    "phase0 = dm.Phase(ode_class=BrachistochroneODE, transcription=tx0)\n",
+    "\n",
+    "p.model.add_subsystem('phase0', phase0)\n",
+    "\n",
+    "phase0.set_time_options(fix_initial=True, duration_bounds=(.5, 10))\n",
+    "\n",
+    "phase0.add_state('x', fix_initial=True, fix_final=False)\n",
+    "\n",
+    "phase0.add_state('y', fix_initial=True, fix_final=False)\n",
+    "\n",
+    "phase0.add_state('v', fix_initial=True, fix_final=False)\n",
+    "\n",
+    "phase0.add_control('theta', continuity=True, rate_continuity=True,\n",
+    "                   units='deg', lower=0.01, upper=179.9)\n",
+    "\n",
+    "phase0.add_parameter('g', units='m/s**2', val=9.80665)\n",
+    "\n",
+    "phase0.add_boundary_constraint('x', loc='final', equals=10)\n",
+    "phase0.add_boundary_constraint('y', loc='final', equals=5)\n",
+    "\n",
+    "# Add alternative timeseries output to provide control inputs for the next phase\n",
+    "phase0.add_timeseries('timeseries2', transcription=tx1, subset='control_input')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The ODE for integrating the arc-length of the wire"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class BrachistochroneArclengthODE(om.ExplicitComponent):\n",
+    "\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('v', val=np.zeros(nn), desc='velocity', units='m/s')\n",
+    "        self.add_input('theta', val=np.zeros(nn), desc='angle of wire', units='rad')\n",
+    "        self.add_output('Sdot', val=np.zeros(nn), desc='rate of change of arclength', units='m/s')\n",
+    "\n",
+    "        # Setup partials\n",
+    "        arange = np.arange(nn)\n",
+    "\n",
+    "        self.declare_partials(of='Sdot', wrt='v', rows=arange, cols=arange)\n",
+    "        self.declare_partials(of='Sdot', wrt='theta', rows=arange, cols=arange)\n",
+    "\n",
+    "    def compute(self, inputs, outputs):\n",
+    "        theta = inputs['theta']\n",
+    "        v = inputs['v']\n",
+    "        outputs['Sdot'] = np.sqrt(1.0 + (1.0/np.tan(theta))**2) * v * np.sin(theta)\n",
+    "\n",
+    "    def compute_partials(self, inputs, jacobian):\n",
+    "        theta = inputs['theta']\n",
+    "        v = inputs['v']\n",
+    "        cos_theta = np.cos(theta)\n",
+    "        sin_theta = np.sin(theta)\n",
+    "        tan_theta = np.tan(theta)\n",
+    "        cot_theta = 1.0 / tan_theta\n",
+    "        csc_theta = 1.0 / sin_theta\n",
+    "\n",
+    "        jacobian['Sdot', 'v'] = sin_theta * np.sqrt(1.0 + cot_theta**2)\n",
+    "        jacobian['Sdot', 'theta'] = v * (cos_theta * (cot_theta**2 + 1) - cot_theta * csc_theta) / \\\n",
+    "            (np.sqrt(1 + cot_theta**2))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The second phase to integrate the arclength of the wire.\n",
+    "\n",
+    "```{admonition} Things to note about this phase\n",
+    "- Initial time and duration are input from those of the first phase (note the connections on line 16).\n",
+    "- _theta_ and _v_ are time-varying but determined by the first phase. Note the setting of `opt=False` and the connections on lines 18 and 19.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "phase1 = dm.Phase(ode_class=BrachistochroneArclengthODE, transcription=tx1)\n",
+    "\n",
+    "p.model.add_subsystem('phase1', phase1)\n",
+    "\n",
+    "phase1.set_time_options(fix_initial=True, input_duration=True)\n",
+    "\n",
+    "phase1.add_state('S', fix_initial=True, fix_final=False,\n",
+    "                 rate_source='Sdot', units='m')\n",
+    "\n",
+    "phase1.add_control('theta', opt=False, units='deg', targets='theta')\n",
+    "phase1.add_control('v', opt=False, units='m/s', targets='v')\n",
+    "\n",
+    "#\n",
+    "# Connect the two phases\n",
+    "#\n",
+    "p.model.connect('phase0.t_duration', 'phase1.t_duration')\n",
+    "\n",
+    "p.model.connect('phase0.timeseries2.controls:theta', 'phase1.controls:theta')\n",
+    "p.model.connect('phase0.timeseries2.states:v', 'phase1.controls:v')\n",
+    "\n",
+    "# Minimize time\n",
+    "phase1.add_objective('time', loc='final', ref=1)\n",
+    "\n",
+    "# Constraint the length of the wire.\n",
+    "phase1.add_boundary_constraint('S', loc='final', upper=11.5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Setup and run"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p.model.linear_solver = om.DirectSolver()\n",
+    "p.setup()\n",
+    "\n",
+    "p['phase0.t_initial'] = 0.0\n",
+    "p['phase0.t_duration'] = 2.0\n",
+    "\n",
+    "p.set_val('phase0.states:x', phase0.interp('x', ys=[0, 10]))\n",
+    "p.set_val('phase0.states:y', phase0.interp('y', ys=[10, 5]))\n",
+    "p.set_val('phase0.states:v', phase0.interp('v', ys=[0, 9.9]))\n",
+    "p.set_val('phase0.controls:theta', phase0.interp('theta', ys=[5, 100]))\n",
+    "\n",
+    "p['phase0.parameters:g'] = 9.80665\n",
+    "\n",
+    "p['phase1.states:S'] = 0.0\n",
+    "\n",
+    "res = dm.run_problem(p)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "remove-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from openmdao.utils.assert_utils import assert_near_equal\n",
+    "\n",
+    "assert not res\n",
+    "assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.85266, tolerance=1.0E-3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plots\n",
+    "\n",
+    "The following plots show the trajectory of the x, y, and v states (the top plot) and the trajectory of the arclength state (the bottom plot).  Note that these plots are linked, but use different grid spacings - the arclength is integrated on a significantly more dense grid.  This is enabled by the secondary timeseries output `timeseries2` in the first phase."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from bokeh.plotting import figure, show, output_notebook, output_file, save\n",
+    "\n",
+    "from bokeh.palettes import d3\n",
+    "from bokeh.resources import INLINE\n",
+    "from bokeh.models import Legend\n",
+    "\n",
+    "from IPython.display import HTML\n",
+    "\n",
+    "c = d3['Category10'][10]\n",
+    "i = np.array(0)\n",
+    "legend_contents = []\n",
+    "\n",
+    "sol_case = om.CaseReader('dymos_solution.db').get_case('final')\n",
+    "\n",
+    "sol_x = sol_case.get_val('phase0.timeseries.states:x')\n",
+    "sol_y = sol_case.get_val('phase0.timeseries.states:y')\n",
+    "sol_v = sol_case.get_val('phase0.timeseries.states:v')\n",
+    "sol_t0 = sol_case.get_val('phase0.timeseries.time')\n",
+    "sol_t1 = sol_case.get_val('phase1.timeseries.time')\n",
+    "sol_s = sol_case.get_val('phase1.timeseries.states:S')\n",
+    "\n",
+    "def add_plot(p, x, y, label, i):\n",
+    "    circle = p.circle(x.ravel(), y.ravel(), color=c[i], size=5)\n",
+    "    line = p.line(x.ravel(), y.ravel(), color=c[i])\n",
+    "    legend_contents.append((label, [circle, line]))\n",
+    "    i += 1\n",
+    "    \n",
+    "p1 = figure(plot_width=800, plot_height=300)\n",
+    "add_plot(p1, sol_t0, sol_x, 'x (m)', i)\n",
+    "add_plot(p1, sol_t0, sol_y, 'y (m)', i)\n",
+    "add_plot(p1, sol_t0, sol_v, 'v', i)\n",
+    "add_plot(p1, sol_t1, sol_s, 'arclength', i)\n",
+    "\n",
+    "p1.add_layout(Legend(items=legend_contents), 'right')\n",
+    "\n",
+    "p1.xaxis.axis_label = 'time (s)'\n",
+    "p1.yaxis.axis_label = 'state value'\n",
+    "p1.legend.location = 'bottom_right'\n",
+    "     \n",
+    "output_file('plot.html', mode='inline')\n",
+    "plot_file = save(p1)\n",
+    "\n",
+    "HTML(filename=plot_file)"
+   ]
+  }
+ ],
+ "metadata": {
+  "celltoolbar": "Slideshow",
+  "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.10.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}