Added interpolation of results to arbitrary nodes, removed use of `sc…

…ipy.interpolate.interp1d` from `phase.interp`. (#1129)

* Added ability to retrieve interpolated values in phase on an arbitrary set of nodes.

* Use make_interp_spline for backward compatibility

* Fixed sorting issue when using np.unique to remove duplicate time values in interpolation.

* optional [all] seems to fix no_snopt

* raise on invalid desvars in orbit raise problem

* reset invalid_desvar_behavior

* bumping up number of testflo procs to 4.

* removed old upgrade guide

* balancing lagrange multipliers

* water rocket for docs failing latest.

* removed birkhoff from racecar test due to slowness

* change interpolant to make_interp_spline

* switch implicit duration example to use BoundsEnforceLS

* testing implicit duration with SLSQP

*  two burn orbit raise mpi example radau transcription instead of gl

.github/workflows/dymos_tests_workflow.yml

@@ -126,7 +126,8 @@ jobs:
             PETSc: 3.21.0
             PYOPTSPARSE: 'v2.11.0'
             OPENMDAO: 'latest'
-            OPTIONAL: '[test]'
+            OPTIONAL: '[all]'
+            JAX: '0.4.28'
             EXCLUDE: ${{ github.event_name == 'workflow_dispatch'  && ! inputs.no_snopt }}
 
           # baseline versions except no MPI/PETSc
@@ -149,7 +150,7 @@ jobs:
             PYOPTSPARSE: 'latest'
             SNOPT: 7.7
             OPENMDAO: 'dev'
-            OPTIONAL: '[test]'
+            OPTIONAL: '[all]'
             JAX: 'latest'
             EXCLUDE: ${{ github.event_name == 'workflow_dispatch'  && ! inputs.latest }}
 
@@ -442,11 +443,10 @@ jobs:
           testflo -b benchmark --pre_announce
           cd $HOME
           if [[ "${{ matrix.NAME }}" != "latest" ]]; then
-            testflo dymos -n 2 --pre_announce --show_skipped --durations 20 --coverage --coverpkg dymos
+            testflo dymos -n 4 --pre_announce --show_skipped --durations 20 --coverage --coverpkg dymos
           else
-            testflo dymos -n 2 --pre_announce --show_skipped --durations 20
+            testflo dymos -n 4 --pre_announce --show_skipped --durations 20
           fi
-
       - name: Submit coverage
         if: github.event_name != 'workflow_dispatch'
         continue-on-error: true

docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb

@@ -100,9 +100,9 @@
    "outputs": [],
    "source": [
     "import numpy as np\n",
+    "from scipy.interpolate import make_interp_spline\n",
     "\n",
     "import openmdao.api as om\n",
-    "from openmdao.components.interp_util.interp import InterpND\n",
     "\n",
     "import dymos as dm\n",
     "from dymos.models.atmosphere.atmos_1976 import USatm1976Data\n",
@@ -205,12 +205,11 @@
     "\n",
     "        alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]\n",
     "        rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]\n",
-    "        self.rho_interp = InterpND(points=np.array(alt_data),\n",
-    "                                   values=np.array(rho_data),\n",
-    "                                   method='slinear')\n",
-    "        \n",
-    "    def compute(self, inputs, outputs):\n",
+    "        self.rho_interp = make_interp_spline(np.array(alt_data),\n",
+    "                                             np.array(rho_data),\n",
+    "                                             k=1)\n",
     "\n",
+    "    def compute(self, inputs, outputs):\n",
     "        gam = inputs['gam']\n",
     "        v = inputs['v']\n",
     "        h = inputs['h']\n",
@@ -220,9 +219,9 @@
     "\n",
     "        GRAVITY = 9.80665  # m/s**2\n",
     "\n",
-    "        rho = self.rho_interp.interpolate(h)\n",
+    "        rho = self.rho_interp(h)\n",
     "\n",
-    "        q = 0.5*rho*inputs['v']**2\n",
+    "        q = 0.5*rho*v**2\n",
     "        qS = q * S\n",
     "        D = qS * CD\n",
     "        cgam = np.cos(gam)\n",
@@ -324,7 +323,7 @@
     "# The units of the states are automatically inferred by multiplying the units\n",
     "# of those rates by the time units.\n",
     "phase.add_state('r', fix_initial=True, solve_segments='forward')\n",
-    "phase.add_state('h', fix_initial=True, solve_segments='forward')\n",
+    "phase.add_state('h', fix_initial=True, lower=0.0, solve_segments='forward')\n",
     "phase.add_state('gam', fix_initial=False, solve_segments='forward')\n",
     "phase.add_state('v', fix_initial=False, solve_segments='forward')\n",
     "\n",
@@ -453,7 +452,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.6"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,

docs/dymos_book/examples/water_rocket/water_rocket.ipynb

@@ -601,9 +601,11 @@
     "traj = p.model.add_subsystem('traj', traj)\n",
     "\n",
     "p.driver = om.pyOptSparseDriver(optimizer='IPOPT', print_results=False)\n",
-    "p.driver.opt_settings['print_level'] = 4\n",
+    "p.driver.opt_settings['print_level'] = 5\n",
     "p.driver.opt_settings['max_iter'] = 1000\n",
-    "p.driver.opt_settings['mu_strategy'] = 'monotone'\n",
+    "p.driver.opt_settings['mu_strategy'] = 'adaptive'\n",
+    "p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'  # for faster convergence\n",
+    "p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'\n",
     "p.driver.declare_coloring(tol=1.0E-12)\n",
     "\n",
     "# Finish Problem Setup\n",
@@ -718,9 +720,11 @@
     "traj = p.model.add_subsystem('traj', traj)\n",
     "\n",
     "p.driver = om.pyOptSparseDriver(optimizer='IPOPT')\n",
-    "p.driver.opt_settings['print_level'] = 4\n",
+    "p.driver.opt_settings['print_level'] = 5\n",
     "p.driver.opt_settings['max_iter'] = 1000\n",
-    "p.driver.opt_settings['mu_strategy'] = 'monotone'\n",
+    "p.driver.opt_settings['mu_strategy'] = 'adaptive'\n",
+    "p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'  # for faster convergence\n",
+    "p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'\n",
     "p.driver.declare_coloring(tol=1.0E-12)\n",
     "\n",
     "# Finish Problem Setup\n",
@@ -822,13 +826,6 @@
     ":filter: docname in docnames\n",
     "```"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -842,7 +839,7 @@
    }
   },
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "py311_2",
    "language": "python",
    "name": "python3"
   },
