From f05f58cacdc97dfbf9407f815a58bb742bd9409d Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Fri, 15 Nov 2024 12:13:28 -0500 Subject: [PATCH] Added interpolation of results to arbitrary nodes, removed use of `scipy.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 | 10 +- .../cannonball_implicit_duration.ipynb | 19 +- .../examples/water_rocket/water_rocket.ipynb | 23 +- .../features/phases/timeseries.ipynb | 158 +++- .../test_doc_cannonball_implicit_duration.py | 30 +- .../finite_burn_orbit_raise_problem.py | 37 +- .../test/test_ex_two_burn_orbit_raise_mpi.py | 4 +- .../test/test_finite_burn_eom.py | 2 +- .../test/test_multi_phase_restart.py | 22 +- .../test/test_oscillator_vector_states.py | 2 +- dymos/examples/racecar/test/test_racecar.py | 3 +- dymos/examples/water_rocket/phases.py | 8 +- .../water_rocket/test/test_water_rocket.py | 16 +- .../water_rocket/water_engine_comp.py | 4 +- dymos/phase/phase.py | 123 +-- dymos/phase/test/test_interp.py | 9 + dymos/test/test_upgrade_guide.py | 777 ------------------ .../common/test/test_control_interp_comp.py | 8 +- .../test/test_explicit_shooting.py | 16 +- .../test/test_collocation_defect_sol_opt.py | 8 +- .../test/test_control_endpoint_defect_comp.py | 2 +- 21 files changed, 311 insertions(+), 970 deletions(-) delete mode 100644 dymos/test/test_upgrade_guide.py diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml index e6b374099..7540efb15 100644 --- a/.github/workflows/dymos_tests_workflow.yml +++ b/.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 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 index f34180b4c..89c99afb2 100644 --- a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb +++ b/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, diff --git a/docs/dymos_book/examples/water_rocket/water_rocket.ipynb b/docs/dymos_book/examples/water_rocket/water_rocket.ipynb index 80756dc4a..14c74f170 100644 --- a/docs/dymos_book/examples/water_rocket/water_rocket.ipynb +++ b/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, diff --git a/docs/dymos_book/features/phases/timeseries.ipynb b/docs/dymos_book/features/phases/timeseries.ipynb index 8594c1957..14886f01b 100644 --- a/docs/dymos_book/features/phases/timeseries.ipynb +++ b/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, diff --git a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py index ddb0f9b55..5c49828c9 100644 --- a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py +++ b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py @@ -25,7 +25,7 @@ class TestTwoPhaseCannonballForDocs(unittest.TestCase): @require_pyoptsparse(optimizer='IPOPT') def test_two_phase_cannonball_for_docs(self): import openmdao.api as om - from openmdao.components.interp_util.interp import InterpND + from scipy.interpolate import make_interp_spline from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -113,9 +113,7 @@ def setup(self): alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0] rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0] - self.rho_interp = InterpND(points=np.array(alt_data), - values=np.array(rho_data), - method='slinear') + self.rho_interp = make_interp_spline(alt_data, rho_data, k=1) def compute(self, inputs, outputs): @@ -128,7 +126,7 @@ def compute(self, inputs, outputs): GRAVITY = 9.80665 # m/s**2 - rho = self.rho_interp.interpolate(h) + rho = self.rho_interp(h) q = 0.5*rho*v**2 qS = q * S @@ -147,19 +145,9 @@ def compute(self, inputs, outputs): p = om.Problem(model=om.Group()) - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'IPOPT' + p.driver = om.ScipyOptimizeDriver() p.driver.declare_coloring() - p.driver.opt_settings['derivative_test'] = 'first-order' - p.driver.opt_settings['mu_strategy'] = 'monotone' - p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas' - p.driver.opt_settings['bound_mult_init_method'] = 'mu-based' - p.driver.opt_settings['mu_init'] = 0.01 - p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' - - p.set_solver_print(level=0, depth=99) - p.model.add_subsystem('size_comp', CannonballSizeComp(), promotes_inputs=['radius', 'dens']) p.model.set_input_defaults('dens', val=7.87, units='g/cm**3') @@ -168,7 +156,7 @@ def compute(self, inputs, outputs): traj = p.model.add_subsystem('traj', dm.Trajectory()) - transcription = dm.Radau(num_segments=25, order=3, compressed=False) + transcription = dm.Radau(num_segments=20, order=3, compressed=False) phase = dm.Phase(ode_class=CannonballODE, transcription=transcription) phase = traj.add_phase('phase', phase) @@ -181,7 +169,7 @@ def compute(self, inputs, outputs): phase.set_time_options(fix_initial=True, duration_bounds=(1, 100), units='s') phase.add_state('r', fix_initial=True, solve_segments='forward') phase.add_state('h', fix_initial=True, solve_segments='forward') - phase.add_state('gam', fix_initial=False, solve_segments='forward') + phase.add_state('gam', fix_initial=False, initial_bounds=(0, np.pi/2), solve_segments='forward') phase.add_state('v', fix_initial=False, solve_segments='forward') phase.add_parameter('S', units='m**2', static_target=True) @@ -200,8 +188,10 @@ def compute(self, inputs, outputs): # In this problem, the default ArmijoGoldsteinLS has issues with extrapolating # the states and causes the optimization to fail. # Using the default linesearch or BoundsEnforceLS work better here. - phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, iprint=2) - phase.nonlinear_solver.linesearch = None + phase.nonlinear_solver = om.NewtonSolver(iprint=0, solve_subsystems=True, + maxiter=100, stall_limit=3) + phase.nonlinear_solver.linesearch = om.BoundsEnforceLS() + phase.linear_solver = om.DirectSolver() phase.add_objective('r', loc='final', scaler=-1.0) diff --git a/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py b/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py index a94e4e792..c0d67ac31 100644 --- a/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py +++ b/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py @@ -51,7 +51,7 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F burn1.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='TU') burn1.add_state('r', fix_initial=True, fix_final=False, defect_scaler=1.0, rate_source='r_dot', units='DU') - burn1.add_state('theta', fix_initial=True, fix_final=False, defect_scaler=1.0, + burn1.add_state('theta', fix_initial=True, fix_final=False, defect_scaler=1.0E-3, rate_source='theta_dot', units='rad') burn1.add_state('vr', fix_initial=True, fix_final=False, defect_scaler=1.0, rate_source='vr_dot', units='DU/TU') @@ -61,7 +61,7 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F rate_source='at_dot', units='DU/TU**2') burn1.add_state('deltav', fix_initial=True, fix_final=False, rate_source='deltav_dot', units='DU/TU') - burn1.add_control('u1', rate_continuity=True, rate2_continuity=True, units='deg', scaler=0.01, + burn1.add_control('u1', rate_continuity=False, rate2_continuity=False, units='deg', scaler=0.01, rate_continuity_scaler=0.001, rate2_continuity_scaler=0.001, lower=-30, upper=30) # Second Phase (Coast) @@ -70,7 +70,7 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F coast.set_time_options(initial_bounds=(0.5, 20), duration_bounds=(.5, 50), duration_ref=50, units='TU') coast.add_state('r', fix_initial=False, fix_final=False, defect_scaler=1.0, rate_source='r_dot', units='DU') - coast.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=1.0, + coast.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=1.0E-3, rate_source='theta_dot', units='rad') coast.add_state('vr', fix_initial=False, fix_final=False, defect_scaler=1.0, rate_source='vr_dot', units='DU/TU') @@ -94,11 +94,11 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F burn2.set_time_options(initial_bounds=(1.0, 60), duration_bounds=(-10.0, -0.5), initial_ref=10, units='TU') - burn2.add_state('r', fix_initial=True, fix_final=False, defect_scaler=1.0, + burn2.add_state('r', fix_initial=True, fix_final=False, defect_scaler=1.0E-2, rate_source='r_dot', units='DU') - burn2.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=1.0, + burn2.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=1000.0, rate_source='theta_dot', units='rad') - burn2.add_state('vr', fix_initial=True, fix_final=False, defect_scaler=1.0, + burn2.add_state('vr', fix_initial=True, fix_final=False, defect_scaler=1.0E-2, rate_source='vr_dot', units='DU/TU') burn2.add_state('vt', fix_initial=True, fix_final=False, defect_scaler=1.0, rate_source='vt_dot', units='DU/TU') @@ -109,16 +109,16 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F burn2.add_objective('deltav', loc='initial', scaler=100.0) - burn2.add_control('u1', rate_continuity=True, rate2_continuity=True, units='deg', + burn2.add_control('u1', rate_continuity=False, rate2_continuity=False, units='deg', scaler=0.01, lower=-180, upper=180) else: burn2.set_time_options(initial_bounds=(0.5, 50), duration_bounds=(.5, 10), initial_ref=10, units='TU') - burn2.add_state('r', fix_initial=False, fix_final=True, defect_scaler=1.0, + burn2.add_state('r', fix_initial=False, fix_final=True, defect_scaler=1.0E-2, rate_source='r_dot', targets=['r'], units='DU') - burn2.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=1.0, + burn2.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=1000.0, rate_source='theta_dot', targets=['theta'], units='rad') - burn2.add_state('vr', fix_initial=False, fix_final=True, defect_scaler=1.0, + burn2.add_state('vr', fix_initial=False, fix_final=True, defect_scaler=1.0E-2, rate_source='vr_dot', targets=['vr'], units='DU/TU') burn2.add_state('vt', fix_initial=False, fix_final=True, defect_scaler=1.0, rate_source='vt_dot', targets=['vt'], units='DU/TU') @@ -129,7 +129,7 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F burn2.add_objective('deltav', loc='final', scaler=100.0) - burn2.add_control('u1', rate_continuity=True, rate2_continuity=True, units='deg', + burn2.add_control('u1', rate_continuity=False, rate2_continuity=False, units='deg', scaler=0.01, lower=-180, upper=180, targets=['u1']) burn1.add_timeseries_output('pos_x') @@ -227,8 +227,6 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas' p.driver.opt_settings['mu_strategy'] = 'monotone' p.driver.opt_settings['derivative_test'] = 'first-order' - if show_output: - p.driver.opt_settings['print_level'] = 0 traj = make_traj(transcription=transcription, transcription_order=transcription_order, compressed=compressed, connected=connected, @@ -239,14 +237,14 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP # Needed to move the direct solver down into the phases for use with MPI. # - After moving down, used fewer iterations (about 30 less) - p.setup(check=True) + p.setup(check=False) # Set Initial Guesses traj.set_parameter_val('c', val=1.5, units='DU/TU') - burn1 = p.model.traj.phases.burn1 - burn2 = p.model.traj.phases.burn2 - coast = p.model.traj.phases.coast + burn1 = p.model._get_subsystem('traj.phases.burn1') + burn2 = p.model._get_subsystem('traj.phases.burn2') + coast = p.model._get_subsystem('traj.phases.coast') if burn1 in p.model.traj.phases._subsystems_myproc: burn1.set_time_val(initial=0.0, duration=2.25) @@ -288,6 +286,11 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP burn2.set_state_val('deltav', [0.1, 0.2]) burn2.set_control_val('u1', [0, 0]) + p.final_setup() + if burn2 in p.model.traj.phases._subsystems_myproc: + print(burn2.get_val('t_initial')) + print(burn2.get_val('t_duration')) + if run_driver or simulate: dm.run_problem(p, run_driver=run_driver, simulate=simulate, restart=restart, make_plots=True, solution_record_file=solution_record_file, simulation_record_file=simulation_record_file) diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py index 2902f24af..1793a16a7 100644 --- a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py +++ b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_mpi.py @@ -20,7 +20,7 @@ def test_ex_two_burn_orbit_raise_mpi(self): CONNECTED = False - p = two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3, + p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3, compressed=False, optimizer=optimizer, simulate=True, connected=CONNECTED, show_output=False) @@ -45,7 +45,7 @@ def test_ex_two_burn_orbit_raise_connected_mpi(self): CONNECTED = True - p = two_burn_orbit_raise_problem(transcription='gauss-lobatto', transcription_order=3, + p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3, compressed=False, optimizer=optimizer, simulate=True, connected=CONNECTED, show_output=False) diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_finite_burn_eom.py b/dymos/examples/finite_burn_orbit_raise/test/test_finite_burn_eom.py index 15107ccf6..79d6b481b 100644 --- a/dymos/examples/finite_burn_orbit_raise/test/test_finite_burn_eom.py +++ b/dymos/examples/finite_burn_orbit_raise/test/test_finite_burn_eom.py @@ -78,7 +78,7 @@ def test_outputs(self): assert_near_equal(p['at_dot'], p['accel']**2 / p['c']) def test_partials(self): - cpd = self.p.check_partials(compact_print=False) + cpd = self.p.check_partials(compact_print=False, out_stream=None) assert_check_partials(cpd, atol=1.0E-5, rtol=2.0) diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py b/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py index 766c84742..1598b7297 100644 --- a/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py +++ b/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py @@ -53,7 +53,7 @@ def test_ex_two_burn_orbit_raise_connected(self): sim_case2 = om.CaseReader(sim_db).get_case('final') # Verify that the second case has the same inputs and outputs - assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8) + assert_cases_equal(sim_case1, sim_case2, tol=1.0E-7) def test_restart_from_solution_radau(self): optimizer = 'IPOPT' @@ -86,7 +86,7 @@ def test_restart_from_solution_radau(self): sim_case2 = om.CaseReader(sim_db).get_case('final') # Verify that the second case has the same inputs and outputs - assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8) + assert_cases_equal(sim_case1, sim_case2, tol=1.0E-7) @require_pyoptsparse(optimizer='IPOPT') @@ -152,8 +152,8 @@ def test_ex_two_burn_orbit_raise_connected(self): sim_case2 = om.CaseReader(sim_db).get_case('final') # Verify that the second case has the same inputs and outputs - assert_cases_equal(case1, case2, tol=1.0E-8) - assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8) + assert_cases_equal(case1, case2, tol=1.0E-7) + assert_cases_equal(sim_case1, sim_case2, tol=1.0E-7) @unittest.skipUnless(MPI, "MPI is required.") def test_restart_from_solution_radau_to_connected(self): @@ -178,17 +178,17 @@ def test_restart_from_solution_radau_to_connected(self): sim_case1 = om.CaseReader(sim_db).get_case('final') # Run again without an actual optimizer - p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3, - compressed=False, optimizer=optimizer, run_driver=False, - show_output=False, restart=sol_db, - connected=True, solution_record_file='dymos_solution2.db', - simulation_record_file='dymos_simulation2.db') + p2 = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3, + compressed=False, optimizer=optimizer, run_driver=False, + show_output=False, restart=sol_db, + connected=True, solution_record_file='dymos_solution2.db', + simulation_record_file='dymos_simulation2.db') sol_db = 'dymos_solution2.db' sim_db = 'dymos_simulation2.db' if om_version()[0] > (3, 34, 2): - sol_db = p.get_outputs_dir() / sol_db - sim_db = p.model.traj.sim_prob.get_outputs_dir() / sim_db + sol_db = p2.get_outputs_dir() / sol_db + sim_db = p2.model.traj.sim_prob.get_outputs_dir() / sim_db case2 = om.CaseReader(sol_db).get_case('final') sim_case2 = om.CaseReader(sim_db).get_case('final') diff --git a/dymos/examples/oscillator/test/test_oscillator_vector_states.py b/dymos/examples/oscillator/test/test_oscillator_vector_states.py index 924d19afa..d56d4e4d6 100644 --- a/dymos/examples/oscillator/test/test_oscillator_vector_states.py +++ b/dymos/examples/oscillator/test/test_oscillator_vector_states.py @@ -7,7 +7,7 @@ import numpy as np -# @use_tempdirs +@use_tempdirs class TestDocOscillator(unittest.TestCase): def test_matrix_param(self): diff --git a/dymos/examples/racecar/test/test_racecar.py b/dymos/examples/racecar/test/test_racecar.py index df5ecffad..612be00ea 100644 --- a/dymos/examples/racecar/test/test_racecar.py +++ b/dymos/examples/racecar/test/test_racecar.py @@ -35,8 +35,7 @@ def test_racecar_for_docs(self): s_final = track.get_total_length() 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(num_nodes=100)} + 'gauss-lobatto': dm.GaussLobatto(num_segments=50, order=3, compressed=True)} for tx_name, tx in txs.items(): with self.subTest(tx_name): diff --git a/dymos/examples/water_rocket/phases.py b/dymos/examples/water_rocket/phases.py index 43cd63ef6..59fc4ce72 100644 --- a/dymos/examples/water_rocket/phases.py +++ b/dymos/examples/water_rocket/phases.py @@ -16,7 +16,7 @@ def new_propelled_ascent_phase(transcription): propelled_ascent.add_state('r', units='m', rate_source='eom.r_dot', fix_initial=True, fix_final=False, ref=100.0, defect_ref=100.0) propelled_ascent.add_state('h', units='m', rate_source='eom.h_dot', targets=['atmos.h'], - fix_initial=True, fix_final=False, ref=100.0, defect_ref=100.0) + fix_initial=True, fix_final=False, ref=10.0, defect_ref=10.0) propelled_ascent.add_state('gam', units='deg', rate_source='eom.gam_dot', targets=['eom.gam'], fix_initial=False, fix_final=False, lower=0, upper=85.0, ref=90) propelled_ascent.add_state('v', units='m/s', rate_source='eom.v_dot', targets=['dynamic_pressure.v', 'eom.v'], @@ -27,7 +27,7 @@ def new_propelled_ascent_phase(transcription): lower=1.02) propelled_ascent.add_state('V_w', units='L', rate_source='water_engine.Vdot', targets=['water_engine.V_w', 'mass_adder.V_w'], - fix_initial=False, fix_final=True, ref=10, defect_ref=10, + fix_initial=False, fix_final=True, ref=10., defect_ref=0.1, lower=0) propelled_ascent.add_parameter( @@ -94,7 +94,7 @@ def new_descent_phase(transcription): def new_water_rocket_trajectory(objective): - tx_prop = dm.Radau(num_segments=50, order=3, compressed=True) + tx_prop = dm.Radau(num_segments=20, order=3, compressed=True) tx_bal = dm.Radau(num_segments=10, order=3, compressed=True) tx_desc = dm.Radau(num_segments=10, order=3, compressed=True) traj = dm.Trajectory() @@ -112,7 +112,7 @@ def new_water_rocket_trajectory(objective): if objective == 'height': ballistic_ascent.add_objective('h', loc='final', ref=-1.0) elif objective == 'range': - descent.add_objective('r', loc='final', ref=-100.0) + descent.add_objective('r', loc='final', ref=-1.0) else: raise ValueError(f"objective='{objective}' is not defined. Try using 'height' or 'range'") diff --git a/dymos/examples/water_rocket/test/test_water_rocket.py b/dymos/examples/water_rocket/test/test_water_rocket.py index a8329e276..3882a46e0 100644 --- a/dymos/examples/water_rocket/test/test_water_rocket.py +++ b/dymos/examples/water_rocket/test/test_water_rocket.py @@ -31,8 +31,10 @@ def test_water_rocket_height_for_docs(self): p.driver = om.pyOptSparseDriver(optimizer='IPOPT', print_results=False) p.driver.opt_settings['print_level'] = 0 - p.driver.opt_settings['max_iter'] = 500 - p.driver.opt_settings['mu_strategy'] = 'monotone' + p.driver.opt_settings['max_iter'] = 1000 + p.driver.opt_settings['mu_strategy'] = 'adaptive' + p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' # for faster convergence + p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas' p.driver.declare_coloring(tol=1.0E-12) # Finish Problem Setup @@ -41,7 +43,7 @@ def test_water_rocket_height_for_docs(self): p.setup() set_sane_initial_guesses(phases) - dm.run_problem(p, run_driver=True, simulate=True) + dm.run_problem(p, run_driver=True, simulate=True, make_plots=True) summary = summarize_results(p) for key, entry in summary.items(): @@ -54,7 +56,7 @@ def test_water_rocket_height_for_docs(self): plot_trajectory(p, exp_out) plot_states(p, exp_out) - # plt.show() + plt.show() # Check results (tolerance is relative unless value is zero) assert_near_equal(summary['Launch angle'].value, 85, 0.01) @@ -71,8 +73,10 @@ def test_water_rocket_range_for_docs(self): p.driver = om.pyOptSparseDriver(optimizer='IPOPT') p.driver.opt_settings['print_level'] = 0 - p.driver.opt_settings['max_iter'] = 300 - p.driver.opt_settings['mu_strategy'] = 'monotone' + p.driver.opt_settings['max_iter'] = 1000 + p.driver.opt_settings['mu_strategy'] = 'adaptive' + p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' # for faster convergence + p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas' p.driver.declare_coloring(tol=1.0E-12) # Finish Problem Setup diff --git a/dymos/examples/water_rocket/water_engine_comp.py b/dymos/examples/water_rocket/water_engine_comp.py index 9eae3d5c9..5f109d71a 100644 --- a/dymos/examples/water_rocket/water_engine_comp.py +++ b/dymos/examples/water_rocket/water_engine_comp.py @@ -34,6 +34,9 @@ def setup(self): subsys=_WaterThrust(num_nodes=nn), promotes=['rho_w', 'A_out', 'F']) + self.set_input_defaults('A_out', val=np.pi*13e-3**2/4*np.ones(nn)) + self.set_input_defaults('p', val=6.5e5*np.ones(nn)) + self.connect('water_exhaust_speed.v_out', 'water_flow_rate.v_out') self.connect('water_exhaust_speed.v_out', 'water_thrust.v_out') @@ -187,5 +190,4 @@ def compute_partials(self, inputs, partials): p = om.Problem() p.model = WaterEngine(num_nodes=1) p.setup() - p.check_config(checks=['unconnected_inputs'], out_file=None) p.final_setup() diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index e0abbbaee..b83fa50ee 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -2350,73 +2350,6 @@ def _check_parameter_options(self): f"phase '{self.name}': {', '.join(invalid_options)}", RuntimeWarning) - def interpolate(self, xs=None, ys=None, nodes='all', kind='linear', axis=0): - """ - Return an array of values on interpolated to the given node subset of the phase. - - Parameters - ---------- - xs : ndarray or Sequence or None - Array of integration variable values. - ys : ndarray or Sequence or None - Array of control/state/parameter values. - nodes : str or None - The name of the node subset. - kind : str - Specifies the kind of interpolation, as per the scipy.interpolate package. - One of ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic' - where 'zero', 'slinear', 'quadratic' and 'cubic' refer to a spline - interpolation of zeroth, first, second or third order) or as an - integer specifying the order of the spline interpolator to use. - Default is 'linear'. - axis : int - Specifies the axis along which interpolation should be performed. Default is - the first axis (0). - - Returns - ------- - np.array - The values of y interpolated at nodes of the specified type. - """ - om.issue_warning('phase.interpolate has been deprecated and will be removed from Dymos ' - '2.0.0. Use phase.interp instead, which uses a different order for the ' - 'arguments but is more terse and can interpolate polynomial control ' - 'values.', category=om.OMDeprecationWarning) - - if not isinstance(ys, Iterable): - raise ValueError('ys must be provided as an Iterable of length at least 2.') - if nodes not in ('col', 'all', 'state_disc', 'state_input', 'control_disc', - 'control_input', 'segment_ends'): - raise ValueError("nodes must be one of 'col', 'all', 'state_disc', " - "'state_input', 'control_disc', 'control_input', or 'segment_ends'") - if xs is None: - if len(ys) != 2: - raise ValueError('xs may only be unspecified when len(ys)=2') - if kind != 'linear': - raise ValueError('kind must be linear when xs is unspecified.') - xs = [-1, 1] - elif len(xs) != np.prod(np.asarray(xs).shape): - raise ValueError('xs must be viewable as a 1D array') - - gd = self.options['transcription'].grid_data - - if gd is None: - raise RuntimeError('interpolate cannot be called until the associated ' - 'problem has been setup') - - node_locations = gd.node_ptau[gd.subset_node_indices[nodes]] - # Affine transform xs into tau space [-1, 1] - _xs = np.asarray(xs).ravel() - m = 2.0 / (_xs[-1] - _xs[0]) - b = 1.0 - (m * _xs[-1]) - taus = m * _xs + b - interpfunc = interpolate.interp1d(taus, ys, axis=axis, kind=kind, - bounds_error=False, fill_value='extrapolate') - res = np.atleast_2d(interpfunc(node_locations)) - if res.shape[0] == 1: - res = res.T - return res - def interp(self, name=None, ys=None, xs=None, nodes=None, kind='linear', axis=0): """ Interpolate values onto the given subset of nodes in the phase. @@ -2435,14 +2368,12 @@ def interp(self, name=None, ys=None, xs=None, nodes=None, kind='linear', axis=0) Array of control/state/parameter values. xs : ndarray or Sequence or None Array of integration variable values. - nodes : str or None - The name of the node subset or None (default). + nodes : str or Sequence or None + The name of the node subset, a set of nodes in phase tau space ([-1, 1] across the phase), or None (default). kind : str Specifies the kind of interpolation, as per the scipy.interpolate package. - One of ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic' - where 'zero', 'slinear', 'quadratic' and 'cubic' refer to a spline - interpolation of zeroth, first, second or third order) or as an - integer specifying the order of the spline interpolator to use. + One of ('linear', 'quadratic', 'cubic'). Any other value will result in + linear interpolation and is allowed for backward compatibility. Default is 'linear'. axis : int Specifies the axis along which interpolation should be performed. Default is @@ -2455,11 +2386,24 @@ def interp(self, name=None, ys=None, xs=None, nodes=None, kind='linear', axis=0) """ if not isinstance(ys, Iterable): raise ValueError('ys must be provided as an Iterable of length at least 2.') - if nodes not in ('col', 'all', 'state_disc', 'state_input', 'control_disc', - 'control_input', 'segment_ends', None): - raise ValueError("nodes must be one of 'col', 'all', 'state_disc', " - "'state_input', 'control_disc', 'control_input', 'segment_ends', or " - "None.") + if isinstance(nodes, str): + if nodes not in ('col', 'all', 'state_disc', 'state_input', 'control_disc', + 'control_input', 'segment_ends', None): + raise ValueError("nodes must be one of 'col', 'all', 'state_disc', " + "'state_input', 'control_disc', 'control_input', 'segment_ends', or " + "None.") + + order_map = {'linear': 1, + 'quadratic': 2, + 'cubic': 3} + + if isinstance(kind, str): + _kind = order_map.get(kind, 1) + elif isinstance(kind, int): + _kind = kind + else: + raise ValueError("kind must be an integer of the spline order or one of 'linear', 'quadratic', or 'cubic'. " + " Any other string value will result in linear interpolation.") if xs is None: if len(ys) != 2: @@ -2492,17 +2436,28 @@ def interp(self, name=None, ys=None, xs=None, nodes=None, kind='linear', axis=0) raise ValueError('Could not find a state, control, or polynomial control named ' f'{name} to be interpolated.\nPlease explicitly specify the ' f'node subset onto which this value should be interpolated.') - else: + elif isinstance(nodes, str): node_locations = gd.node_ptau[gd.subset_node_indices[nodes]] + else: + node_locations = nodes # Affine transform xs into tau space [-1, 1] _xs = np.asarray(xs).ravel() - m = 2.0 / (_xs[-1] - _xs[0]) - b = 1.0 - (m * _xs[-1]) - taus = m * _xs + b - interpfunc = interpolate.interp1d(taus, ys, axis=axis, kind=kind, - bounds_error=False, fill_value='extrapolate') + reversed_xs = _xs[0] > _xs[1] + _xs, unique_idxs = np.unique(_xs, return_index=True) + if reversed_xs: + _xs = _xs[::-1] + unique_idxs = unique_idxs[::-1] + _ys = np.asarray(ys) + _ys = np.atleast_2d(_ys).T if len(_ys.shape) == 1 else _ys + + dptau_dt = 2.0 / (_xs[-1] - _xs[0]) + b = 1.0 - (dptau_dt * _xs[-1]) + taus = dptau_dt * _xs + b + + interpfunc = interpolate.make_interp_spline(taus, _ys[unique_idxs, ...], k=_kind) res = np.atleast_2d(interpfunc(node_locations)) + if res.shape[0] == 1: res = res.T return res diff --git a/dymos/phase/test/test_interp.py b/dymos/phase/test/test_interp.py index 9411dbe40..a0e0849ac 100644 --- a/dymos/phase/test/test_interp.py +++ b/dymos/phase/test/test_interp.py @@ -130,6 +130,15 @@ def test_underspecified(self): self.assertEqual(str(e.exception), expected) + def test_nodes_sequence(self): + tx = dm.GaussLobatto(num_segments=6, order=3, compressed=True) + phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx) + + interp_vals = phase.interp(ys=[0, 5], nodes=np.linspace(-1, 1, 100)) + check_vals = np.atleast_2d(np.linspace(0, 5, 100)).T + + assert_near_equal(interp_vals, check_vals) + if __name__ == '__main__': # pragma: no cover unittest.main() diff --git a/dymos/test/test_upgrade_guide.py b/dymos/test/test_upgrade_guide.py deleted file mode 100644 index b40b02fc4..000000000 --- a/dymos/test/test_upgrade_guide.py +++ /dev/null @@ -1,777 +0,0 @@ -import os -import unittest - -import openmdao.api as om - -import dymos as dm - -from openmdao.utils.assert_utils import assert_near_equal -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse -from dymos.utils.misc import om_version - - -@use_tempdirs -class TestUpgrade_0_16_0(unittest.TestCase): - - @require_pyoptsparse(optimizer='SLSQP') - def test_parameters(self): - """ - # upgrade_doc: begin set_val - p.set_val('traj.phase0.design_parameters:thrust', 2.1, units='MN') - # upgrade_doc: end set_val - # upgrade_doc: begin parameter_timeseries - thrust = p.get_val('traj.phase0.timeseries.design_parameters:thrust') - # upgrade_doc: end parameter_timeseries - """ - import openmdao.api as om - import dymos as dm - - # - # Setup and solve the optimal control problem - # - p = om.Problem(model=om.Group()) - p.driver = om.pyOptSparseDriver() - p.driver.declare_coloring(tol=1.0E-12) - - from dymos.examples.ssto.launch_vehicle_ode import LaunchVehicleODE - - # - # Initialize our Trajectory and Phase - # - traj = dm.Trajectory() - - phase = dm.Phase(ode_class=LaunchVehicleODE, - transcription=dm.GaussLobatto(num_segments=12, order=3, compressed=False)) - - traj.add_phase('phase0', phase) - p.model.add_subsystem('traj', traj) - - # - # Set the options for the variables - # - phase.set_time_options(fix_initial=True, duration_bounds=(10, 500)) - - phase.add_state('x', fix_initial=True, ref=1.0E5, defect_ref=10000.0, - rate_source='xdot') - phase.add_state('y', fix_initial=True, ref=1.0E5, defect_ref=10000.0, - rate_source='ydot') - phase.add_state('vx', fix_initial=True, ref=1.0E3, defect_ref=1000.0, - rate_source='vxdot') - phase.add_state('vy', fix_initial=True, ref=1.0E3, defect_ref=1000.0, - rate_source='vydot') - phase.add_state('m', fix_initial=True, ref=1.0E3, defect_ref=100.0, - rate_source='mdot') - - phase.add_control('theta', units='rad', lower=-1.57, upper=1.57) - phase.add_parameter('thrust', units='N', opt=False, val=2100000.0) - - # - # Set the options for our constraints and objective - # - phase.add_boundary_constraint('y', loc='final', equals=1.85E5, linear=True) - phase.add_boundary_constraint('vx', loc='final', equals=7796.6961) - phase.add_boundary_constraint('vy', loc='final', equals=0) - - phase.add_objective('time', loc='final', scaler=0.01) - - p.model.linear_solver = om.DirectSolver() - - # - # Setup and set initial values - # - p.setup(check=True) - - p.set_val('traj.phase0.t_initial', 0.0) - p.set_val('traj.phase0.t_duration', 150.0) - p.set_val('traj.phase0.states:x', phase.interpolate(ys=[0, 1.15E5], nodes='state_input')) - p.set_val('traj.phase0.states:y', phase.interpolate(ys=[0, 1.85E5], nodes='state_input')) - p.set_val('traj.phase0.states:vx', phase.interpolate(ys=[0, 7796.6961], nodes='state_input')) - p.set_val('traj.phase0.states:vy', phase.interpolate(ys=[1.0E-6, 0], nodes='state_input')) - p.set_val('traj.phase0.states:m', phase.interpolate(ys=[117000, 1163], nodes='state_input')) - p.set_val('traj.phase0.controls:theta', phase.interpolate(ys=[1.5, -0.76], nodes='control_input')) - - # upgrade_doc: begin set_val - p.set_val('traj.phase0.parameters:thrust', 2.1, units='MN') - # upgrade_doc: end set_val - - # - # Solve the Problem - # - dm.run_problem(p) - - # - # Check the results. - # - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 143, tolerance=0.05) - assert_near_equal(p.get_val('traj.phase0.timeseries.y')[-1], 1.85E5, 1e-4) - assert_near_equal(p.get_val('traj.phase0.timeseries.vx')[-1], 7796.6961, 1e-4) - assert_near_equal(p.get_val('traj.phase0.timeseries.vy')[-1], 0, 1e-4) - - # upgrade_doc: begin parameter_timeseries - thrust = p.get_val('traj.phase0.parameter_vals:thrust') - # upgrade_doc: end parameter_timeseries - assert_near_equal(thrust.ravel(), 2.1E6) - - def test_parameter_no_timeseries(self): - import openmdao.api as om - import dymos as dm - from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE - - p = om.Problem(model=om.Group()) - - p.driver = om.ScipyOptimizeDriver() - p.driver.declare_coloring() - - phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.GaussLobatto(num_segments=8, order=3, compressed=True)) - - p.model.add_subsystem('phase0', phase) - - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - - phase.add_state('x', fix_initial=True, fix_final=True, solve_segments=False) - - phase.add_state('y', fix_initial=True, fix_final=True, solve_segments=False) - - phase.add_state('v', fix_initial=True, fix_final=False, solve_segments=False) - - phase.add_control('theta', continuity=True, rate_continuity=True, opt=True, - units='deg', lower=0.01, upper=179.9, ref=1, ref0=0) - - # upgrade_doc: begin parameter_no_timeseries - phase.add_parameter('g', opt=False, units='m/s**2', include_timeseries=False) - # upgrade_doc: end parameter_no_timeseries - - # Minimize time at the end of the phase - phase.add_objective('time_phase', loc='final', scaler=10) - - p.model.linear_solver = om.DirectSolver() - p.setup(check=True) - - p['phase0.t_initial'] = 0.0 - p['phase0.t_duration'] = 2.0 - - p['phase0.states:x'] = phase.interpolate(ys=[0, 10], nodes='state_input') - p['phase0.states:y'] = phase.interpolate(ys=[10, 5], nodes='state_input') - p['phase0.states:v'] = phase.interpolate(ys=[0, 9.9], nodes='state_input') - p['phase0.controls:theta'] = phase.interpolate(ys=[5, 100], nodes='control_input') - p['phase0.parameters:g'] = 9.80665 - - p.run_driver() - - with self.assertRaises(KeyError): - p.get_val('phase0.timeseries.g}') - - def test_simplified_ode_timeseries_output(self): - """ - # upgrade_doc: begin simplified_ode_output_timeseries - phase.add_timeseries_output('tas_comp.TAS', shape=(1,), units='m/s') - # upgrade_doc: end simplified_ode_output_timeseries - """ - import openmdao.api as om - import dymos as dm - from dymos.examples.aircraft_steady_flight.aircraft_ode import AircraftODE - - p = om.Problem(model=om.Group()) - p.driver = om.ScipyOptimizeDriver() - p.driver.declare_coloring() - - # Pass Reference Area from an external source - assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp()) - assumptions.add_output('S', val=427.8, units='m**2') - assumptions.add_output('mass_empty', val=1.0, units='kg') - assumptions.add_output('mass_payload', val=1.0, units='kg') - - transcription = dm.GaussLobatto(num_segments=1, - order=13, - compressed=False) - phase = dm.Phase(ode_class=AircraftODE, transcription=transcription) - p.model.add_subsystem('phase0', phase) - - phase.set_time_options(fix_initial=True, - duration_bounds=(3600, 3600), - duration_ref=3600) - - phase.set_time_options(fix_initial=True, - duration_bounds=(3600, 3600), - duration_ref=3600) - - phase.add_state('range', units='km', fix_initial=True, fix_final=False, scaler=0.01, - rate_source='range_rate_comp.dXdt:range', - defect_scaler=0.01) - phase.add_state('mass_fuel', units='kg', fix_final=True, upper=20000.0, lower=0.0, - rate_source='propulsion.dXdt:mass_fuel', - scaler=1.0E-4, defect_scaler=1.0E-2) - phase.add_state('alt', - rate_source='climb_rate', - units='km', fix_initial=True) - - phase.add_control('mach', targets=['tas_comp.mach', 'aero.mach'], units=None, opt=False) - - phase.add_control('climb_rate', targets=['gam_comp.climb_rate'], units='m/s', opt=False) - - phase.add_parameter('S', - targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'], - units='m**2') - - phase.add_parameter('mass_empty', targets=['mass_comp.mass_empty'], units='kg') - phase.add_parameter('mass_payload', targets=['mass_comp.mass_payload'], units='kg') - - phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0, shape=(1,)) - - # upgrade_doc: begin simplified_ode_output_timeseries - phase.add_timeseries_output('tas_comp.TAS') - # upgrade_doc: end simplified_ode_output_timeseries - - p.model.connect('assumptions.S', 'phase0.parameters:S') - p.model.connect('assumptions.mass_empty', 'phase0.parameters:mass_empty') - p.model.connect('assumptions.mass_payload', 'phase0.parameters:mass_payload') - - phase.add_objective('time', loc='final', ref=3600) - - p.model.linear_solver = om.DirectSolver() - - p.setup() - - p['phase0.t_initial'] = 0.0 - p['phase0.t_duration'] = 1.515132 * 3600.0 - p['phase0.states:range'] = phase.interpolate(ys=(0, 1296.4), nodes='state_input') - p['phase0.states:mass_fuel'] = phase.interpolate(ys=(12236.594555, 0), nodes='state_input') - p['phase0.states:alt'] = 5.0 - p['phase0.controls:mach'] = 0.8 - p['phase0.controls:climb_rate'] = 0.0 - - p['assumptions.S'] = 427.8 - p['assumptions.mass_empty'] = 0.15E6 - p['assumptions.mass_payload'] = 84.02869 * 400 - - dm.run_problem(p) - - time = p.get_val('phase0.timeseries.time') - tas = p.get_val('phase0.timeseries.TAS', units='km/s') - range = p.get_val('phase0.timeseries.range') - - assert_near_equal(range, tas*time, tolerance=1.0E-4) - - exp_out = phase.simulate() - - time = exp_out.get_val('phase0.timeseries.time') - tas = exp_out.get_val('phase0.timeseries.TAS', units='km/s') - range = exp_out.get_val('phase0.timeseries.range') - - assert_near_equal(range, tas*time, tolerance=1.0E-4) - - @require_pyoptsparse(optimizer='SLSQP') - def test_glob_timeseries_outputs(self): - """ - # upgrade_doc: begin glob_timeseries_outputs - phase.add_timeseries_output('aero.mach', shape=(1,), units=None) - phase.add_timeseries_output('aero.CD0', shape=(1,), units=None) - phase.add_timeseries_output('aero.kappa', shape=(1,), units=None) - phase.add_timeseries_output('aero.CLa', shape=(1,), units=None) - phase.add_timeseries_output('aero.CL', shape=(1,), units=None) - phase.add_timeseries_output('aero.CD', shape=(1,), units=None) - phase.add_timeseries_output('aero.q', shape=(1,), units='N/m**2') - phase.add_timeseries_output('aero.f_lift', shape=(1,), units='N') - phase.add_timeseries_output('aero.f_drag', shape=(1,), units='N') - # upgrade_doc: end glob_timeseries_outputs - """ - from dymos.examples.min_time_climb.min_time_climb_ode import MinTimeClimbODE - - p = om.Problem(model=om.Group()) - - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'SLSQP' - p.driver.declare_coloring(tol=1.0E-12) - - traj = dm.Trajectory() - - phase = dm.Phase(ode_class=MinTimeClimbODE, - transcription=dm.GaussLobatto(num_segments=15, compressed=True)) - - traj.add_phase('phase0', phase) - - p.model.add_subsystem('traj', traj) - - phase.set_time_options(fix_initial=True, duration_bounds=(50, 400), - duration_ref=100.0) - - phase.add_state('r', fix_initial=True, lower=0, upper=1.0E6, - ref=1.0E3, defect_ref=1.0E3, units='m', - rate_source='flight_dynamics.r_dot') - - phase.add_state('h', fix_initial=True, lower=0, upper=20000.0, - ref=1.0E2, defect_ref=1.0E2, units='m', - rate_source='flight_dynamics.h_dot') - - phase.add_state('v', fix_initial=True, lower=10.0, - ref=1.0E2, defect_ref=1.0E2, units='m/s', - rate_source='flight_dynamics.v_dot') - - phase.add_state('gam', fix_initial=True, lower=-1.5, upper=1.5, - ref=1.0, defect_ref=1.0, units='rad', - rate_source='flight_dynamics.gam_dot') - - phase.add_state('m', fix_initial=True, lower=10.0, upper=1.0E5, - ref=1.0E3, defect_ref=1.0E3, units='kg', - rate_source='prop.m_dot') - - phase.add_control('alpha', units='deg', lower=-8.0, upper=8.0, scaler=1.0, - rate_continuity=True, rate_continuity_scaler=100.0, - rate2_continuity=False, targets=['alpha']) - - phase.add_parameter('S', val=49.2386, units='m**2', opt=False, targets=['S']) - phase.add_parameter('Isp', val=1600.0, units='s', opt=False, targets=['Isp']) - phase.add_parameter('throttle', val=1.0, opt=False, targets=['throttle']) - - phase.add_boundary_constraint('h', loc='final', equals=20000, scaler=1.0E-3, units='m') - phase.add_boundary_constraint('aero.mach', loc='final', equals=1.0, shape=(1,)) - phase.add_boundary_constraint('gam', loc='final', equals=0.0, units='rad') - - phase.add_path_constraint(name='h', lower=100.0, upper=20000, ref=20000) - phase.add_path_constraint(name='aero.mach', lower=0.1, upper=1.8, shape=(1,)) - - phase.add_objective('time', loc='final', ref=1.0) - - p.model.linear_solver = om.DirectSolver() - - # upgrade_doc: begin glob_timeseries_outputs - phase.add_timeseries_output('aero.*') - # upgrade_doc: end glob_timeseries_outputs - - p.setup(check=True) - - p['traj.phase0.t_initial'] = 0.0 - p['traj.phase0.t_duration'] = 500 - - p['traj.phase0.states:r'] = phase.interpolate(ys=[0.0, 50000.0], nodes='state_input') - p['traj.phase0.states:h'] = phase.interpolate(ys=[100.0, 20000.0], nodes='state_input') - p['traj.phase0.states:v'] = phase.interpolate(ys=[135.964, 283.159], nodes='state_input') - p['traj.phase0.states:gam'] = phase.interpolate(ys=[0.0, 0.0], nodes='state_input') - p['traj.phase0.states:m'] = phase.interpolate(ys=[19030.468, 10000.], nodes='state_input') - p['traj.phase0.controls:alpha'] = phase.interpolate(ys=[0.0, 0.0], nodes='control_input') - - # - # Solve for the optimal trajectory - # - p.run_model() - - outputs = p.model.list_outputs(units=True, out_stream=None, prom_name=True) - op_dict = {options['prom_name']: options['units'] for abs_name, options in outputs} - - for name, units in [('mach', None), ('CD0', None), ('kappa', None), ('CLa', None), - ('CL', None), ('CD', None), ('q', 'N/m**2'), ('f_lift', 'N'), - ('f_drag', 'N')]: - self.assertEqual(op_dict[f'traj.phase0.timeseries.{name}'], units) - - @require_pyoptsparse(optimizer='SLSQP') - def test_sequence_timeseries_outputs(self): - """ - # upgrade_doc: begin sequence_timeseries_outputs - phase.add_timeseries_output('aero.mach', shape=(1,), units=None) - phase.add_timeseries_output('aero.CD0', shape=(1,), units=None) - phase.add_timeseries_output('aero.kappa', shape=(1,), units=None) - phase.add_timeseries_output('aero.CLa', shape=(1,), units=None) - phase.add_timeseries_output('aero.CL', shape=(1,), units=None) - phase.add_timeseries_output('aero.CD', shape=(1,), units=None) - phase.add_timeseries_output('aero.q', shape=(1,), units='N/m**2') - phase.add_timeseries_output('aero.f_lift', shape=(1,), units='lbf') - phase.add_timeseries_output('aero.f_drag', shape=(1,), units='N') - phase.add_timeseries_output('prop.thrust', shape=(1,), units='lbf') - # upgrade_doc: end sequence_timeseries_outputs - # upgrade_doc: begin state_endpoint_values - final_range = p.get_val('traj.phase0.final_conditions.states:x0++') - final_alpha = p.get_val('traj.phase0.final_conditions.controls:alpha++') - # upgrade_doc: end state_endpoint_values - """ - from dymos.examples.min_time_climb.min_time_climb_ode import MinTimeClimbODE - - p = om.Problem(model=om.Group()) - - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'SLSQP' - p.driver.declare_coloring(tol=1.0E-12) - - traj = dm.Trajectory() - - phase = dm.Phase(ode_class=MinTimeClimbODE, - transcription=dm.GaussLobatto(num_segments=15, compressed=True)) - - traj.add_phase('phase0', phase) - - p.model.add_subsystem('traj', traj) - - phase.set_time_options(fix_initial=True, duration_bounds=(50, 400), - duration_ref=100.0) - - phase.add_state('r', fix_initial=True, lower=0, upper=1.0E6, - ref=1.0E3, defect_ref=1.0E3, units='m', - rate_source='flight_dynamics.r_dot') - - phase.add_state('h', fix_initial=True, lower=0, upper=20000.0, - ref=1.0E2, defect_ref=1.0E2, units='m', - rate_source='flight_dynamics.h_dot') - - phase.add_state('v', fix_initial=True, lower=10.0, - ref=1.0E2, defect_ref=1.0E2, units='m/s', - rate_source='flight_dynamics.v_dot') - - phase.add_state('gam', fix_initial=True, lower=-1.5, upper=1.5, - ref=1.0, defect_ref=1.0, units='rad', - rate_source='flight_dynamics.gam_dot') - - phase.add_state('m', fix_initial=True, lower=10.0, upper=1.0E5, - ref=1.0E3, defect_ref=1.0E3, units='kg', - rate_source='prop.m_dot') - - phase.add_control('alpha', units='deg', lower=-8.0, upper=8.0, scaler=1.0, - rate_continuity=True, rate_continuity_scaler=100.0, - rate2_continuity=False, targets=['alpha']) - - phase.add_parameter('S', val=49.2386, units='m**2', opt=False, targets=['S']) - phase.add_parameter('Isp', val=1600.0, units='s', opt=False, targets=['Isp']) - phase.add_parameter('throttle', val=1.0, opt=False, targets=['throttle']) - - phase.add_boundary_constraint('h', loc='final', equals=20000, scaler=1.0E-3, units='m') - phase.add_boundary_constraint('aero.mach', loc='final', equals=1.0, shape=(1,)) - phase.add_boundary_constraint('gam', loc='final', equals=0.0, units='rad') - - phase.add_path_constraint(name='h', lower=100.0, upper=20000, ref=20000) - phase.add_path_constraint(name='aero.mach', lower=0.1, upper=1.8, shape=(1,)) - - phase.add_objective('time', loc='final', ref=1.0) - - p.model.linear_solver = om.DirectSolver() - - # upgrade_doc: begin sequence_timeseries_outputs - phase.add_timeseries_output(['aero.*', 'prop.thrust'], - units={'aero.f_lift': 'lbf', 'prop.thrust': 'lbf'}) - # upgrade_doc: end sequence_timeseries_outputs - - p.setup(check=True) - - p['traj.phase0.t_initial'] = 0.0 - p['traj.phase0.t_duration'] = 500 - - p['traj.phase0.states:r'] = phase.interpolate(ys=[0.0, 50000.0], nodes='state_input') - p['traj.phase0.states:h'] = phase.interpolate(ys=[100.0, 20000.0], nodes='state_input') - p['traj.phase0.states:v'] = phase.interpolate(ys=[135.964, 283.159], nodes='state_input') - p['traj.phase0.states:gam'] = phase.interpolate(ys=[0.0, 0.0], nodes='state_input') - p['traj.phase0.states:m'] = phase.interpolate(ys=[19030.468, 10000.], nodes='state_input') - p['traj.phase0.controls:alpha'] = phase.interpolate(ys=[0.0, 0.0], nodes='control_input') - - p.run_model() - - # upgrade_doc: begin state_endpoint_values - final_range = p.get_val('traj.phase0.timeseries.r')[-1, ...] - final_alpha = p.get_val('traj.phase0.timeseries.alpha')[-1, ...] - # upgrade_doc: end state_endpoint_values - self.assertEqual(final_range, 50000.0) - self.assertEqual(final_alpha, 0.0) - - outputs = p.model.list_outputs(units=True, out_stream=None, prom_name=True) - op_dict = {options['prom_name']: options['units'] for abs_name, options in outputs} - - for name, units in [('mach', None), ('CD0', None), ('kappa', None), ('CLa', None), - ('CL', None), ('CD', None), ('q', 'N/m**2'), ('f_lift', 'lbf'), - ('f_drag', 'N'), ('thrust', 'lbf')]: - self.assertEqual(op_dict[f'traj.phase0.timeseries.{name}'], units) - - -@use_tempdirs -class TestUpgrade_0_17_0(unittest.TestCase): - - def test_tags(self): - """ - # upgrade_doc: begin tag_rate_source - self.add_output('xdot', val=np.zeros(nn), desc='velocity component in x', units='m/s') - self.add_output('ydot', val=np.zeros(nn), desc='velocity component in y', units='m/s') - self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2') - # upgrade_doc: end tag_rate_source - # upgrade_doc: begin declare_rate_source - phase.add_state('x', rate_source='xdot', fix_initial=True, fix_final=True) - phase.add_state('y', rate_source='ydot', fix_initial=True, fix_final=True) - phase.add_state('v', rate_source='vdot', fix_initial=True, fix_final=False) - # upgrade_doc: end declare_rate_source - """ - import numpy as np - import openmdao.api as om - - class BrachistochroneODE(om.ExplicitComponent): - - def initialize(self): - self.options.declare('num_nodes', types=int) - self.options.declare('static_gravity', types=(bool,), default=False, - desc='If True, treat gravity as a static (scalar) input, rather than ' - 'having different values at each node.') - - def setup(self): - nn = self.options['num_nodes'] - g_default_val = 9.80665 if self.options['static_gravity'] else 9.80665 * np.ones(nn) - - # Inputs - self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s') - - self.add_input('g', val=g_default_val, desc='grav. acceleration', units='m/s/s') - - self.add_input('theta', val=np.ones(nn), desc='angle of wire', units='rad') - - # upgrade_doc: begin tag_rate_source - self.add_output('xdot', val=np.zeros(nn), desc='velocity component in x', units='m/s', - tags=['dymos.state_rate_source:x', 'dymos.state_units:m']) - - self.add_output('ydot', val=np.zeros(nn), desc='velocity component in y', units='m/s', - tags=['dymos.state_rate_source:y', 'dymos.state_units:m']) - - self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2', - tags=['dymos.state_rate_source:v', 'dymos.state_units:m/s']) - # upgrade_doc: end tag_rate_source - - self.add_output('check', val=np.zeros(nn), desc='check solution: v/sin(theta) = constant', - units='m/s') - - # Setup partials - arange = np.arange(self.options['num_nodes']) - self.declare_partials(of='vdot', wrt='theta', rows=arange, cols=arange) - - self.declare_partials(of='xdot', wrt='v', rows=arange, cols=arange) - self.declare_partials(of='xdot', wrt='theta', rows=arange, cols=arange) - - self.declare_partials(of='ydot', wrt='v', rows=arange, cols=arange) - self.declare_partials(of='ydot', wrt='theta', rows=arange, cols=arange) - - self.declare_partials(of='check', wrt='v', rows=arange, cols=arange) - self.declare_partials(of='check', wrt='theta', rows=arange, cols=arange) - - if self.options['static_gravity']: - c = np.zeros(self.options['num_nodes']) - self.declare_partials(of='vdot', wrt='g', rows=arange, cols=c) - else: - self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange) - - def compute(self, inputs, outputs): - theta = inputs['theta'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - g = inputs['g'] - v = inputs['v'] - - outputs['vdot'] = g * cos_theta - outputs['xdot'] = v * sin_theta - outputs['ydot'] = -v * cos_theta - outputs['check'] = v / sin_theta - - def compute_partials(self, inputs, partials): - theta = inputs['theta'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - g = inputs['g'] - v = inputs['v'] - - partials['vdot', 'g'] = cos_theta - partials['vdot', 'theta'] = -g * sin_theta - - partials['xdot', 'v'] = sin_theta - partials['xdot', 'theta'] = v * cos_theta - - partials['ydot', 'v'] = -cos_theta - partials['ydot', 'theta'] = v * sin_theta - - partials['check', 'v'] = 1 / sin_theta - partials['check', 'theta'] = -v * cos_theta / sin_theta ** 2 - - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - - # - # Initialize the Problem and the optimization driver - # - p = om.Problem(model=om.Group()) - p.driver = om.ScipyOptimizeDriver() - p.driver.declare_coloring() - - # - # Create a trajectory and add a phase to it - # - traj = p.model.add_subsystem('traj', dm.Trajectory()) - - phase = traj.add_phase('phase0', - dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.GaussLobatto(num_segments=10))) - - # - # Set the variables - # - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - - # upgrade_doc: begin declare_rate_source - phase.add_state('x', fix_initial=True, fix_final=True) - phase.add_state('y', fix_initial=True, fix_final=True) - phase.add_state('v', fix_initial=True, fix_final=False) - # upgrade_doc: end declare_rate_source - - phase.add_control('theta', continuity=True, rate_continuity=True, - units='deg', lower=0.01, upper=179.9) - - phase.add_parameter('g', units='m/s**2', val=9.80665) - # - # Minimize time at the end of the phase - # - phase.add_objective('time', loc='final', scaler=10) - - p.model.linear_solver = om.DirectSolver() - - # - # Setup the Problem - # - p.setup() - - # - # Set the initial values - # - p['traj.phase0.t_initial'] = 0.0 - p['traj.phase0.t_duration'] = 2.0 - - p['traj.phase0.states:x'] = phase.interpolate(ys=[0, 10], nodes='state_input') - p['traj.phase0.states:y'] = phase.interpolate(ys=[10, 5], nodes='state_input') - p['traj.phase0.states:v'] = phase.interpolate(ys=[0, 9.9], nodes='state_input') - p['traj.phase0.controls:theta'] = phase.interpolate(ys=[5, 100.5], nodes='control_input') - - # - # Solve for the optimal trajectory - # - dm.run_problem(p) - - # Test the results - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3) - - -@use_tempdirs -class TestUpgrade_0_19_0(unittest.TestCase): - - @require_pyoptsparse(optimizer='SLSQP') - def _make_problem(self, transcription='gauss-lobatto', num_segments=8, transcription_order=3, - compressed=True, optimizer='SLSQP', run_driver=True, - force_alloc_complex=False, - solve_segments=False): - - p = om.Problem(model=om.Group()) - - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = optimizer - p.driver.declare_coloring(tol=1.0E-12) - - if transcription == 'gauss-lobatto': - t = dm.GaussLobatto(num_segments=num_segments, - order=transcription_order, - compressed=compressed) - elif transcription == 'radau-ps': - t = dm.Radau(num_segments=num_segments, - order=transcription_order, - compressed=compressed) - - # upgrade_doc: begin exec_comp_ode - def ode(num_nodes): - return om.ExecComp(['vdot = g * cos(theta)', - 'xdot = v * sin(theta)', - 'ydot = -v * cos(theta)'], - g={'val': 9.80665, 'units': 'm/s**2'}, - v={'shape': (num_nodes,), 'units': 'm/s'}, - theta={'shape': (num_nodes,), 'units': 'rad'}, - vdot={'shape': (num_nodes,), 'units': 'm/s**2'}, - xdot={'shape': (num_nodes,), 'units': 'm/s'}, - ydot={'shape': (num_nodes,), 'units': 'm/s'}, - has_diag_partials=True) - - phase = dm.Phase(ode_class=ode, transcription=t) - # upgrade_doc: end declare_rate_source - traj = dm.Trajectory() - p.model.add_subsystem('traj0', traj) - traj.add_phase('phase0', phase) - - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - - phase.add_state('x', fix_initial=True, fix_final=False, solve_segments=solve_segments, - rate_source='xdot') - phase.add_state('y', fix_initial=True, fix_final=False, solve_segments=solve_segments, - rate_source='ydot') - - # Note that by omitting the targets here Dymos will automatically attempt to connect - # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found. - phase.add_state('v', fix_initial=True, fix_final=False, solve_segments=solve_segments, - rate_source='vdot') - - phase.add_control('theta', - continuity=True, rate_continuity=True, - units='deg', lower=0.01, upper=179.9) - - phase.add_parameter('g', units='m/s**2', static_target=True) - - phase.add_boundary_constraint('x', loc='final', equals=10) - phase.add_boundary_constraint('y', loc='final', equals=5) - # Minimize time at the end of the phase - phase.add_objective('time_phase', loc='final', scaler=10) - - p.setup(check=['unconnected_inputs'], force_alloc_complex=force_alloc_complex) - - p['traj0.phase0.t_initial'] = 0.0 - p['traj0.phase0.t_duration'] = 2.0 - - p['traj0.phase0.states:x'] = phase.interpolate(ys=[0, 10], nodes='state_input') - p['traj0.phase0.states:y'] = phase.interpolate(ys=[10, 5], nodes='state_input') - p['traj0.phase0.states:v'] = phase.interpolate(ys=[0, 9.9], nodes='state_input') - p['traj0.phase0.controls:theta'] = phase.interpolate(ys=[5, 100], nodes='control_input') - p['traj0.phase0.parameters:g'] = 9.80665 - - dm.run_problem(p, run_driver=run_driver, simulate=True) - - return p - - def run_asserts(self, p): - - sol_db = 'dymos_solution.db' - sim_db = 'dymos_simulation.db' - if om_version()[0] > (3, 34, 2): - sol_db = p.get_outputs_dir() / sol_db - sim_db = p.model.traj0.sim_prob.get_outputs_dir() / sim_db - - for db in (sol_db, sim_db): - p = om.CaseReader(db).get_case('final') - - t_initial = p.get_val('traj0.phase0.timeseries.time')[0] - tf = p.get_val('traj0.phase0.timeseries.time')[-1] - - x0 = p.get_val('traj0.phase0.timeseries.x')[0] - xf = p.get_val('traj0.phase0.timeseries.x')[-1] - - y0 = p.get_val('traj0.phase0.timeseries.y')[0] - yf = p.get_val('traj0.phase0.timeseries.y')[-1] - - v0 = p.get_val('traj0.phase0.timeseries.v')[0] - vf = p.get_val('traj0.phase0.timeseries.v')[-1] - - g = p.get_val('traj0.phase0.parameter_vals:g')[0] - - thetaf = p.get_val('traj0.phase0.timeseries.theta')[-1] - - assert_near_equal(t_initial, 0.0) - assert_near_equal(x0, 0.0) - assert_near_equal(y0, 10.0) - assert_near_equal(v0, 0.0) - - assert_near_equal(tf, 1.8016, tolerance=0.01) - assert_near_equal(xf, 10.0, tolerance=0.01) - assert_near_equal(yf, 5.0, tolerance=0.01) - assert_near_equal(vf, 9.902, tolerance=0.01) - assert_near_equal(g, 9.80665, tolerance=0.01) - - assert_near_equal(thetaf, 100.12, tolerance=0.01) - - def test_ex_brachistochrone_radau_uncompressed(self): - p = self._make_problem(transcription='radau-ps', compressed=False) - self.run_asserts(p) - - def test_ex_brachistochrone_gl_uncompressed(self): - p = self._make_problem(transcription='gauss-lobatto', compressed=False) - self.run_asserts(p) diff --git a/dymos/transcriptions/common/test/test_control_interp_comp.py b/dymos/transcriptions/common/test/test_control_interp_comp.py index 147e8391a..9b19588f3 100644 --- a/dymos/transcriptions/common/test/test_control_interp_comp.py +++ b/dymos/transcriptions/common/test/test_control_interp_comp.py @@ -164,7 +164,7 @@ def test_control_interp_scalar(self): np.atleast_2d(b_rate2_expected).T) np.set_printoptions(linewidth=1024) - cpd = p.check_partials(compact_print=False, method='cs') + cpd = p.check_partials(compact_print=False, method='cs', out_stream=None) assert_check_partials(cpd) def test_control_interp_vector(self, transcription='gauss-lobatto', compressed=True): @@ -263,7 +263,7 @@ def test_control_interp_vector(self, transcription='gauss-lobatto', compressed=T a2_rate2_expected) np.set_printoptions(linewidth=1024) - cpd = p.check_partials(compact_print=True, method='cs') + cpd = p.check_partials(compact_print=True, method='cs', out_stream=None) assert_check_partials(cpd) @@ -362,7 +362,7 @@ def test_control_interp_matrix_3x1(self, transcription='gauss-lobatto', compress a2_rate2_expected) np.set_printoptions(linewidth=1024) - cpd = p.check_partials(compact_print=True, method='cs') + cpd = p.check_partials(compact_print=True, method='cs', out_stream=None) assert_check_partials(cpd) @@ -474,7 +474,7 @@ def test_control_interp_matrix_2x2(self, transcription='gauss-lobatto', compress a3_rate2_expected) with np.printoptions(linewidth=100000, edgeitems=100000): - cpd = p.check_partials(compact_print=True, method='cs') + cpd = p.check_partials(compact_print=True, method='cs', out_stream=None) assert_check_partials(cpd) diff --git a/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py b/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py index 73954efa7..fad5d62c3 100644 --- a/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py +++ b/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py @@ -111,7 +111,7 @@ def test_1_state_run_model(self): prob.run_model() with np.printoptions(linewidth=1024): - cpd = prob.check_partials(compact_print=True) + cpd = prob.check_partials(compact_print=True, out_stream=None) assert_check_partials(cpd, rtol=1.0E-5) @@ -157,7 +157,7 @@ def test_2_states_run_model(self): assert_near_equal(y_f[-1, ...], 0.1691691, tolerance=1.0E-5) with np.printoptions(linewidth=1024): - cpd = prob.check_partials(compact_print=True, method='fd') + cpd = prob.check_partials(compact_print=True, method='fd', out_stream=None) assert_check_partials(cpd, atol=1.0E-5, rtol=1.0E-5) dymos_options['include_check_partials'] = False @@ -226,10 +226,14 @@ def test_brachistochrone_explicit_shooting(self): assert_near_equal(t[-1, ...], 1.8016, tolerance=1.0E-2) with np.printoptions(linewidth=1024): - cpd = prob.check_partials(compact_print=True, method='cs', excludes=['traj0.phases.phase0.integrator']) + cpd = prob.check_partials(compact_print=True, method='cs', + excludes=['traj0.phases.phase0.integrator'], + out_stream=None) assert_check_partials(cpd, atol=1.0E-5, rtol=1.0E-5) - cpd = prob.check_partials(compact_print=True, method='fd', includes=['traj0.phases.phase0.integrator']) + cpd = prob.check_partials(compact_print=True, method='fd', + includes=['traj0.phases.phase0.integrator'], + out_stream=None) assert_check_partials(cpd, atol=1.0E-5, rtol=1.0E-5) dymos_options['include_check_partials'] = False @@ -373,11 +377,11 @@ def test_brachistochrone_explicit_shooting_path_constraint_polynomial_control(se assert_near_equal(t[-1, ...], 1.807379, tolerance=1.0E-5) with np.printoptions(linewidth=1024): - cpd = prob.check_partials(compact_print=True, method='cs', + cpd = prob.check_partials(compact_print=True, method='cs', out_stream=None, excludes=['traj0.phases.phase0.integrator']) assert_check_partials(cpd, atol=1.0E-5, rtol=1.0E-5) - cpd = prob.check_partials(compact_print=True, method='fd', + cpd = prob.check_partials(compact_print=True, method='fd', out_stream=None, includes=['traj0.phases.phase0.integrator']) assert_check_partials(cpd, atol=1.0E-5, rtol=1.0E-5) diff --git a/dymos/transcriptions/pseudospectral/components/test/test_collocation_defect_sol_opt.py b/dymos/transcriptions/pseudospectral/components/test/test_collocation_defect_sol_opt.py index aac1cafa4..50c88684e 100644 --- a/dymos/transcriptions/pseudospectral/components/test/test_collocation_defect_sol_opt.py +++ b/dymos/transcriptions/pseudospectral/components/test/test_collocation_defect_sol_opt.py @@ -185,22 +185,22 @@ def test_partials(self): np.set_printoptions(linewidth=1024, edgeitems=1e1000) p = self.make_prob('radau-ps', n_segs=2, order=5, compressed=False) - cpd = p.check_partials(compact_print=True, method='fd') + cpd = p.check_partials(compact_print=True, method='fd', out_stream=None) del cpd['state_indep'] assert_check_partials(cpd) p = self.make_prob('radau-ps', n_segs=2, order=5, compressed=True) - cpd = p.check_partials(compact_print=True, method='fd') + cpd = p.check_partials(compact_print=True, method='fd', out_stream=None) del cpd['state_indep'] assert_check_partials(cpd) p = self.make_prob('gauss-lobatto', n_segs=3, order=5, compressed=False) - cpd = p.check_partials(compact_print=True, method='fd') + cpd = p.check_partials(compact_print=True, method='fd', out_stream=None) del cpd['state_indep'] assert_check_partials(cpd) p = self.make_prob('gauss-lobatto', n_segs=4, order=3, compressed=True) - cpd = p.check_partials(compact_print=True, method='fd') + cpd = p.check_partials(compact_print=True, method='fd', out_stream=None) del cpd['state_indep'] assert_check_partials(cpd) diff --git a/dymos/transcriptions/pseudospectral/components/test/test_control_endpoint_defect_comp.py b/dymos/transcriptions/pseudospectral/components/test/test_control_endpoint_defect_comp.py index 87db3c314..8c710099f 100644 --- a/dymos/transcriptions/pseudospectral/components/test/test_control_endpoint_defect_comp.py +++ b/dymos/transcriptions/pseudospectral/components/test/test_control_endpoint_defect_comp.py @@ -62,7 +62,7 @@ def test_results(self): tolerance=1.0E-9) def test_partials(self): - cpd = self.p.check_partials(compact_print=False, method='cs') + cpd = self.p.check_partials(compact_print=False, method='cs', out_stream=None) assert_check_partials(cpd)