diff --git a/docs/dymos_book/_toc.yml b/docs/dymos_book/_toc.yml index 9d3147daf..c1c2fb882 100644 --- 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 index 1a3b1569a..4bb2c42fe 100644 --- 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 new file mode 100644 index 000000000..2b69b3136 --- /dev/null +++ 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 +} diff --git a/docs/dymos_book/getting_started/collocation.ipynb b/docs/dymos_book/getting_started/collocation.ipynb index 422b4ab91..f27c1ab80 100644 --- a/docs/dymos_book/getting_started/collocation.ipynb +++ b/docs/dymos_book/getting_started/collocation.ipynb @@ -106,7 +106,7 @@ "\n", "_How long does it take the ball to reach the ground from an initial height of 100 m?_\n", "\n", - "We will explain how this problem is solved using the two implicit simulation techniques in Dymos: high-order Gauss-Lobatto collocation, and the Radau Pseudospectral Method.\n", + "We will explain how this problem is solved using the three implicit simulation techniques in Dymos: high-order Gauss-Lobatto collocation, the Radau Pseudospectral Method, and the Birkhoff Pseudospectral Method.\n", "\n", "## Finding the trajectory using high-order Legendre-Gauss-Lobatto collocation\n", "\n", @@ -207,11 +207,44 @@ "\n", "![Screenshot](scripts/lgr_solution_5.png)\n", "\n", + "## Finding the trajectory using the Birkhoff Pseudospectral Method\n", + "\n", + "As before, there are two states and two 3rd-order polynomials representing the time history of each states.\n", + "In the Birkhoff pseudospectral method, all states and state rates are used as free variables.\n", + "Additionally, the initial and final nodes are represented independently and may also be free variables depending on the problem.\n", + "So each state eight pieces of information to be fully solved.\n", + "In this problem, the initial and final points of the position and the initial point of the velocity are fixed.\n", + "The final velocity is left open to be determined by the terminal boundary constraint.\n", + "\n", + "![Screenshot](scripts/birkhoff_animation_1.png)\n", + "\n", + "The state variables constrained such that the end points are equal to the initial and final nodes.\n", + "For illustrative purposes, the state variables are initialized as the average of the guess for initial and final states.\n", + "\n", + "![Screenshot](scripts/birkhoff_animation_2.png)\n", + "\n", + "We can use the initial state value, the state rates, and an integration matrix to compute the state value at each subsequent node.\n", + "These values can be compared to obtain the two state defects.\n", + "\n", + "![Screenshot](scripts/birkhoff_animation_3.png)\n", + "\n", + "The state values can be used to evaluate the ODE.\n", + "The state rates computed from the ODE can be compared against the state rate free variables to obtain the three state rate defects.\n", + "\n", + "![Screenshot](scripts/birkhoff_animation_4.png)\n", + "\n", + "\n", + "The free variables are iterated until the defects are zero and the boundary constraints are satisfied.\n", + "\n", + "![Screenshot](scripts/birkhoff_animation_5.png)\n", + "\n", + "\n", "## Multiple Segments\n", "\n", "The simple examples above use one single 3rd-order polynomial to represent the time history of the states.\n", "With more complex states, a single 3rd-order polynomial would not be capable of replicating the time history.\n", - "In practice, we can use more polynomial segments _and_ higher-order polynomials to better match the dynamics.\n", + "In practice, we can use more polynomial segments _and_ higher-order polynomials to better match the dynamics for the Gauss-Lobatto and Radau Pseudospectral methods.\n", + "For the Birkhoff pseudospectral method we only use one segment but use higher-order polynomials.\n", "\n", "In the case of adding more segments, one can choose whether the segment boundaries are shared by the segments, or independent.\n", "If they are shared, fewer design variables are necessary, but it may be more difficult for the optimizer to search for a solution.\n", @@ -227,6 +260,16 @@ "In addition, the interpolation step used by the Gauss-Lobatto method can result in interpolated state values falling well outside the user's expected range if the initial guess is not sufficiently accurate.\n", "This can cause convergence issues if there are nonlinear solvers within the ODE that rely on a reasonable guess to achieve convergence.\n", "\n", + "## Birkhoff Pseudospectral Method vs LGL and LGR Collocation\n", + "\n", + "There are two notable differences between the birkhoff pseudospectral method and the two other methods.\n", + "Firstly, the collocation constraints are applied in _integral_ form rather than _differential_ form.\n", + "Secondly, an additional layer of slack is introduced to the NLP by introducing _virtual variables_ as stand-ins for the state rates.\n", + "These virtual variables are constrained to equal the state rate obtained by evaluating the ODE and are used to compute state values and apply the state defects.\n", + "These two factors combine to produce an NLP that is better conditioned and whose solution is more equivalent to the solution of the original optimal control problem.\n", + "The improved conditioning allows for the use of very large segments in the discretization of the state and therefore only one large segment is utilized when using this transcription.\n", + "The stronger optimality conditions come with the cost of additional iterations of the NLP solver but the resulting solutions typically have lower cost functions at the expense of additional compute time.\n", + "\n", "## Satisfying Defects with a Nonlinear Solver\n", "\n", "Typically, collocation problems are solved by posing the defects as constraints to an optimizer.\n", @@ -271,7 +314,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.11.0" } }, "nbformat": 4, diff --git a/docs/dymos_book/getting_started/scripts/birkhoff_animation_1.png b/docs/dymos_book/getting_started/scripts/birkhoff_animation_1.png new file mode 100644 index 000000000..05e0f58c0 Binary files /dev/null and b/docs/dymos_book/getting_started/scripts/birkhoff_animation_1.png differ diff --git a/docs/dymos_book/getting_started/scripts/birkhoff_animation_2.png b/docs/dymos_book/getting_started/scripts/birkhoff_animation_2.png new file mode 100644 index 000000000..c4bc159ff Binary files /dev/null and b/docs/dymos_book/getting_started/scripts/birkhoff_animation_2.png differ diff --git a/docs/dymos_book/getting_started/scripts/birkhoff_animation_3.png b/docs/dymos_book/getting_started/scripts/birkhoff_animation_3.png new file mode 100644 index 000000000..b105f6e94 Binary files /dev/null and b/docs/dymos_book/getting_started/scripts/birkhoff_animation_3.png differ diff --git a/docs/dymos_book/getting_started/scripts/birkhoff_animation_4.png b/docs/dymos_book/getting_started/scripts/birkhoff_animation_4.png new file mode 100644 index 000000000..1f34094e6 Binary files /dev/null and b/docs/dymos_book/getting_started/scripts/birkhoff_animation_4.png differ diff --git a/docs/dymos_book/getting_started/scripts/birkhoff_animation_5.png b/docs/dymos_book/getting_started/scripts/birkhoff_animation_5.png new file mode 100644 index 000000000..3fadf7b3a Binary files /dev/null and b/docs/dymos_book/getting_started/scripts/birkhoff_animation_5.png differ diff --git a/docs/dymos_book/references.bib b/docs/dymos_book/references.bib index 381796caa..09cbe4d2e 100644 --- a/docs/dymos_book/references.bib +++ b/docs/dymos_book/references.bib @@ -595,3 +595,15 @@ @book{hull2003oct year={2003}, publisher={Springer} } + +@article{Meditch1964, + author={Meditch, J.}, + journal={IEEE Transactions on Automatic Control}, + title={On the problem of optimal thrust programming for a lunar soft landing}, + year={1964}, + volume={9}, + number={4}, + pages={477-484}, + keywords={Moon;Optimal control;Control systems;Variable speed drives;Fuels;Vehicles;History;Feedback control;Books;Kalman filters}, + doi={10.1109/TAC.1964.1105758}} + diff --git a/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py b/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py index 50d5d8f11..b78e4b609 100644 --- a/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py +++ b/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py @@ -227,9 +227,7 @@ def test_cruise_results_birkhoff(self): assumptions.add_output('mass_empty', val=1.0, units='kg') assumptions.add_output('mass_payload', val=1.0, units='kg') - grid = dm.BirkhoffGrid(num_nodes=50, grid_type='lgl') - - transcription = dm.Birkhoff(grid=grid, solve_segments='forward') + transcription = dm.Birkhoff(num_nodes=50, grid_type='lgl', solve_segments='forward') phase = dm.Phase(ode_class=AircraftODE, transcription=transcription) traj = dm.Trajectory() traj.add_phase('phase0', phase) diff --git a/dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py b/dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py index 9aa715da4..d51b4c9ca 100644 --- a/dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py +++ b/dymos/examples/aircraft_steady_flight/test/test_ex_aircraft_steady_flight.py @@ -135,7 +135,7 @@ def test_ex_aircraft_steady_flight_opt_radau(self): @unittest.skipIf(om_version < (3, 32, 2), 'Test requires OpenMDAO 3.32.2 or later') def test_ex_aircraft_steady_flight_opt_birkhoff(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=50, grid_type='cgl'), + tx = dm.Birkhoff(num_nodes=50, grid_type='cgl', solve_segments=False) p = ex_aircraft_steady_flight(transcription=tx, optimizer='IPOPT') assert_near_equal(p.get_val('traj.phase0.timeseries.range', units='NM')[-1], diff --git a/dymos/examples/brachistochrone/test/ex_brachistochrone.py b/dymos/examples/brachistochrone/test/ex_brachistochrone.py index 03ed2c850..62ff993fc 100644 --- a/dymos/examples/brachistochrone/test/ex_brachistochrone.py +++ b/dymos/examples/brachistochrone/test/ex_brachistochrone.py @@ -41,9 +41,8 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran compressed=compressed) t = dm.ExplicitShooting(grid=grid) elif transcription == 'birkhoff': - grid = dm.BirkhoffGrid(num_nodes=transcription_order + 1, - grid_type='cgl') - t = dm.Birkhoff(grid=grid) + t = dm.Birkhoff(num_nodes=transcription_order + 1, + grid_type='cgl') # phase = dm.ImplicitPhase(ode_class=BrachistochroneODE, num_nodes=11) traj = dm.Trajectory() diff --git a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py index c5ea6ba09..8840492ba 100644 --- a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py +++ b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py @@ -34,8 +34,8 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran compressed=compressed) fix_final = not solve_segments elif transcription == 'birkhoff': - gd = dm.BirkhoffGrid(num_nodes=25, grid_type=grid_type) - transcription = dm.Birkhoff(grid=gd) + transcription = dm.Birkhoff(num_nodes=25, + grid_type=grid_type) fix_final = not solve_segments traj = dm.Trajectory() diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_birkhoff_constraints.py b/dymos/examples/brachistochrone/test/test_brachistochrone_birkhoff_constraints.py index 5ed178422..222c568d5 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_birkhoff_constraints.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_birkhoff_constraints.py @@ -19,8 +19,7 @@ def test_brachistochrone_control_prefix(self): p.driver.options['optimizer'] = 'SLSQP' p.driver.declare_coloring(tol=1.0E-12) - grid = dm.BirkhoffGrid(num_nodes=25, grid_type='lgl') - tx = dm.Birkhoff(grid=grid) + tx = dm.Birkhoff(num_nodes=25, grid_type='lgl') traj = dm.Trajectory() phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx) @@ -75,8 +74,7 @@ def test_brachistochrone_control_no_prefix(self): p.driver.options['optimizer'] = 'SLSQP' p.driver.declare_coloring(tol=1.0E-12) - grid = dm.BirkhoffGrid(num_nodes=25, grid_type='lgl') - tx = dm.Birkhoff(grid=grid) + tx = dm.Birkhoff(num_nodes=25, grid_type='lgl') traj = dm.Trajectory() phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx) @@ -131,8 +129,7 @@ def test_brachistochrone_ode_prefix(self): p.driver.options['optimizer'] = 'SLSQP' p.driver.declare_coloring(tol=1.0E-12) - grid = dm.BirkhoffGrid(num_nodes=25, grid_type='lgl') - tx = dm.Birkhoff(grid=grid) + tx = dm.Birkhoff(num_nodes=25, grid_type='lgl') traj = dm.Trajectory() phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx) @@ -186,8 +183,7 @@ def test_brachistochrone_ode_no_prefix(self): p.driver.options['optimizer'] = 'SLSQP' p.driver.declare_coloring(tol=1.0E-12) - grid = dm.BirkhoffGrid(num_nodes=25, grid_type='lgl') - tx = dm.Birkhoff(grid=grid) + tx = dm.Birkhoff(num_nodes=25, grid_type='lgl') traj = dm.Trajectory() phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx) diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py b/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py index e1c388e22..595847fe8 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py @@ -31,9 +31,8 @@ def _make_problem(self, transcription='gauss-lobatto', num_segments=8, transcrip order=transcription_order, compressed=compressed) elif transcription == 'birkhoff': - grid = dm.BirkhoffGrid(num_nodes=transcription_order+1, - grid_type='cgl') - t = dm.Birkhoff(grid=grid) + t = dm.Birkhoff(num_nodes=transcription_order+1, + grid_type='cgl') def ode(num_nodes): return om.ExecComp(['vdot = g * cos(theta)', diff --git a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py index 166a95c72..d817c6a25 100644 --- a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py +++ b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py @@ -198,7 +198,7 @@ def test_brachistochrone_polynomial_control_birkhoff(self): p.driver.declare_coloring() phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=15))) + transcription=dm.Birkhoff(num_nodes=15)) p.model.add_subsystem('phase0', phase) @@ -464,7 +464,7 @@ def test_brachistochrone_polynomial_control_birkhoff(self): p.driver.declare_coloring() phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=15))) + transcription=dm.Birkhoff(num_nodes=15)) p.model.add_subsystem('phase0', phase) @@ -729,7 +729,7 @@ def test_brachistochrone_polynomial_control_birkhoff(self): p.driver.declare_coloring() phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=15))) + transcription=dm.Birkhoff(num_nodes=15)) p.model.add_subsystem('phase0', phase) @@ -997,7 +997,7 @@ def test_brachistochrone_polynomial_control_birkhoff(self): p.driver.declare_coloring() phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=15))) + transcription=dm.Birkhoff(num_nodes=15)) p.model.add_subsystem('phase0', phase) @@ -1265,7 +1265,7 @@ def test_brachistochrone_polynomial_control_birkhoff(self): p.driver.declare_coloring() phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=15))) + transcription=dm.Birkhoff(num_nodes=15)) p.model.add_subsystem('phase0', phase) @@ -1469,7 +1469,7 @@ def test_brachistochrone_polynomial_control_birkhoff(self): p.driver.declare_coloring() phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=15))) + transcription=dm.Birkhoff(num_nodes=15)) p.model.add_subsystem('phase0', phase) diff --git a/dymos/examples/brachistochrone/test/test_ex_brachistochrone_refine_grid.py b/dymos/examples/brachistochrone/test/test_ex_brachistochrone_refine_grid.py index 137e6deb1..1a3d01211 100644 --- a/dymos/examples/brachistochrone/test/test_ex_brachistochrone_refine_grid.py +++ b/dymos/examples/brachistochrone/test/test_ex_brachistochrone_refine_grid.py @@ -31,8 +31,8 @@ def make_problem(self, transcription='radau-ps', num_segments=5, transcription_o order=transcription_order, compressed=compressed) elif transcription == 'birkhoff': - grid = dm.BirkhoffGrid(num_nodes=transcription_order+1, grid_type='cgl') - t = dm.Birkhoff(grid=grid) + t = dm.Birkhoff(num_nodes=transcription_order+1, + grid_type='cgl') traj = dm.Trajectory() phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t) diff --git a/dymos/examples/cannonball/test/test_cannonball_matrix_state.py b/dymos/examples/cannonball/test/test_cannonball_matrix_state.py index cd20cb848..df2826f57 100644 --- a/dymos/examples/cannonball/test/test_cannonball_matrix_state.py +++ b/dymos/examples/cannonball/test/test_cannonball_matrix_state.py @@ -101,7 +101,7 @@ def test_cannonball_matrix_state_gl(self): @require_pyoptsparse(optimizer='IPOPT') def test_cannonball_matrix_state_birkhoff_lgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=20, grid_type='lgl')) + tx = dm.Birkhoff(num_nodes=20, grid_type='lgl') p = self._make_problem(tx) @@ -119,7 +119,7 @@ def test_cannonball_matrix_state_birkhoff_lgl(self): @require_pyoptsparse(optimizer='IPOPT') def test_cannonball_matrix_state_birkhoff_cgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=20, grid_type='cgl')) + tx = dm.Birkhoff(num_nodes=20, grid_type='cgl') p = self._make_problem(tx) diff --git a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py index cd2b3097d..1f8e3b87a 100644 --- a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py +++ b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py @@ -148,7 +148,7 @@ def make_problem(self, connected=False): traj = dm.Trajectory() p.model.add_subsystem('traj', traj) - transcription = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=10)) + transcription = dm.Birkhoff(num_nodes=10) ascent = dm.Phase(ode_class=CannonballODE, transcription=transcription) ascent = traj.add_phase('ascent', ascent) @@ -175,7 +175,7 @@ def make_problem(self, connected=False): upper=400000, lower=0, ref=100000) # Second Phase (descent) - transcription = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=10)) + transcription = dm.Birkhoff(num_nodes=10) descent = dm.Phase(ode_class=CannonballODE, transcription=transcription) traj.add_phase('descent', descent) diff --git a/dymos/examples/double_integrator/test/test_double_integrator.py b/dymos/examples/double_integrator/test/test_double_integrator.py index 02b2c5e97..c3071a244 100644 --- a/dymos/examples/double_integrator/test/test_double_integrator.py +++ b/dymos/examples/double_integrator/test/test_double_integrator.py @@ -29,7 +29,7 @@ def double_integrator(transcription='gauss-lobatto', compressed=True, grid_type= elif transcription == "radau-ps": t = dm.Radau(num_segments=30, order=3, compressed=compressed) elif transcription == 'birkhoff': - t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=100, grid_type=grid_type)) + t = dm.Birkhoff(num_nodes=100, grid_type=grid_type) else: raise ValueError('invalid transcription') diff --git a/dymos/examples/double_integrator/test/test_double_integrator_breakwell.py b/dymos/examples/double_integrator/test/test_double_integrator_breakwell.py index c314b3fb3..6b1e2755c 100644 --- a/dymos/examples/double_integrator/test/test_double_integrator_breakwell.py +++ b/dymos/examples/double_integrator/test/test_double_integrator_breakwell.py @@ -37,8 +37,7 @@ def double_integrator_direct_collocation(transcription='gauss-lobatto', compress elif transcription == "radau-ps": t = dm.Radau(num_segments=30, order=3, compressed=compressed) elif transcription == 'birkhoff': - grid = dm.BirkhoffGrid(num_nodes=51) - t = dm.Birkhoff(grid=grid) + t = dm.Birkhoff(num_nodes=51) else: raise ValueError('invalid transcription') diff --git a/dymos/examples/flying_robot/test/test_flying_robot.py b/dymos/examples/flying_robot/test/test_flying_robot.py index 985dbc176..4cf6e8844 100644 --- a/dymos/examples/flying_robot/test/test_flying_robot.py +++ b/dymos/examples/flying_robot/test/test_flying_robot.py @@ -23,8 +23,7 @@ def flying_robot_direct_collocation(transcription='gauss-lobatto', compressed=Tr elif transcription == "radau-ps": t = dm.Radau(num_segments=8, order=5, compressed=compressed) elif transcription == 'birkhoff': - grid = dm.BirkhoffGrid(num_nodes=30) - t = dm.Birkhoff(grid=grid) + t = dm.Birkhoff(num_nodes=30) else: raise ValueError('invalid transcription') diff --git a/dymos/examples/goddard_rocket_problem/test/test_goddard_rocket_problem.py b/dymos/examples/goddard_rocket_problem/test/test_goddard_rocket_problem.py index fa1538898..8d3f3a539 100644 --- a/dymos/examples/goddard_rocket_problem/test/test_goddard_rocket_problem.py +++ b/dymos/examples/goddard_rocket_problem/test/test_goddard_rocket_problem.py @@ -29,7 +29,7 @@ def goddard_rocket_direct_collocation(grid_type='lgl'): p.driver.opt_settings['Major optimality tolerance'] = 1.0E-3 p.driver.opt_settings['iSumm'] = 6 - t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=150, grid_type=grid_type)) + t = dm.Birkhoff(num_nodes=150, grid_type=grid_type) traj = p.model.add_subsystem('traj', dm.Trajectory()) diff --git a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py index b6b960478..8bc437f2d 100644 --- a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py +++ b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py @@ -56,7 +56,7 @@ def make_problem(self, transcription=dm.GaussLobatto, optimizer='SLSQP', numseg= elif transcription == 'radau-ps': t = dm.Radau(num_segments=numseg, order=3) elif transcription == 'birkhoff': - t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=order + 1, grid_type='lgl'), + t = dm.Birkhoff(num_nodes=order + 1, grid_type='lgl', solve_segments=solve_segments) traj = p.model.add_subsystem('traj', dm.Trajectory()) diff --git a/dymos/examples/low_thrust_spiral/test/test_low_thrust_spiral.py b/dymos/examples/low_thrust_spiral/test/test_low_thrust_spiral.py index 57ecf3d88..f63f8b443 100644 --- a/dymos/examples/low_thrust_spiral/test/test_low_thrust_spiral.py +++ b/dymos/examples/low_thrust_spiral/test/test_low_thrust_spiral.py @@ -23,7 +23,7 @@ def low_thrust_spiral_direct_collocation(grid_type='lgl'): p.driver.opt_settings['Major optimality tolerance'] = 5.0E-4 p.driver.opt_settings['iSumm'] = 6 - t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=101, grid_type=grid_type)) + t = dm.Birkhoff(num_nodes=101, grid_type=grid_type) traj = p.model.add_subsystem('traj', dm.Trajectory()) diff --git a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py index bb6edfbdb..6fd096779 100644 --- a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py +++ b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py @@ -51,7 +51,7 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', elif transcription == 'radau-ps': tx = dm.Radau(num_segments=num_seg, order=transcription_order) elif transcription == 'birkhoff': - tx = dm.Birkhoff(order=transcription_order) + tx = dm.Birkhoff(num_nodes=transcription_order) traj = dm.Trajectory() diff --git a/dymos/examples/moon_landing/test/test_moon_landing_problem.py b/dymos/examples/moon_landing/test/test_moon_landing_problem.py index 0ef628c3c..e3ab11092 100644 --- a/dymos/examples/moon_landing/test/test_moon_landing_problem.py +++ b/dymos/examples/moon_landing/test/test_moon_landing_problem.py @@ -22,7 +22,7 @@ def make_problem(self, grid_type): self.p.driver.opt_settings['linear_solver'] = 'mumps' self.p.driver.declare_coloring() - t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=20, grid_type=grid_type)) + t = dm.Birkhoff(num_nodes=20, grid_type=grid_type) traj = self.p.model.add_subsystem('traj', dm.Trajectory()) phase = dm.Phase(ode_class=MoonLandingProblemODE, diff --git a/dymos/examples/racecar/test/test_racecar.py b/dymos/examples/racecar/test/test_racecar.py index 32540919e..88d6bf5ee 100644 --- a/dymos/examples/racecar/test/test_racecar.py +++ b/dymos/examples/racecar/test/test_racecar.py @@ -42,7 +42,7 @@ def test_racecar_for_docs(self): txs = {'radau': dm.Radau(num_segments=50, order=3, compressed=True), 'gauss-lobatto': dm.GaussLobatto(num_segments=50, order=3, compressed=True), - 'birkhoff': dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=100))} + 'birkhoff': dm.Birkhoff(num_nodes=100)} for tx_name, tx in txs.items(): with self.subTest(tx_name): diff --git a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py index b25c9d9fa..270dce8cc 100644 --- a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py +++ b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py @@ -96,7 +96,7 @@ def test_robot_arm_gl(self): @require_pyoptsparse(optimizer='IPOPT') def test_robot_arm_birkhoff_lgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=30, grid_type='lgl')) + tx = dm.Birkhoff(num_nodes=30, grid_type='lgl') p = self.make_problem(tx=tx, optimizer='IPOPT') dm.run_problem(p) @@ -105,7 +105,7 @@ def test_robot_arm_birkhoff_lgl(self): @require_pyoptsparse(optimizer='IPOPT') def test_robot_arm_birkhoff_cgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=30, grid_type='cgl')) + tx = dm.Birkhoff(num_nodes=30, grid_type='cgl') p = self.make_problem(tx=tx, optimizer='IPOPT') dm.run_problem(p) diff --git a/dymos/examples/robot_arm/test/test_robot_arm.py b/dymos/examples/robot_arm/test/test_robot_arm.py index 9c9e3d661..424e179ee 100644 --- a/dymos/examples/robot_arm/test/test_robot_arm.py +++ b/dymos/examples/robot_arm/test/test_robot_arm.py @@ -155,7 +155,7 @@ def test_robot_arm_gl(self): @require_pyoptsparse(optimizer='IPOPT') def test_robot_arm_birkhoff_lgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=30, grid_type='lgl')) + tx = dm.Birkhoff(num_nodes=30, grid_type='lgl') p = self.make_problem(tx=tx, optimizer='IPOPT') dm.run_problem(p) @@ -181,7 +181,7 @@ def test_robot_arm_birkhoff_lgl(self): @require_pyoptsparse(optimizer='IPOPT') def test_robot_arm_birkhoff_cgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=30, grid_type='cgl')) + tx = dm.Birkhoff(num_nodes=30, grid_type='cgl') p = self.make_problem(tx=tx, optimizer='IPOPT') dm.run_problem(p) diff --git a/dymos/examples/shuttle_reentry/test/test_reentry.py b/dymos/examples/shuttle_reentry/test/test_reentry.py index b79b0c9e1..05be78634 100644 --- a/dymos/examples/shuttle_reentry/test/test_reentry.py +++ b/dymos/examples/shuttle_reentry/test/test_reentry.py @@ -121,7 +121,7 @@ def test_reentry_constrained_gauss_lobatto(self): @require_pyoptsparse(optimizer='IPOPT') def test_reentry_constrained_birkhoff_lgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=30, grid_type='lgl')) + tx = dm.Birkhoff(num_nodes=30, grid_type='lgl') p = self.make_problem(constrained=True, transcription=tx, optimizer='IPOPT') p.run_driver() @@ -135,7 +135,7 @@ def test_reentry_constrained_birkhoff_lgl(self): @require_pyoptsparse(optimizer='IPOPT') def test_reentry_constrained_birkhoff_cgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=30, grid_type='cgl')) + tx = dm.Birkhoff(num_nodes=30, grid_type='cgl') p = self.make_problem(constrained=True, transcription=tx, optimizer='IPOPT') p.run_driver() @@ -177,7 +177,7 @@ def test_reentry_unconstrained_gauss_lobatto(self): @require_pyoptsparse(optimizer='IPOPT') def test_reentry_unconstrained_birkhoff_lgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=60, grid_type='lgl')) + tx = dm.Birkhoff(num_nodes=60, grid_type='lgl') p = self.make_problem(constrained=False, transcription=tx, optimizer='IPOPT') p.run_driver() @@ -191,7 +191,7 @@ def test_reentry_unconstrained_birkhoff_lgl(self): @require_pyoptsparse(optimizer='IPOPT') def test_reentry_unconstrained_birkhoff_cgl(self): - tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=60, grid_type='cgl')) + tx = dm.Birkhoff(num_nodes=60, grid_type='cgl') p = self.make_problem(constrained=False, transcription=tx, optimizer='IPOPT') p.run_driver() diff --git a/dymos/test/test_load_case.py b/dymos/test/test_load_case.py index 1f07393ad..3aa4826db 100644 --- a/dymos/test/test_load_case.py +++ b/dymos/test/test_load_case.py @@ -198,7 +198,7 @@ def test_load_case_radau_to_birkhoff(self): case = om.CaseReader('dymos_solution.db').get_case('final') # create a problem with a different transcription with a different number of variables - q = setup_problem(dm.Birkhoff(grid=dm.BirkhoffGrid(num_nodes=50))) + q = setup_problem(dm.Birkhoff(num_nodes=50)) # Fill q with junk so that we can be sure load_case worked q['phase0.t_initial'] = -88 diff --git a/dymos/transcriptions/pseudospectral/birkhoff.py b/dymos/transcriptions/pseudospectral/birkhoff.py index 9ba383c79..78706835e 100644 --- a/dymos/transcriptions/pseudospectral/birkhoff.py +++ b/dymos/transcriptions/pseudospectral/birkhoff.py @@ -34,10 +34,13 @@ def initialize(self): """ Declare transcription options. """ - self.options.declare('grid', types=(BirkhoffGrid, str), - allow_none=True, default=None, - desc='The grid distribution used to layout the control inputs and provide the default ' - 'output nodes. Use either "cgl" or "lgl".') + self.options.declare('num_nodes', types=int, default=3, + desc='The number of nodes in the grid') + + self.options.declare('grid_type', values=('cgl', 'lgl'), default='cgl', + desc='Specifies which type of grid is to be used. ' + 'Options are Chebyshev-Gauss-Lobatto ("cgl") ' + 'and Legendre-Gauss-Lobatto ("lgl")') self.options.declare(name='solve_segments', default=False, values=(False, 'forward', 'backward'), @@ -55,11 +58,8 @@ def init_grid(self): """ Setup the GridData object for the Transcription. """ - if self.options['grid'] in ('cgl', None): - self.grid_data = BirkhoffGrid(num_nodes=self.options['order'], - grid_type='cgl') - else: - self.grid_data = self.options['grid'] + self.grid_data = BirkhoffGrid(num_nodes=self.options['num_nodes'], + grid_type=self.options['grid_type']) def setup_time(self, phase): """