diff --git a/.github/workflows/dymos_docs_workflow.yml b/.github/workflows/dymos_docs_workflow.yml index ef0c85f09..204b97b49 100644 --- a/.github/workflows/dymos_docs_workflow.yml +++ b/.github/workflows/dymos_docs_workflow.yml @@ -15,7 +15,7 @@ on: jobs: test_ubuntu: - runs-on: ubuntu-latest + runs-on: ubuntu-22.04 timeout-minutes: 90 @@ -102,7 +102,12 @@ jobs: echo "Install jax" echo "=============================================================" if [[ "${{ matrix.JAX }}" == "latest" ]]; then - python -m pip install jaxlib jax + { + python -m pip install jaxlib + } || { + echo "Failed to install jaxlib..." + } + python -m pip install jax else python -m pip install jaxlib==${{ matrix.JAX }} jax==${{ matrix.JAX }} fi diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml index 9c4fcff46..0d2503498 100644 --- a/.github/workflows/dymos_tests_workflow.yml +++ b/.github/workflows/dymos_tests_workflow.yml @@ -77,7 +77,7 @@ run-name: ${{ inputs.run_name }} jobs: test_ubuntu: - runs-on: ubuntu-latest + runs-on: ubuntu-22.04 timeout-minutes: 90 @@ -241,7 +241,12 @@ jobs: echo "Install jax" echo "=============================================================" if [[ "${{ matrix.JAX }}" == "latest" ]]; then - python -m pip install jaxlib jax + { + python -m pip install jaxlib + } || { + echo "Failed to install jaxlib..." + } + python -m pip install jax else python -m pip install jaxlib==${{ matrix.JAX }} jax==${{ matrix.JAX }} fi diff --git a/docs/dymos_book/_toc.yml b/docs/dymos_book/_toc.yml index c1c2fb882..8a3984034 100644 --- a/docs/dymos_book/_toc.yml +++ b/docs/dymos_book/_toc.yml @@ -42,7 +42,7 @@ parts: - file: examples/water_rocket/water_rocket - file: examples/cart_pole/cart_pole - file: examples/brachistochrone/brachistochrone_tandem_phases - - file: examples/cannonball_implicit_duration/cannonball_implicit_duration + - file: examples/implicit_duration/brachistochrone_implicit_duration - caption: Feature Reference chapters: - file: features/phases/phases_list diff --git a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb deleted file mode 100644 index 89c99afb2..000000000 --- a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb +++ /dev/null @@ -1,460 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "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 # noqa: F401\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 # noqa: F401\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 # noqa: F401\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": [ - "# Cannonball Implicit Duration\n", - "\n", - "This example demonstrates determining the duration of a phase using a \n", - "nonlinear solver to satisfy an endpoint condition. As in the multi-phase\n", - "cannonball problem, we will simultaneously find the optimal design for \n", - "the cannonball (its radius) and the optimal flight profile (its launch\n", - "angle). However, here we will use a single phase that terminates at the\n", - "condition\n", - "\n", - "\\begin{align}\n", - " h(t_f) = 0\n", - "\\end{align}\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## The cannonball sizing component\n", - "\n", - "This component computes the area (needed to compute drag) and mass (needed to compute energy) of a cannonball of a given radius and density.\n", - "\n", - "This component sits upstream of the trajectory model and feeds its outputs to the trajectory as parameters." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "active-ipynb", - "remove-input", - "remove-output" - ] - }, - "outputs": [], - "source": [ - "%matplotlib inline" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": [ - "output_scroll" - ] - }, - "outputs": [], - "source": [ - "import numpy as np\n", - "from scipy.interpolate import make_interp_spline\n", - "\n", - "import openmdao.api as om\n", - "\n", - "import dymos as dm\n", - "from dymos.models.atmosphere.atmos_1976 import USatm1976Data\n", - "\n", - "#############################################\n", - "# Component for the design part of the model\n", - "#############################################\n", - "class CannonballSizeComp(om.ExplicitComponent):\n", - " \"\"\"\n", - " Compute the area and mass of a cannonball with a given radius and density.\n", - "\n", - " Notes\n", - " -----\n", - " This component is not vectorized with 'num_nodes' as is the usual way\n", - " with Dymos, but is instead intended to compute a scalar mass and reference\n", - " area from scalar radius and density inputs. This component does not reside\n", - " in the ODE but instead its outputs are connected to the trajectory via\n", - " input design parameters.\n", - " \"\"\"\n", - " def setup(self):\n", - " self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m')\n", - " self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3')\n", - "\n", - " self.add_output(name='mass', shape=(1,), desc='cannonball mass', units='kg')\n", - " self.add_output(name='S', shape=(1,), desc='aerodynamic reference area', units='m**2')\n", - "\n", - " self.declare_partials(of='mass', wrt='dens')\n", - " self.declare_partials(of='mass', wrt='radius')\n", - "\n", - " self.declare_partials(of='S', wrt='radius')\n", - "\n", - " def compute(self, inputs, outputs):\n", - " radius = inputs['radius']\n", - " dens = inputs['dens']\n", - "\n", - " outputs['mass'] = (4/3.) * dens * np.pi * radius ** 3\n", - " outputs['S'] = np.pi * radius ** 2\n", - "\n", - " def compute_partials(self, inputs, partials):\n", - " radius = inputs['radius']\n", - " dens = inputs['dens']\n", - "\n", - " partials['mass', 'dens'] = (4/3.) * np.pi * radius ** 3\n", - " partials['mass', 'radius'] = 4. * dens * np.pi * radius ** 2\n", - "\n", - " partials['S', 'radius'] = 2 * np.pi * radius\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## The cannonball ODE component\n", - "\n", - "This component computes the state rates and the kinetic energy of the cannonball.\n", - "By calling the `declare_coloring` method wrt all inputs and using method `'fd'`, we're telling OpenMDAO to automatically determine the sparsity pattern of the outputs with respect to the inputs, **and** to automatically compute those outputs using a finite-difference approximation." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "class CannonballODE(om.ExplicitComponent):\n", - " \"\"\"\n", - " Cannonball ODE assuming flat earth and accounting for air resistance\n", - " \"\"\"\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", - " # static parameters\n", - " self.add_input('m', units='kg')\n", - " self.add_input('S', units='m**2')\n", - " # 0.5 good assumption for a sphere\n", - " self.add_input('CD', 0.5)\n", - "\n", - " # time varying inputs\n", - " self.add_input('h', units='m', shape=nn)\n", - " self.add_input('v', units='m/s', shape=nn)\n", - " self.add_input('gam', units='rad', shape=nn)\n", - "\n", - " # state rates\n", - " self.add_output('v_dot', shape=nn, units='m/s**2', tags=['dymos.state_rate_source:v'])\n", - " self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam'])\n", - " self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h'])\n", - " self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])\n", - " self.add_output('ke', shape=nn, units='J')\n", - "\n", - " # Ask OpenMDAO to compute the partial derivatives using finite-difference\n", - " # with a partial coloring algorithm for improved performance, and use\n", - " # a graph coloring algorithm to automatically detect the sparsity pattern.\n", - " self.declare_coloring(wrt='*', method='fd')\n", - "\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 = 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", - " m = inputs['m']\n", - " S = inputs['S']\n", - " CD = inputs['CD']\n", - "\n", - " GRAVITY = 9.80665 # m/s**2\n", - "\n", - " rho = self.rho_interp(h)\n", - "\n", - " q = 0.5*rho*v**2\n", - " qS = q * S\n", - " D = qS * CD\n", - " cgam = np.cos(gam)\n", - " sgam = np.sin(gam)\n", - " outputs['v_dot'] = - D/m-GRAVITY*sgam\n", - " outputs['gam_dot'] = -(GRAVITY/v)*cgam\n", - " outputs['h_dot'] = v*sgam\n", - " outputs['r_dot'] = v*cgam\n", - " outputs['ke'] = 0.5*m*v**2" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Generating an initial guess\n", - "\n", - "Since we will be using a nonlinear solver to converge the dynamics and the duration, a simple linear initial guess will not be sufficient. The initial guess we will use is generated from the kinematics equations that assume no atmospheric drag. We compute the state values for some duration and initial flight path angle\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "def initial_guess(t_dur, gam0=np.pi/3):\n", - " t = np.linspace(0, t_dur, num=100)\n", - " g = -9.80665\n", - " v0 = -0.5*g*t_dur/np.sin(gam0)\n", - " r = v0*np.cos(gam0)*t\n", - " h = v0*np.sin(gam0)*t + 0.5*g*t**2\n", - " v = np.sqrt(v0*np.cos(gam0)**2 + (v0*np.sin(gam0) + g*t)**2)\n", - " gam = np.arctan2(v0*np.sin(gam0) + g*t, v0*np.cos(gam0)**2)\n", - "\n", - " guess = {'t': t, 'r': r, 'h': h, 'v': v, 'gam': gam}\n", - "\n", - " return guess\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Building and running the problem\n", - "\n", - "The following code defines the components for the physical\n", - "cannonball calculations and ODE problem, sets up trajectory using a\n", - "single phase, and specifies the stopping condition for the phase.\n", - "The initial flight path angle is free, since 45 degrees is not\n", - "necessarily optimal once air resistance is taken into account." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "p = om.Problem(name='cannonball_implicit_duration', model=om.Group())\n", - "\n", - "p.driver = om.pyOptSparseDriver()\n", - "p.driver.options['optimizer'] = 'IPOPT'\n", - "p.driver.declare_coloring()\n", - "\n", - "p.driver.opt_settings['derivative_test'] = 'first-order'\n", - "p.driver.opt_settings['mu_strategy'] = 'monotone'\n", - "p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'\n", - "p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'\n", - "p.driver.opt_settings['mu_init'] = 0.01\n", - "p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'\n", - "\n", - "p.set_solver_print(level=0, depth=99)\n", - "\n", - "p.model.add_subsystem('size_comp', CannonballSizeComp(),\n", - " promotes_inputs=['radius', 'dens'])\n", - "p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')\n", - "p.model.add_design_var('radius', lower=0.01, upper=0.10,\n", - " ref0=0.01, ref=0.10, units='m')\n", - "\n", - "traj = p.model.add_subsystem('traj', dm.Trajectory())\n", - "\n", - "transcription = dm.Radau(num_segments=25, order=3, compressed=False)\n", - "phase = dm.Phase(ode_class=CannonballODE, transcription=transcription)\n", - "\n", - "phase = traj.add_phase('phase', phase)\n", - "\n", - "# Duration is bounded to be greater than 1 to ensure the nonlinear solve\n", - "# does not converge to the initial point.\n", - "phase.set_time_options(fix_initial=True, duration_bounds=(1, 100), units='s')\n", - "\n", - "# All initial states except flight path angle are fixed\n", - "# The output of the ODE which provides the rate source for each state\n", - "# is obtained from the tags used on those outputs in the ODE.\n", - "# 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, 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", - "phase.add_parameter('S', units='m**2', static_target=True)\n", - "phase.add_parameter('m', units='kg', static_target=True)\n", - "phase.add_parameter('CD', val=0.5, opt=False, static_target=True)\n", - "\n", - "phase.add_boundary_constraint('ke', loc='initial',\n", - " upper=400000, lower=0, ref=100000)\n", - "\n", - "# A duration balance is added setting altitude to zero.\n", - "# A nonlinear solver is used to find the duration of required to satisfy the condition.\n", - "phase.set_duration_balance('h', val=0.0)\n", - "\n", - "# Prior to setup, this phase has the default NonlinearRunOnce and LinearRunOnce solvers.\n", - "# Because of the implicit behavior introduced by set_duration_balance, it will automatically\n", - "# use a NewtonSolver with an Armijo Goldstein linesearch.\n", - "# In this case that linesearch causes extrapolation of the atmospheric model into invalid regions.\n", - "# To account for this, we explicitly assign a NewtonSolver with a different linesearch.\n", - "phase.nonlinear_solver = om.NewtonSolver(iprint=0, solve_subsystems=True)\n", - "phase.nonlinear_solver.linesearch = None\n", - "\n", - "phase.add_objective('r', loc='final', scaler=-1.0)\n", - "\n", - "p.model.connect('size_comp.mass', 'traj.phase.parameters:m')\n", - "p.model.connect('size_comp.S', 'traj.phase.parameters:S')\n", - "\n", - "# Finish Problem Setup\n", - "p.setup()\n", - "\n", - "#############################################\n", - "# Set constants and initial guesses\n", - "#############################################\n", - "p.set_val('radius', 0.05, units='m')\n", - "p.set_val('dens', 7.87, units='g/cm**3')\n", - "\n", - "p.set_val('traj.phase.parameters:CD', 0.5)\n", - "\n", - "guess = initial_guess(t_dur=30.0)\n", - "\n", - "phase.set_time_val(initial=0.0, duration=30.0)\n", - "\n", - "phase.set_state_val('r', vals=guess['r'], time_vals=guess['t'])\n", - "phase.set_state_val('h', vals=guess['h'], time_vals=guess['t'])\n", - "phase.set_state_val('v', vals=guess['v'], time_vals=guess['t'])\n", - "phase.set_state_val('gam', vals=guess['gam'], time_vals=guess['t'], units='rad')\n", - "\n", - "#####################################################\n", - "# Run the optimization and final explicit simulation\n", - "#####################################################\n", - "dm.run_problem(p, simulate=True, make_plots=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Displaying the timeseries data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "from IPython.display import HTML\n", - "\n", - "# Define the path to the HTML file\n", - "html_file_path = p.get_reports_dir() / 'traj_results_report.html'\n", - "html_content = html_file_path.read_text()\n", - "\n", - "# Inject CSS to control the output cell height and avoid scrollbars\n", - "html_with_custom_height = f\"\"\"\n", - "
\n", - " {html_content}\n", - "
\n", - "\"\"\"\n", - "\n", - "HTML(html_with_custom_height)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "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.phase.timeseries.r')[-1],\n", - " 3183.25, tolerance=1.0E-2)" - ] - } - ], - "metadata": { - "celltoolbar": "Tags", - "jupytext": { - "cell_metadata_filter": "-all", - "notebook_metadata_filter": "-all", - "text_representation": { - "extension": ".md", - "format_name": "markdown" - } - }, - "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.4" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/dymos_book/examples/implicit_duration/brachistochrone_implicit_duration.ipynb b/docs/dymos_book/examples/implicit_duration/brachistochrone_implicit_duration.ipynb new file mode 100644 index 000000000..0814c8555 --- /dev/null +++ b/docs/dymos_book/examples/implicit_duration/brachistochrone_implicit_duration.ipynb @@ -0,0 +1,497 @@ +{ + "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 # noqa: F401\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 # noqa: F401\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 # noqa: F401\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": [ + "# Implicit Duration\n", + "\n", + "In some cases, using dymos without the need for an optimizer is desirable.\n", + "Sometimes in the absense of a control, you may just want to propagate the problem.\n", + "Solve segments allows us to do this for a specified amount of time, but if the duration of the phase is itself an implicit output, then we need dymos to understand the residual to associate with the phase duration.\n", + "\n", + "The requirement for this procedure to work is that the dynamics have to be obeyed at each iteration.\n", + "A time-stepping transcription (`ExplicitShooting`) or a pseudospectral transcription with `solve_segments=True` is required.\n", + "\n", + "One important note: In the pseudospectral transcriptions, using `solve_segments` with option `compressed=False`, results in a \"multiple shooting\" approach, whereby the state discontinuitites between segments are treated as constraints to be nulled out by the optimizer. This approach should therefore only be used with \"single shooting\" approaches (`compressed=True` in the case of the pseudospectral transcriptions) \n", + "\n", + "## Solving the brachistochrone without an optimizer.\n", + "\n", + "Consider the brachistochrone problem. Below we implement the brachistochrone equations of motion using OpenMDAO's `JaxExplicitComponent` so that we can use automatic differentiation to provide the partials.\n", + "\n", + "The implicit solution to the brachistochrone is very robust, but what if we wanted to simplify it?\n", + "For instance, we leave the duration of the bead's trajectory as a design variable that we want to minimize, but we also require the enforcement of the final position of the brachistochrone to some (x, y) coordinate pair." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import jax.numpy as jnp\n", + "import openmdao.api as om\n", + "\n", + "\n", + "class BrachistochroneODE(om.JaxExplicitComponent):\n", + " \"\"\"\n", + " The brachistochrone EOM assuming \n", + " \"\"\"\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", + " # static parameters\n", + " self.add_input('v', desc='speed of the bead on the wire', shape=(nn,), units='m/s')\n", + " self.add_input('theta', desc='angle between wire tangent and the nadir', shape=(nn,), units='rad')\n", + " self.add_input('g', desc='gravitational acceleration', tags=['dymos.static_target'], units='m/s**2')\n", + "\n", + " self.add_output('xdot', desc='velocity component in x', shape=(nn,), units='m/s',\n", + " tags=['dymos.state_rate_source:x', 'dymos.state_units:m'])\n", + "\n", + " self.add_output('ydot', desc='velocity component in y', shape=(nn,), units='m/s',\n", + " tags=['dymos.state_rate_source:y', 'dymos.state_units:m'])\n", + "\n", + " self.add_output('vdot', desc='acceleration magnitude', shape=(nn,), units='m/s**2',\n", + " tags=['dymos.state_rate_source:v', 'dymos.state_units:m/s'])\n", + "\n", + " # Setup partials\n", + " ar = jnp.arange(nn, dtype=int)\n", + "\n", + " self.declare_partials(of='vdot', wrt='v', rows=ar, cols=ar)\n", + " self.declare_partials(of='vdot', wrt='theta', rows=ar, cols=ar)\n", + " self.declare_partials(of='vdot', wrt='g')\n", + "\n", + " self.declare_partials(of='xdot', wrt='v', rows=ar, cols=ar)\n", + " self.declare_partials(of='xdot', wrt='theta', rows=ar, cols=ar)\n", + "\n", + " self.declare_partials(of='ydot', wrt='v', rows=ar, cols=ar)\n", + " self.declare_partials(of='ydot', wrt='theta', rows=ar, cols=ar)\n", + "\n", + " def compute(self, v, theta, g):\n", + " sin_theta = jnp.sin(theta)\n", + " cos_theta = jnp.cos(theta)\n", + "\n", + " vdot = g * cos_theta\n", + " xdot = v * sin_theta\n", + " ydot = -v * cos_theta\n", + "\n", + " return xdot, ydot, vdot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Propagating the equations of motion using solve_segments\n", + "\n", + "In the following example, we set `solve_segments='forward'` such that all of the collocation defects associate with the Radau transcription are satisfied by a Newton solver rather than the optimizer.\n", + "This is the first step to solving the problem without an optimizer.\n", + "This ensures that the equations of motion are satisfied to the extent possible assuming the transcription grid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the OpenMDAO problem\n", + "p = om.Problem(model=om.Group())\n", + "\n", + "# Define a Trajectory object\n", + "traj = dm.Trajectory()\n", + "p.model.add_subsystem('traj', subsys=traj)\n", + "\n", + "# Define a Dymos Phase object with Radau Transcription\n", + "phase = dm.Phase(ode_class=BrachistochroneODE,\n", + " transcription=dm.Radau(num_segments=5, order=3, solve_segments='forward'))\n", + "traj.add_phase(name='phase0', phase=phase)\n", + "\n", + "# Set the time options\n", + "phase.set_time_options(fix_initial=True,\n", + " duration_bounds=(0.5, 10.0))\n", + "\n", + "# Set the state options\n", + "phase.add_state('x', rate_source='xdot',\n", + " fix_initial=True)\n", + "phase.add_state('y', rate_source='ydot',\n", + " fix_initial=True)\n", + "phase.add_state('v', rate_source='vdot',\n", + " fix_initial=True)\n", + "\n", + "phase.add_control('theta', lower=0.0, upper=1.5, units='rad')\n", + "phase.add_parameter('g', opt=False, val=9.80665, units='m/s**2')\n", + "\n", + "# Minimize final time.\n", + "phase.add_objective('time', loc='final')\n", + "\n", + "# Set the driver.\n", + "p.driver = om.ScipyOptimizeDriver()\n", + "\n", + "# Allow OpenMDAO to automatically determine our sparsity pattern.\n", + "# Doing so can significant speed up the execution of dymos.\n", + "p.driver.declare_coloring()\n", + "\n", + "phase.nonlinear_solver = om.NewtonSolver(maxiter=1000, solve_subsystems=True, stall_limit=3, iprint=2)\n", + "phase.linear_solver = om.DirectSolver()\n", + "\n", + "# Setup the problem\n", + "p.setup()\n", + "\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", + "phase.set_control_val('theta', vals=[0.0, 90.0], units='deg')\n", + "phase.set_parameter_val('g', val=9.80665)\n", + "\n", + "# Run the driver to solve the problem\n", + "p.run_model()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we plot the resulting trajectory of the bead, we notice that our guesses for time and the control history didn't bring the bead to the desired target at (10, 5):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "\n", + "x = p.get_val('traj.phase0.timeseries.x')\n", + "y = p.get_val('traj.phase0.timeseries.y')\n", + "\n", + "plt.plot(0.0, 10.0, 'ko')\n", + "plt.plot(10.0, 5.0, 'ko')\n", + "plt.plot(x, y)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Stopping the propagation at the desired time.\n", + "\n", + "We can utilize the `set_implicit_duration` method on `Phase` to turn `t_duration` into an implicit output and provide it with a residual.\n", + "\n", + "In our case, we set `phase.set_duration_balance('x', 10.0)`\n", + "\n", + "This specifies that we want `t_duration` to be associated with the following residual:\n", + "\n", + "\\begin{align}\n", + " \\mathcal{R}(t_d) = x - 10 = 0\n", + "\\end{align}\n", + "\n", + "Note that we limit the angle $\\theta$ to be between 0 and 180, the bead must move to the right and we can just terminate the propagation when x=10." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the OpenMDAO problem\n", + "p = om.Problem(model=om.Group())\n", + "\n", + "# Define a Trajectory object\n", + "traj = dm.Trajectory()\n", + "p.model.add_subsystem('traj', subsys=traj)\n", + "\n", + "# Define a Dymos Phase object with Radau Transcription\n", + "phase = dm.Phase(ode_class=BrachistochroneODE,\n", + " transcription=dm.Radau(num_segments=5, order=3, solve_segments='forward', compressed=True))\n", + "traj.add_phase(name='phase0', phase=phase)\n", + "\n", + "# Set the time options\n", + "phase.set_time_options(fix_initial=True,\n", + " duration_bounds=(0.5, 10.0))\n", + "phase.set_duration_balance('x', 10.0)\n", + "\n", + "# Set the state options\n", + "phase.add_state('x', rate_source='xdot',\n", + " fix_initial=True)\n", + "phase.add_state('y', rate_source='ydot',\n", + " fix_initial=True)\n", + "phase.add_state('v', rate_source='vdot',\n", + " fix_initial=True)\n", + "\n", + "phase.add_control('theta', lower=0.0, upper=1.5, units='rad')\n", + "phase.add_parameter('g', opt=False, val=9.80665, units='m/s**2')\n", + "\n", + "# Minimize final time.\n", + "phase.add_objective('time', loc='final')\n", + "\n", + "# Set the driver.\n", + "p.driver = om.ScipyOptimizeDriver()\n", + "\n", + "# Allow OpenMDAO to automatically determine our sparsity pattern.\n", + "# Doing so can significant speed up the execution of dymos.\n", + "p.driver.declare_coloring()\n", + "\n", + "phase.nonlinear_solver = om.NewtonSolver(maxiter=1000, solve_subsystems=True, stall_limit=3)\n", + "phase.nonlinear_solver.linesearch = None\n", + "phase.linear_solver = om.DirectSolver()\n", + "\n", + "# Setup the problem\n", + "p.setup()\n", + "\n", + "phase.set_time_val(initial=0.0, duration=1.8)\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", + "phase.set_control_val('theta', vals=[0.0, 90.0], units='deg')\n", + "phase.set_parameter_val('g', val=9.80665)\n", + "\n", + "# Run the driver to solve the problem\n", + "p.run_model()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "\n", + "x = p.get_val('traj.phase0.timeseries.x')\n", + "y = p.get_val('traj.phase0.timeseries.y')\n", + "v = p.get_val('traj.phase0.timeseries.y')\n", + "theta = p.get_val('traj.phase0.timeseries.theta')\n", + "\n", + "plt.plot(0.0, 10.0, 'ko')\n", + "plt.plot(10.0, 5.0, 'ko')\n", + "plt.plot(x, y)\n", + "plt.show()\n", + "\n", + "print('final time', p.get_val('traj.phase0.timeseries.time')[-1])\n", + "print('final y', y[-1])\n", + "print('final theta', jnp.degrees(theta[-1]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving the brachistochrone without an optimizer\n", + "\n", + "The brachistochrone has an analytic solution such that rate of change of the tangent to the wire (angle $\\theta) is constant.\n", + "We can use this information to pose the brachistochrone as a boundary value problem and solve as a shooting problem.\n", + "\n", + "We add another implicit output that provides the rate of change as a constant value throughout the trajectory.\n", + "That means that `theta_rate` is a parameter as far as dymos is concerned.\n", + "We then make `theta` a state variable whose rate is provided by the `theta_rate` parameter.\n", + "We'll set the initial value of `theta` to zero degrees, since we know that is approximately correct for this problem.\n", + "\n", + "\\begin{align}\n", + " \\mathcal{R}(\\dot{\\theta}) &= y_f - 5\n", + "\\end{align}\n", + "\n", + "Dymos currently doesn't support making parameters implicit outputs, so we're going to do this the manual OpenMDAO way." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the OpenMDAO problem\n", + "p = om.Problem(model=om.Group())\n", + "\n", + "# Define a Trajectory object\n", + "traj = dm.Trajectory()\n", + "p.model.add_subsystem('traj', subsys=traj)\n", + "\n", + "# Define a Dymos Phase object with Radau Transcription\n", + "phase = dm.Phase(ode_class=BrachistochroneODE,\n", + " transcription=dm.Radau(num_segments=10, order=3, solve_segments='forward', compressed=True))\n", + "traj.add_phase(name='phase0', phase=phase)\n", + "\n", + "# Set the time options\n", + "phase.set_time_options(fix_initial=True)\n", + "phase.set_duration_balance('x', 10.0)\n", + "\n", + "# Set the state options\n", + "phase.add_state('x', rate_source='xdot',\n", + " fix_initial=True)\n", + "phase.add_state('y', rate_source='ydot',\n", + " fix_initial=True)\n", + "phase.add_state('v', rate_source='vdot',\n", + " fix_initial=True)\n", + "phase.add_state('theta', rate_source='theta_rate',\n", + " fix_initial=True, units='rad')\n", + "\n", + "phase.add_parameter('theta_rate', opt=False, units='rad/s')\n", + "phase.add_parameter('g', opt=False, val=9.80665, units='m/s**2')\n", + "\n", + "# Minimize final time.\n", + "phase.add_objective('time', loc='final')\n", + "\n", + "# Set the driver.\n", + "p.driver = om.ScipyOptimizeDriver()\n", + "\n", + "# Allow OpenMDAO to automatically determine our sparsity pattern.\n", + "# Doing so can significant speed up the execution of dymos.\n", + "p.driver.declare_coloring()\n", + "\n", + "phase.nonlinear_solver = om.NewtonSolver(maxiter=1000, solve_subsystems=True, stall_limit=3,\n", + " iprint=0, atol=1.0E-8, rtol=1.0E-8, debug_print=True)\n", + "phase.nonlinear_solver.linesearch = None\n", + "phase.linear_solver = om.DirectSolver()\n", + "\n", + "theta_rate_balance = p.model.add_subsystem('theta_rate_balance', om.BalanceComp())\n", + "\n", + "theta_rate_balance.add_balance('theta_rate', units='rad/s', eq_units='m', lhs_name='y_final', rhs_val=5.0)\n", + "p.model.connect('theta_rate_balance.theta_rate', 'traj.phase0.parameters:theta_rate')\n", + "p.model.connect('traj.phase0.timeseries.y', 'theta_rate_balance.y_final', src_indices=[-1])\n", + "\n", + "# Now we need a solver to converge the loop around the entire problem.\n", + "p.model.nonlinear_solver = om.NewtonSolver(maxiter=100, solve_subsystems=True, stall_limit=3,\n", + " iprint=0, atol=1.0E-8, rtol=1.0E-8, debug_print=True)\n", + "p.model.nonlinear_solver.linesearch = None\n", + "p.model.linear_solver = om.DirectSolver()\n", + "\n", + "# Setup the problem\n", + "p.setup()\n", + "\n", + "phase.set_time_val(initial=0.0, duration=1.8)\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", + "phase.set_state_val('theta', vals=[0.0, 45], units='deg')\n", + "phase.set_parameter_val('g', val=9.80665)\n", + "\n", + "# Note that we set theta_rate at the balance comp, which is then passed into the phase as a parameter.\n", + "p.model.set_val('theta_rate_balance.theta_rate', val=0.5, units='rad/s')\n", + "\n", + "# Run the driver to solve the problem\n", + "p.run_model()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "\n", + "x = p.get_val('traj.phase0.timeseries.x')\n", + "y = p.get_val('traj.phase0.timeseries.y')\n", + "v = p.get_val('traj.phase0.timeseries.y')\n", + "theta = p.get_val('traj.phase0.timeseries.theta')\n", + "\n", + "plt.plot(0.0, 10.0, 'ko')\n", + "plt.plot(10.0, 5.0, 'ko')\n", + "plt.plot(x, y)\n", + "plt.show()\n", + "\n", + "print('final time', p.get_val('traj.phase0.timeseries.time')[-1])\n", + "print('final y', y[-1])\n", + "print('final theta', jnp.degrees(theta[-1]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using a combination of solve_segments, implicit phase duration, and an outer loop to make a phase parameter an implicit output, we've managed to solve a two point boundary value problem without the need of an optimizer. Since we're solving the segments in an inner loop here, each converged phase provides a physically valid traejctory, though it may not satisfy our desired final conditions unless the solvers are converged at all levels.\n", + "\n", + "Note that solving problems in this way can be filled with pitfalls.\n", + "Here we have a formulation of the brachistochrone that is robust to a variety of values for the parameter `theta_rate`.\n", + "An alternative formulation is to keep the ratio $\\frac{\\sin{\\theta}}{v}$ constant. However, using the inverse sine in that formulation can easily result in domain issues which makes convergence difficult.\n", + "\n", + "Even in this case, setting the `solve_subsystems` option on the outermost Newton solver to `False` seems to be too challenging and the solver fails to converge." + ] + } + ], + "metadata": { + "celltoolbar": "Tags", + "jupytext": { + "cell_metadata_filter": "-all", + "notebook_metadata_filter": "-all", + "text_representation": { + "extension": ".md", + "format_name": "markdown" + } + }, + "kernelspec": { + "display_name": "py311_2", + "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.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/dymos/phase/test/test_add_boundary_constraint.py b/dymos/phase/test/test_add_boundary_constraint.py index 732c7420f..4e2b003e2 100644 --- a/dymos/phase/test/test_add_boundary_constraint.py +++ b/dymos/phase/test/test_add_boundary_constraint.py @@ -5,6 +5,7 @@ import openmdao.api as om from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse from openmdao.utils.assert_utils import assert_near_equal +from dymos.utils.misc import om_version import dymos as dm @@ -30,7 +31,6 @@ def setup(self): self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange) self.declare_partials(of='vdot', wrt='theta', rows=arange, cols=arange) - self.declare_partials('*', '*', method='cs') self.declare_coloring(wrt='*', method='cs', show_sparsity=False) def compute(self, inputs, outputs): @@ -50,6 +50,8 @@ def compute(self, inputs, outputs): @use_tempdirs class TestAddBoundaryConstraint(unittest.TestCase): + @unittest.skipIf(om_version()[0] >= (3, 36, 1) and om_version()[0] < (3, 37, 0), + reason='Test is broken in OpenMDAO 3.37.0 interim development.') @require_pyoptsparse(optimizer='SLSQP') def test_simple_no_exception(self): p = om.Problem(model=om.Group()) diff --git a/dymos/phase/test/test_implicit_duration.py b/dymos/phase/test/test_implicit_duration.py index b3ea1e392..ba244c100 100644 --- a/dymos/phase/test/test_implicit_duration.py +++ b/dymos/phase/test/test_implicit_duration.py @@ -203,3 +203,7 @@ def test_implicit_duration_gl_expr_condition(self): assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.4735192, tolerance=1E-6) assert_near_equal(p.get_val('traj.phase.timeseries.h')[-1], 0.0, tolerance=1E-6) + + +if __name__ == '__main__': + unittest.main() diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py index e360fa5eb..3df618fc7 100644 --- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py +++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py @@ -466,6 +466,8 @@ def configure_duration_balance(self, phase): duration_balance_comp = phase._get_subsystem('t_duration_balance_comp') configure_duration_balance_introspection(phase) options = phase.time_options['t_duration_balance_options'] + lower, upper = phase.time_options['duration_bounds'] + duration_val = phase.time_options['duration_val'] if options['mult_val'] is None: use_mult = False @@ -476,9 +478,9 @@ def configure_duration_balance(self, phase): src_idx = [-1] if options['index'] is None else [[-1]]+options['index'] - duration_balance_comp.add_balance('t_duration', val=options['val'], + duration_balance_comp.add_balance('t_duration', val=duration_val, lower=lower, upper=upper, eq_units=options['units'], units=phase.time_options['units'], - use_mult=use_mult, mult_val=mult_val) + rhs_val=options['val'], use_mult=use_mult, mult_val=mult_val) phase.connect('t_duration_balance_comp.t_duration', 't_duration')