Refactored the Radau transcription to avoid the use of StateIndependentsComp

.github/workflows/dymos_tests_workflow.yml

@@ -290,7 +290,7 @@ jobs:
             if [[ "${{ matrix.SNOPT }}" ]]; then
               echo "SNOPT ${{ matrix.SNOPT }} was requested but is not available on conda-forge"
             fi
-
+            conda install -c conda-forge lapack
             conda install -c conda-forge pyoptsparse
           else
             pip install git+https://github.com/OpenMDAO/build_pyoptsparse
@@ -446,9 +446,9 @@ jobs:
           testflo -b benchmark --pre_announce
           cd $HOME
           if [[ "${{ matrix.NAME }}" != "latest" ]]; then
-            testflo dymos -n 4 --pre_announce --show_skipped --durations 20 --coverage --coverpkg dymos
+            testflo dymos -n 1 --pre_announce --show_skipped --durations 20 --coverage --coverpkg dymos
           else
-            testflo dymos -n 4 --pre_announce --show_skipped --durations 20
+            testflo dymos -n 1 --pre_announce --show_skipped --durations 20
           fi
       - name: Submit coverage
         if: github.event_name != 'workflow_dispatch'

docs/dymos_book/examples/brachistochrone/brachistochrone_tandem_phases.ipynb

@@ -264,9 +264,8 @@
     "phase0.set_state_val('v', [0, 9.9])\n",
     "phase0.set_control_val('theta', [5, 100.5])\n",
     "\n",
-    "p['phase0.parameters:g'] = 9.80665\n",
-    "\n",
-    "p['phase1.states:S'] = 0.0\n",
+    "phase0.set_parameter_val('g', 9.80665)\n",
+    "phase1.set_state_val('S', [0.0, 10.0])\n",
     "\n",
     "res = dm.run_problem(p)"
    ]
@@ -284,8 +283,8 @@
    "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)"
+    "assert res.success\n",
+    "assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.85, tolerance=5.0E-3)"
    ]
   },
   {
@@ -320,10 +319,6 @@
     "i = np.array(0)\n",
     "legend_contents = []\n",
     "\n",
-    "print(p.get_outputs_dir())\n",
-    "import os\n",
-    "print(os.listdir(os.getcwd()))\n",
-    "\n",
     "sol_case = om.CaseReader(p.get_outputs_dir() / 'dymos_solution.db').get_case('final')\n",
     "\n",
     "sol_x = sol_case.get_val('phase0.timeseries.x')\n",
@@ -334,7 +329,7 @@
     "sol_s = sol_case.get_val('phase1.timeseries.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",
+    "    circle = p.scatter(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",
@@ -372,7 +367,7 @@
  "metadata": {
   "celltoolbar": "Slideshow",
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "py311_2",
    "language": "python",
    "name": "python3"
   },

docs/dymos_book/examples/hull/hull_problem.ipynb

@@ -97,7 +97,7 @@
     "We seek to minimize the control effort required and minimize the displacement from the origin.\n",
     "\n",
     "\\begin{align}\n",
-    "    \\mathrm{Minimize} \\, J &= 2.5x_f^2 \\, + \\, 0.5 \\int_0^1 u^2 dt\n",
+    "    \\mathrm{Minimize} \\, J &= 2.5x_f^2 \\, + \\, 0.5 \\int_0^{10} u^2 dt\n",
     "\\end{align}\n",
     "\n",
     "Subject to the initial conditions\n",
@@ -182,7 +182,6 @@
     "        nn = self.options['num_nodes']\n",
     "\n",
     "        self.add_input('u', val=np.zeros(nn), desc='control')\n",
-    "\n",
     "        self.add_output('L', val=np.zeros(nn), desc='Lagrangian', units='1/s')\n",
     "\n",
     "        # Setup partials\n",
@@ -195,6 +194,7 @@
     "\n",
     "    def compute_partials(self, inputs, partials, discrete_inputs=None):\n",
     "        u = inputs['u']\n",
+    "\n",
     "        partials['L', 'u'] = u\n"
    ]
   },
@@ -296,7 +296,7 @@
     "# Create a trajectory and add a phase to it\n",
     "#\n",
     "traj = p.model.add_subsystem('traj', dm.Trajectory())\n",
-    "tx = transcription=dm.Radau(num_segments=24)\n",
+    "tx = dm.Birkhoff(num_nodes=60)\n",
     "phase = traj.add_phase('phase0', dm.Phase(ode_class=HullProblemODE, transcription=tx))"
    ]
   },