@@ -856,7 +853,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,

docs/dymos_book/features/phases/timeseries.ipynb

@@ -101,6 +101,162 @@
     "        :noindex:\n",
     "```"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Interpolating Timeseries Outputs\n",
+    "\n",
+    "Outputs from the timeseries are provided at the node resolution used to solve the problem.\n",
+    "Sometimes it is useful to visualize the data at a different resolution - for instance to see oscillatory behavior in the solution that may not be noticeable by only viewing the solution nodes.\n",
+    "\n",
+    "The Phase.interp method allows nodes to be provided as a sequence in the \"tau space\" (nondimensional time) of the phase. In phase tau space,\n",
+    "the time-like variable is -1 at the start of the phase and +1 at the end of the phase.\n",
+    "\n",
+    "For instance, let's say we want to solve the minimum time climb example and plot the control solution at 100 nodes equally spaced throughout the phase."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from dymos.examples.min_time_climb.min_time_climb_ode import MinTimeClimbODE\n",
+    "\n",
+    "p = om.Problem(model=om.Group())\n",
+    "\n",
+    "p.driver = om.pyOptSparseDriver()\n",
+    "p.driver.options['optimizer'] = 'IPOPT'\n",
+    "p.driver.options['print_results'] = 'minimal'\n",
+    "p.driver.declare_coloring()\n",
+    "\n",
+    "p.driver.opt_settings['tol'] = 1.0E-5\n",
+    "p.driver.opt_settings['print_level'] = 0\n",
+    "p.driver.opt_settings['mu_strategy'] = 'monotone'\n",
+    "p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'\n",
+    "p.driver.opt_settings['mu_init'] = 0.01\n",
+    "\n",
+    "tx = dm.Radau(num_segments=8, order=3)\n",
+    "\n",
+    "traj = dm.Trajectory()\n",
+    "\n",
+    "phase = dm.Phase(ode_class=MinTimeClimbODE, transcription=tx)\n",
+    "traj.add_phase('phase0', phase)\n",
+    "\n",
+    "p.model.add_subsystem('traj', traj)\n",
+    "\n",
+    "phase.set_time_options(fix_initial=True, duration_bounds=(50, 400),\n",
+    "                       duration_ref=100.0)\n",
+    "\n",
+    "phase.add_state('r', fix_initial=True, lower=0, upper=1.0E6,\n",
+    "                ref=1.0E3, defect_ref=1.0E3, units='m',\n",
+    "                rate_source='flight_dynamics.r_dot')\n",
+    "\n",
+    "phase.add_state('h', fix_initial=True, lower=0, upper=20000.0,\n",
+    "                ref=20_000, defect_ref=20_000, units='m',\n",
+    "                rate_source='flight_dynamics.h_dot', targets=['h'])\n",
+    "\n",
+    "phase.add_state('v', fix_initial=True, lower=10.0,\n",
+    "                ref=1.0E2, defect_ref=1.0E2, units='m/s',\n",
+    "                rate_source='flight_dynamics.v_dot', targets=['v'])\n",
+    "\n",
+    "phase.add_state('gam', fix_initial=True, lower=-1.5, upper=1.5,\n",
+    "                ref=1.0, defect_ref=1.0, units='rad',\n",
+    "                rate_source='flight_dynamics.gam_dot', targets=['gam'])\n",
+    "\n",
+    "phase.add_state('m', fix_initial=True, lower=10.0, upper=1.0E5,\n",
+    "                ref=10_000, defect_ref=10_000, units='kg',\n",
+    "                rate_source='prop.m_dot', targets=['m'])\n",
+    "\n",
+    "phase.add_control('alpha', units='deg', lower=-8.0, upper=8.0, scaler=1.0,\n",
+    "                    rate_continuity=True, rate_continuity_scaler=100.0,\n",
+    "                    rate2_continuity=False, targets=['alpha'])\n",
+    "\n",
+    "phase.add_parameter('S', val=49.2386, units='m**2', opt=False, targets=['S'])\n",
+    "phase.add_parameter('Isp', val=1600.0, units='s', opt=False, targets=['Isp'])\n",
+    "phase.add_parameter('throttle', val=1.0, opt=False, targets=['throttle'])\n",
+    "\n",
+    "phase.add_boundary_constraint('h', loc='final', equals=20000, scaler=1.0E-3)\n",
+    "phase.add_boundary_constraint('aero.mach', loc='final', equals=1.0)\n",
+    "phase.add_boundary_constraint('gam', loc='final', equals=0.0)\n",
+    "\n",
+    "phase.add_path_constraint(name='h', lower=100.0, upper=20000, ref=20000)\n",
+    "phase.add_path_constraint(name='aero.mach', lower=0.1, upper=1.8)\n",
+    "\n",
+    "# Minimize time at the end of the phase\n",
+    "phase.add_objective('time', loc='final', ref=1.0)\n",
+    "\n",
+    "# test mixing wildcard ODE variable expansion and unit overrides\n",
+    "phase.add_timeseries_output(['aero.*', 'prop.thrust', 'prop.m_dot'],\n",
+    "                            units={'aero.f_lift': 'lbf', 'prop.thrust': 'lbf'})\n",
+    "\n",
+    "p.model.linear_solver = om.DirectSolver()\n",
+    "\n",
+    "p.setup(check=True)\n",
+    "\n",
+    "p['traj.phase0.t_initial'] = 0.0\n",
+    "p['traj.phase0.t_duration'] = 350.0\n",
+    "\n",
+    "phase.set_time_val(initial=0.0, duration=350.0)\n",
+    "phase.set_state_val('r', [0.0, 111319.54])\n",
+    "phase.set_state_val('h', [100.0, 20000.0])\n",
+    "phase.set_state_val('v', [135.964, 283.159])\n",
+    "phase.set_state_val('gam', [0.0, 0.0])\n",
+    "phase.set_state_val('m', [19030.468, 16841.431])\n",
+    "phase.set_control_val('alpha', [0.0, 0.0])\n",
+    "\n",
+    "dm.run_problem(p)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we pull the values of time and alpha from the timeseries.\n",
+    "Remember that the timeseries provides values at the nodes only. The fitting polynomials may exhibit oscillations between these nodes.\n",
+    "\n",
+    "We use cubic interpolation to interpolate time and alpha onto 100 evenly distributed nodes across the phase, and then plot the results."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "\n",
+    "t = p.get_val('traj.phase0.timeseries.time')\n",
+    "alpha = p.get_val('traj.phase0.timeseries.alpha', units='deg')\n",
+    "\n",
+    "t_100 = phase.interp(xs=t, ys=t,\n",
+    "                     nodes=np.linspace(-1, 1, 100),\n",
+    "                     kind='cubic')\n",
+    "alpha_100 = phase.interp(xs=t, ys=alpha,\n",
+    "                         nodes=np.linspace(-1, 1, 100),\n",
+    "                         kind='cubic')\n",
+    "\n",
+    "# Plot the solution nodes\n",
+    "%matplotlib inline\n",
+    "plt.plot(t, alpha, 'o')\n",
+    "plt.plot(t_100, alpha_100, '-')\n",
+    "plt.xlabel('time (s)')\n",
+    "plt.ylabel('angle of attack (deg)')\n",
+    "plt.grid()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Both the `simulate` method and the `interp` method can provide more dense output. The difference is that simulate is effectively finding a different solution for the states, by starting from the initial values and propagating the ODE forward in time using the interpolated control splines as inputs. The `interp` method is making some assumptions about the interpolation (for instance, cubic interpolation in this case), so it may not be a perfectly accurate representation of the underlying interpolating polynomials."
+   ]
   }
  ],
  "metadata": {
@@ -127,7 +283,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.5"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,