@@ -319,7 +319,7 @@
     "We would then need to put a boundary constraint in place to enforce their final values.\n",
     "Feel free to experiment with different ways of enforcing the boundary constraints on this problem and see how it affects performance.\n",
     "\n",
-    "The scaler values (`ref`) are all set to 1 here.\n",
+    "Experimentation has shown that scaling the control up (the reference value is 1/100) and scaling the objective up seems to converge to a beter solution.\n",
     "\n",
     "Also, we don't need to specify targets for any of the variables here because their names _are_ the targets in the top-level of the model.\n",
     "The rate source and units for the states are obtained from the tags in the ODE component we previously defined.\n"
@@ -341,8 +341,10 @@
     "\n",
     "phase.set_time_options(fix_initial=True, fix_duration=True)\n",
     "phase.add_state('x', fix_initial=True, fix_final=False, rate_source='u')\n",
-    "phase.add_state('xL', fix_initial=True, fix_final=False, rate_source='L')\n",
-    "phase.add_control('u', opt=True, targets=['u'])"
+    "phase.add_state('x_L', fix_initial=True, fix_final=False, rate_source='L')\n",
+    "phase.add_control('u', opt=True, ref=0.01, rate_continuity=True, rate2_continuity=True)\n",
+    "\n",
+    "phase.simulate_options['times_per_seg'] = 100"
    ]
   },
   {
@@ -361,7 +363,7 @@
     "#\n",
     "# Minimize time at the end of the phase\n",
     "#\n",
-    "phase.add_objective('J = 2.5*x**2 + xL')\n",
+    "phase.add_objective('J = 2.5 * x**2 + x_L', ref=0.01)\n",
     "\n",
     "#\n",
     "# Setup the Problem\n",
@@ -386,11 +388,7 @@
     "It will automatically record the final solution achieved by the optimizer in case named `'final'` in a file called `dymos_solution.db`.\n",
     "By specifying `simulate=True`, it will automatically follow the solution with an explicit integration using `scipy.solve_ivp`.\n",
     "The results of the simulation are stored in a case named `final` in the file `dymos_simulation.db`.\n",
-    "This explicit simulation demonstrates how the system evolved with the given controls, and serves as a check that we're using a dense enough grid (enough segments and segments of sufficient order) to accurately represent the solution.\n",
-    "\n",
-    "If those two solution didn't agree reasonably well, we could rerun the problem with a more dense grid.\n",
-    "Instead, we're asking Dymos to automatically change the grid if necessary by specifying `refine_method='ph'`.\n",
-    "This will attempt to repeatedly solve the problem and change the number of segments and segment orders until the solution is in reasonable agreement."
+    "This explicit simulation demonstrates how the system evolved with the given controls, and serves as a check that we're using a dense enough grid (enough segments and segments of sufficient order) to accurately represent the solution."
    ]
   },
   {
@@ -413,16 +411,52 @@
     "# Set the initial values\n",
     "#\n",
     "phase.set_state_val('x', [1.5, 1])\n",
-    "phase.set_state_val('xL', [0, 1])\n",
+    "phase.set_state_val('x_L', [0, 1])\n",
     "phase.set_time_val(initial=0.0, duration=10.0)\n",
-    "phase.set_control_val('u', [-7, -0.14])\n",
+    "phase.set_control_val('u', [-0.1, -0.15])\n",
     "\n",
     "#\n",
     "# Solve for the optimal trajectory\n",
     "#\n",
     "dm.run_problem(p, run_driver=True, simulate=True)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5fdd138d",
+   "metadata": {
+    "tags": [
+     "remove-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# # Testing the result\n",
+    "from openmdao.utils.assert_utils import assert_near_equal\n",
+    "c = 5\n",
+    "\n",
+    "def solution(x0, td):\n",
+    "    \"\"\"\n",
+    "    The analytic solution to the Hull problem.\n",
+    "    \"\"\"\n",
+    "    xf = x0 - c * x0 * td / (1 + c * td)\n",
+    "    uf = -c * x0 / (1 + c * td)\n",
+    "\n",
+    "    return xf, uf\n",
+    "\n",
+    "xf, uf = solution(1.5, 10)\n",
+    "\n",
+    "assert_near_equal(p.get_val('traj.phase0.timeseries.x')[-1],\n",
+    "                  xf,\n",
+    "                  tolerance=1e-2)\n",
+    "\n",
+    "assert_near_equal(p.get_val('traj.phase0.timeseries.u')[-1],\n",
+    "                  uf,\n",
+    "                  tolerance=1e-1)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "3f9a5d44",
@@ -442,8 +476,7 @@
    "execution_count": null,
    "id": "beea4c6a",
    "metadata": {
-    "id": "beea4c6a",
-    "scrolled": true
+    "id": "beea4c6a"
    },
    "outputs": [],
    "source": [
@@ -452,12 +485,12 @@
     "\n",
     "t = sol.get_val('traj.phase0.timeseries.time')\n",
     "x = sol.get_val('traj.phase0.timeseries.x')\n",
-    "xL = sol.get_val('traj.phase0.timeseries.xL')\n",
+    "xL = sol.get_val('traj.phase0.timeseries.x_L')\n",
     "u = sol.get_val('traj.phase0.timeseries.u')\n",
     "\n",
     "t_sim = sim.get_val('traj.phase0.timeseries.time')\n",
     "x_sim = sim.get_val('traj.phase0.timeseries.x')\n",
-    "xL_sim = sim.get_val('traj.phase0.timeseries.xL')\n",
+    "xL_sim = sim.get_val('traj.phase0.timeseries.x_L')\n",
     "u_sim = sim.get_val('traj.phase0.timeseries.u')\n",
     "\n",
     "fig = plt.figure(constrained_layout=True, figsize=(12, 4))\n",
@@ -468,7 +501,7 @@
     "u_ax = fig.add_subplot(gs[2, 0])\n",
     "\n",
     "x_ax.set_ylabel('x ($m$)')\n",
-    "xL_ax.set_ylabel('xL ($m/s$)')\n",
+    "xL_ax.set_ylabel('x_L ($m/s$)')\n",
     "u_ax.set_ylabel('u ($m/s^2$)')\n",
     "xL_ax.set_xlabel('t (s)')\n",
     "u_ax.set_xlabel('t (s)')\n",
@@ -480,6 +513,7 @@
     "x_sim_handle, = x_ax.plot(t_sim, x_sim, '-', ms=1)\n",
     "xL_ax.plot(t_sim, xL_sim, '-', ms=1)\n",
     "u_ax.plot(t_sim, u_sim, '-', ms=1)\n",
+    "u_ax.set_ylim(-0.2, -0.1)\n",
     "\n",
     "for ax in [x_ax, xL_ax, u_ax]:\n",
     "    ax.grid(True, alpha=0.2)\n",
@@ -508,7 +542,7 @@
    "provenance": []
   },
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "py311_2",
    "language": "python",
    "name": "python3"
   },
@@ -522,7 +556,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,

docs/dymos_book/examples/implicit_duration/brachistochrone_implicit_duration.ipynb

@@ -74,6 +74,7 @@
    "source": [
     "import jax.numpy as jnp\n",
     "import openmdao.api as om\n",
+    "import dymos as dm\n",
     "\n",
     "\n",
     "class BrachistochroneODE(om.JaxExplicitComponent):\n",
@@ -291,7 +292,7 @@
     "# Setup the problem\n",
     "p.setup()\n",
     "\n",
-    "phase.set_time_val(initial=0.0, duration=1.8)\n",
+    "phase.set_time_val(initial=0.0, duration=2.0)\n",
     "phase.set_state_val('x', vals=[0, 10])\n",
     "phase.set_state_val('y', vals=[10, 5])\n",
     "phase.set_state_val('v', vals=[0.1, 100])\n",
@@ -412,7 +413,7 @@
     "# Setup the problem\n",
     "p.setup()\n",
     "\n",
-    "phase.set_time_val(initial=0.0, duration=1.8)\n",
+    "phase.set_time_val(initial=0.0, duration=2.0)\n",
     "phase.set_state_val('x', vals=[0, 10])\n",
     "phase.set_state_val('y', vals=[10, 5])\n",
     "phase.set_state_val('v', vals=[0., 10])\n",
@@ -450,6 +451,23 @@
     "print('final theta', jnp.degrees(theta[-1]))"
    ]
   },
+  {
+   "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_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)\n",
+    "assert_near_equal(p.get_val('traj.phase0.timeseries.y')[-1], 5.0, tolerance=1.0E-3)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -475,7 +493,7 @@
    }
   },
   "kernelspec": {
-   "display_name": "py311_2",
+   "display_name": "py312",
    "language": "python",
    "name": "python3"
   },
@@ -489,7 +507,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.12.8"
   }
  },
  "nbformat": 4,