From 329dd40649631d0c68fadae850f701f1e3ff467c Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 11 May 2026 13:46:25 -0400 Subject: [PATCH 1/2] Move unique doc tests to test directories and remove redundant doc tests - Move robertson_problem ExplicitShooting/LSODA test to test/ (unique stiff ODE handling) - Move cannonball boundary_balance/implicit duration test to test/ (unique pattern) - Move length_constrained_brachistochrone arc length constraint test to test/ (unique constraint pattern) - Remove 12 redundant doc tests that duplicate notebook examples and existing test/ coverage All tests verified to pass. Reduces maintenance burden and CI time. --- .../doc/test_doc_balanced_field_length.py | 1075 ----------------- .../doc/brachistochrone_ode.py | 73 -- .../doc/test_doc_brachistochrone.py | 403 ------ ...test_doc_brachistochrone_static_gravity.py | 177 --- .../test_doc_brachistochrone_tandem_phases.py | 163 --- ...st_doc_brachistochrone_upstream_control.py | 149 --- ...test_doc_brachistochrone_upstream_state.py | 167 --- dymos/examples/cannonball/doc/__init__.py | 0 .../doc/test_doc_two_phase_cannonball.py | 256 ---- .../test_cannonball_implicit_duration.py} | 0 .../double_integrator/doc/__init__.py | 0 .../doc/test_doc_double_integrator.py | 111 -- .../finite_burn_orbit_raise/doc/__init__.py | 0 .../doc/test_doc_finite_burn_orbit_raise.py | 241 ---- .../examples/hyper_sensitive/doc/__init__.py | 0 .../doc/test_doc_hyper_sensitive.py | 121 -- .../doc/__init__.py | 0 .../test}/__init__.py | 0 ...est_length_constrained_brachistochrone.py} | 0 .../robertson_problem/doc/__init__.py | 0 .../test}/__init__.py | 0 .../test_robertson_problem.py} | 4 +- dymos/examples/robot_arm/doc/__init__.py | 0 .../robot_arm/doc/test_doc_robot_arm.py | 114 -- .../examples/shuttle_reentry/doc/__init__.py | 0 .../shuttle_reentry/doc/test_doc_reentry.py | 110 -- dymos/examples/ssto/doc/__init__.py | 0 .../examples/ssto/doc/test_doc_ssto_earth.py | 142 --- .../test_doc_ssto_linear_tangent_guidance.py | 273 ----- .../doc/test_doc_ssto_polynomial_control.py | 317 ----- dymos/examples/vanderpol/doc/__init__.py | 0 .../vanderpol/doc/test_doc_vanderpol.py | 132 -- 32 files changed, 2 insertions(+), 4026 deletions(-) delete mode 100644 dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py delete mode 100644 dymos/examples/brachistochrone/doc/brachistochrone_ode.py delete mode 100644 dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py delete mode 100644 dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py delete mode 100644 dymos/examples/brachistochrone/doc/test_doc_brachistochrone_tandem_phases.py delete mode 100644 dymos/examples/brachistochrone/doc/test_doc_brachistochrone_upstream_control.py delete mode 100644 dymos/examples/brachistochrone/doc/test_doc_brachistochrone_upstream_state.py delete mode 100644 dymos/examples/cannonball/doc/__init__.py delete mode 100644 dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py rename dymos/examples/cannonball/{doc/test_doc_cannonball_implicit_duration.py => test/test_cannonball_implicit_duration.py} (100%) delete mode 100644 dymos/examples/double_integrator/doc/__init__.py delete mode 100644 dymos/examples/double_integrator/doc/test_doc_double_integrator.py delete mode 100644 dymos/examples/finite_burn_orbit_raise/doc/__init__.py delete mode 100644 dymos/examples/finite_burn_orbit_raise/doc/test_doc_finite_burn_orbit_raise.py delete mode 100644 dymos/examples/hyper_sensitive/doc/__init__.py delete mode 100644 dymos/examples/hyper_sensitive/doc/test_doc_hyper_sensitive.py delete mode 100644 dymos/examples/length_constrained_brachistochrone/doc/__init__.py rename dymos/examples/{balanced_field/doc => length_constrained_brachistochrone/test}/__init__.py (100%) rename dymos/examples/length_constrained_brachistochrone/{doc/test_doc_length_constrained_brachistochrone.py => test/test_length_constrained_brachistochrone.py} (100%) delete mode 100644 dymos/examples/robertson_problem/doc/__init__.py rename dymos/examples/{brachistochrone/doc => robertson_problem/test}/__init__.py (100%) rename dymos/examples/robertson_problem/{doc/test_doc_robertson_problem.py => test/test_robertson_problem.py} (97%) delete mode 100644 dymos/examples/robot_arm/doc/__init__.py delete mode 100644 dymos/examples/robot_arm/doc/test_doc_robot_arm.py delete mode 100644 dymos/examples/shuttle_reentry/doc/__init__.py delete mode 100644 dymos/examples/shuttle_reentry/doc/test_doc_reentry.py delete mode 100644 dymos/examples/ssto/doc/__init__.py delete mode 100644 dymos/examples/ssto/doc/test_doc_ssto_earth.py delete mode 100644 dymos/examples/ssto/doc/test_doc_ssto_linear_tangent_guidance.py delete mode 100644 dymos/examples/ssto/doc/test_doc_ssto_polynomial_control.py delete mode 100644 dymos/examples/vanderpol/doc/__init__.py delete mode 100644 dymos/examples/vanderpol/doc/test_doc_vanderpol.py diff --git a/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py b/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py deleted file mode 100644 index d815cdcaca..0000000000 --- a/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py +++ /dev/null @@ -1,1075 +0,0 @@ -import unittest - -import numpy as np -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse -from openmdao.utils.assert_utils import assert_near_equal - -import openmdao.api as om -from openmdao.utils.general_utils import set_pyoptsparse_opt -import dymos as dm -from dymos.examples.balanced_field.balanced_field_ode import BalancedFieldODEComp - - -regression_data = {} -regression_data[dm.Radau] = { - "constraints": { - "traj.linkages.br_to_v1:time_final|v1_to_vr:time_initial": [0.0], - "traj.linkages.br_to_v1:r_final|v1_to_vr:r_initial": [0.0], - "traj.linkages.br_to_v1:v_final|v1_to_vr:v_initial": [0.0], - "traj.linkages.v1_to_vr:time_final|rotate:time_initial": [0.0], - "traj.linkages.v1_to_vr:r_final|rotate:r_initial": [-1450.0], - "traj.linkages.v1_to_vr:v_final|rotate:v_initial": [30.0], - "traj.linkages.v1_to_vr:alpha_final|rotate:alpha_initial": [0.0], - "traj.linkages.rotate:time_final|climb:time_initial": [0.0], - "traj.linkages.rotate:r_final|climb:r_initial": [276.0], - "traj.linkages.rotate:v_final|climb:v_initial": [2.688888888888883], - "traj.linkages.rotate:alpha_final|climb:alpha_initial": [-5.0], - "traj.linkages.br_to_v1:time_final|rto:time_initial": [0.0], - "traj.linkages.br_to_v1:r_final|rto:r_initial": [0.0], - "traj.linkages.br_to_v1:v_final|rto:v_initial": [-10.0], - "traj.linkages.rto:r_final|climb:r_final": [3.3236], - "traj.br_to_v1.collocation_constraint.defects:r": [ - 0.4166666666666665, - 0.3476289672207841, - 0.2523710327792159, - 0.22222222222222143, - 0.1531845227763402, - 0.05792658833477156, - 0.027777777777776274, - -0.04125992166810405, - -0.1365178561096728, - ], - "traj.br_to_v1.collocation_constraint.defects:v": [ - 0.00649795093524439, - 0.006616445486851315, - 0.007169036357321012, - 0.007437927094391042, - 0.0082239006449162, - 0.009697476299502095, - 0.010257855571831287, - 0.01171130812127473, - 0.01410586855997708, - ], - "traj.climb.h[final]": [1.0], - "traj.climb.gam[final]": [1.0], - "traj.climb.v_over_v_stall[final]": [0.9823389802049667], - "traj.climb.gam[path]": [ - 0.0, - 0.07101020514433644, - 0.1689897948556636, - 0.20000000000000004, - 0.20000000000000004, - 0.27101020514433644, - 0.36898979485566363, - 0.4000000000000001, - 0.4000000000000001, - 0.4710102051443366, - 0.5689897948556638, - 0.6000000000000002, - 0.6000000000000002, - 0.6710102051443366, - 0.7689897948556637, - 0.8000000000000002, - 0.8000000000000002, - 0.8710102051443366, - 0.9689897948556637, - ], - "traj.climb.h[path]": [ - 0.0, - 2.4853571800517744, - 5.914642819948225, - 6.999999999999999, - 6.999999999999999, - 9.485357180051773, - 12.914642819948225, - 13.999999999999998, - 13.999999999999998, - 16.48535718005178, - 19.91464281994823, - 21.0, - 21.0, - 23.48535718005178, - 26.914642819948227, - 27.999999999999996, - 27.999999999999996, - 30.485357180051775, - 33.91464281994823, - ], - "traj.climb.collocation_constraint.defects:r": [ - -0.10822666666666667, - -0.10877224764873458, - -0.10951713737628663, - -0.1097509603715007, - -0.11028284972038041, - -0.11100855808579062, - -0.11123624046202378, - -0.11175394240945394, - -0.11245978793155431, - -0.11268111463694755, - -0.11318413790684131, - -0.11386944544945637, - -0.1140842034756845, - -0.1145720615353963, - -0.11523616265315353, - ], - "traj.climb.collocation_constraint.defects:h": [ - 1.0667999999999995, - 0.29830989982523587, - -0.7731445565941426, - -1.1149253287248142, - -1.902381792745922, - -2.9998640073265532, - -3.349847639357176, - -4.156023740270561, - -5.279186548509148, - -5.637261541533632, - -6.461899026513403, - -7.610379407812043, - -7.976429256188403, - -8.81925842906678, - -9.992677623677165, - ], - "traj.climb.collocation_constraint.defects:v": [ - -0.01461939424215987, - -0.013680311182054173, - -0.012385122427714935, - -0.01197529115802313, - -0.01103693473646186, - -0.009742414304801463, - -0.0093327538694019, - -0.008394772109415277, - -0.007100805084518509, - -0.006691339937664971, - -0.005753853875203694, - -0.00446069768091129, - -0.004051522420108261, - -0.0031147634445885157, - -0.0018227595387125509, - ], - "traj.climb.collocation_constraint.defects:gam": [ - 0.716834308339309, - 0.6879445904659692, - 0.6480853509270524, - 0.6354706790667543, - 0.6065855235688165, - 0.566732541890516, - 0.5541198471763983, - 0.5252392227663774, - 0.48539252338321354, - 0.4727818295474824, - 0.4439058189556195, - 0.4040655795810557, - 0.3914569579135891, - 0.3625857510819955, - 0.32275229358217394, - ], - "traj.climb.continuity_comp.defect_control_rates:alpha_rate": [ - 2.220446049250313e-15, - 2.220446049250313e-15, - 2.2204460492503127e-15, - 2.220446049250314e-15, - ], - "traj.climb.continuity_comp.defect_controls:alpha": [0.0, 0.0, 0.0, 0.0], - "traj.rotate.F_r[final]": [4.989619710033791], - "traj.rotate.collocation_constraint.defects:r": [ - -0.05833333333333443, - -0.058826459757946446, - -0.059506873575385956, - -0.05972222222222309, - -0.0602153486468354, - -0.060895762464275166, - -0.06111111111111143, - -0.06160423753572436, - -0.06228465135316393, - ], - "traj.rotate.collocation_constraint.defects:v": [ - -0.0011079918722422847, - -0.001096507055634512, - -0.00108052139116765, - -0.0010754284124435546, - -0.001063705210478725, - -0.0010473906228746822, - -0.0010421935411017993, - -0.0010302319537806714, - -0.0010135884430394586, - ], - "traj.rto.v[final]": [0.0], - "traj.rto.collocation_constraint.defects:r": [ - -0.22500000000000214, - -0.1490585306095287, - -0.04427480272380285, - -0.011111111111113143, - 0.06483035827935983, - 0.1696140861650864, - 0.2027777777777753, - 0.2787192471682481, - 0.3835029750539747, - ], - "traj.rto.collocation_constraint.defects:v": [ - -0.09358803170582275, - -0.07535584468522619, - -0.05396466020985954, - -0.048104102054440434, - -0.03633156138447398, - -0.023853370440509947, - -0.020813744263609962, - -0.01550084994427481, - -0.011935652531713638, - ], - "traj.v1_to_vr.v_over_v_stall[final]": [0.015444600279553376], - "traj.v1_to_vr.collocation_constraint.defects:r": [ - -0.950000000000002, - -0.9569037699445878, - -0.9664295633887448, - -0.969444444444446, - -0.9763482143890325, - -0.9858740078331893, - -0.9888888888888886, - -0.9957926588334773, - -1.0053184522776337, - ], - "traj.v1_to_vr.collocation_constraint.defects:v": [ - -0.04637708701672346, - -0.04617565837153162, - -0.045893837027591915, - -0.04580370155964456, - -0.04559559812446296, - -0.04530456693268199, - -0.04521151657938215, - -0.044996738354211245, - -0.044696497314589134, - ], - }, - "design_vars": { - "traj.br_to_v1.t_duration": [3.5], - "traj.br_to_v1.states:r": [ - 0.29587585476806844, - 0.7041241452319315, - 0.8333333333333333, - 1.1292091881014017, - 1.5374574785652648, - 1.6666666666666665, - 1.962542521434735, - 2.3707908118985985, - 2.5, - ], - "traj.br_to_v1.states:v": [ - 0.11835034190722737, - 0.2816496580927726, - 0.3333333333333333, - 0.4516836752405607, - 0.6149829914261059, - 0.6666666666666666, - 0.785017008573894, - 0.9483163247594394, - 1.0, - ], - "traj.climb.t_initial": [75.0], - "traj.climb.t_duration": [15.0], - "traj.climb.states:r": [ - 1.524, - 1.5348219552639972, - 1.5497540447360032, - 1.55448, - 1.565301955263997, - 1.5802340447360035, - 1.5849600000000001, - 1.5957819552639971, - 1.6107140447360035, - 1.61544, - 1.6262619552639968, - 1.6411940447360034, - 1.64592, - 1.6567419552639973, - 1.6716740447360035, - 1.6764000000000001, - ], - "traj.climb.states:h": [ - 0.757536868479781, - 1.802783131520219, - 2.1336, - 2.891136868479781, - 3.9363831315202193, - 4.2672, - 5.024736868479782, - 6.069983131520221, - 6.400800000000001, - 7.158336868479783, - 8.20358313152022, - 8.5344, - 9.291936868479782, - 10.33718313152022, - 10.668000000000001, - ], - "traj.climb.states:v": [ - 0.8231111111111112, - 0.8267641916646475, - 0.8318046972242413, - 0.8334, - 0.8370530805535364, - 0.8420935861131303, - 0.8436888888888889, - 0.8473419694424255, - 0.8523824750020192, - 0.853977777777778, - 0.8576308583313142, - 0.8626713638909083, - 0.8642666666666667, - 0.8679197472202032, - 0.872960252779797, - 0.8745555555555556, - ], - "traj.climb.states:gam": [ - 0.12393618822852859, - 0.2949428322501105, - 0.3490658503988659, - 0.4730020386273945, - 0.6440086826489765, - 0.6981317007977318, - 0.8220678890262606, - 0.9930745330478425, - 1.0471975511965979, - 1.1711337394251264, - 1.3421403834467083, - 1.3962634015954636, - 1.5201995898239922, - 1.6912062338455742, - 1.7453292519943295, - ], - "traj.climb.controls:alpha": [ - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - 0.5, - ], - "traj.rotate.t_initial": [70.0], - "traj.rotate.t_duration": [5.0], - "traj.rotate.states:r": [ - 1.75, - 1.7559175170953614, - 1.7640824829046386, - 1.7666666666666668, - 1.772584183762028, - 1.7807491495713053, - 1.7833333333333332, - 1.7892508504286948, - 1.797415816237972, - 1.8, - ], - "traj.rotate.states:v": [ - 0.8, - 0.8059175170953614, - 0.8140824829046388, - 0.8166666666666668, - 0.822584183762028, - 0.8307491495713053, - 0.8333333333333333, - 0.8392508504286947, - 0.847415816237972, - 0.85, - ], - "traj.rotate.controls:alpha": [0.0, 0.0], - "traj.rto.t_initial": [35.0], - "traj.rto.t_duration": [35.0], - "traj.rto.states:r": [ - 2.5, - 2.7958758547680684, - 3.2041241452319316, - 3.3333333333333335, - 3.629209188101402, - 4.037457478565265, - 4.166666666666667, - 4.462542521434735, - 4.870790811898598, - 5.0, - ], - "traj.rto.states:v": [ - 1.1, - 0.96981462390205, - 0.7901853760979503, - 0.7333333333333334, - 0.6031479572353832, - 0.4235187094312835, - 0.3666666666666667, - 0.2364812905687166, - 0.05685204276461679, - 0.0, - ], - "traj.v1_to_vr.t_initial": [35.0], - "traj.v1_to_vr.t_duration": [35.0], - "traj.v1_to_vr.states:r": [ - 2.5, - 2.2396292478040998, - 1.8803707521959006, - 1.7666666666666668, - 1.5062959144707664, - 1.147037418862567, - 1.0333333333333332, - 0.7729625811374332, - 0.4137040855292336, - 0.3, - ], - "traj.v1_to_vr.states:v": [ - 1.0, - 1.0118350341907227, - 1.0281649658092773, - 1.0333333333333334, - 1.0451683675240562, - 1.0614982991426107, - 1.0666666666666667, - 1.0785017008573894, - 1.094831632475944, - 1.1, - ], - }, - "objectives": {"traj.rto.states:r": [5000.0]}, -} - -regression_data[dm.GaussLobatto] = { - "constraints": { - "traj.linkages.br_to_v1:time_final|v1_to_vr:time_initial": [0.0], - "traj.linkages.br_to_v1:r_final|v1_to_vr:r_initial": [0.0], - "traj.linkages.br_to_v1:v_final|v1_to_vr:v_initial": [0.0], - "traj.linkages.v1_to_vr:time_final|rotate:time_initial": [0.0], - "traj.linkages.v1_to_vr:r_final|rotate:r_initial": [-1450.0], - "traj.linkages.v1_to_vr:v_final|rotate:v_initial": [30.0], - "traj.linkages.v1_to_vr:alpha_final|rotate:alpha_initial": [0.0], - "traj.linkages.rotate:time_final|climb:time_initial": [0.0], - "traj.linkages.rotate:r_final|climb:r_initial": [276.0], - "traj.linkages.rotate:v_final|climb:v_initial": [2.688888888888883], - "traj.linkages.rotate:alpha_final|climb:alpha_initial": [-5.0], - "traj.linkages.br_to_v1:time_final|rto:time_initial": [0.0], - "traj.linkages.br_to_v1:r_final|rto:r_initial": [0.0], - "traj.linkages.br_to_v1:v_final|rto:v_initial": [-10.0], - "traj.linkages.rto:r_final|climb:r_final": [3.3236], - "traj.br_to_v1.collocation_constraint.defects:r": [ - 0.47902958681012453, - 0.1870887604303732, - -0.10485206594937782, - ], - "traj.br_to_v1.collocation_constraint.defects:v": [ - 0.010217577615992882, - 0.013042811161855497, - 0.018693278253580733, - ], - "traj.climb.h[final]": [1.0], - "traj.climb.gam[final]": [1.0], - "traj.climb.v_over_v_stall[final]": [0.9823389802049667], - "traj.climb.gam[path]": [ - 0.0, - 0.0883455185920389, - 0.20000000000000004, - 0.20000000000000004, - 0.2883473516820108, - 0.4000000000000001, - 0.4000000000000001, - 0.4883491871897563, - 0.6000000000000002, - 0.6000000000000002, - 0.688351070214837, - 0.8000000000000002, - 0.8000000000000002, - 0.8883530431698609, - ], - "traj.climb.h[path]": [ - 0.0, - 1.7105271253897536, - 6.999999999999999, - 6.999999999999999, - 8.666894430255608, - 13.999999999999998, - 13.999999999999998, - 15.623840303332958, - 21.0, - 21.0, - 22.581391310158487, - 27.999999999999996, - 27.999999999999996, - 29.539573821693008, - ], - "traj.climb.collocation_constraint.defects:r": [ - -0.1635881987944115, - -0.1658477115826216, - -0.16804769765480093, - -0.17018604309853758, - -0.17226068025412108, - ], - "traj.climb.collocation_constraint.defects:h": [ - 0.09618576144272717, - -3.2166033112704633, - -6.608701632651881, - -10.079026976699794, - -13.626448340964746, - ], - "traj.climb.collocation_constraint.defects:v": [ - -0.020091291412078167, - -0.016125860961723488, - -0.01216255501671909, - -0.008201285233091979, - -0.00424289445112884, - ], - "traj.climb.collocation_constraint.defects:gam": [ - 1.0090664482860263, - 0.8871644579785876, - 0.7652774854957041, - 0.643407653401527, - 0.5215561354920641, - ], - "traj.climb.continuity_comp.defect_control_rates:alpha_rate": [0.0, 0.0, 0.0, 0.0], - "traj.rotate.F_r[final]": [4.989619710033791], - "traj.rotate.collocation_constraint.defects:r": [ - -0.08854234507207909, - -0.090625692393153, - -0.09270903971422675, - ], - "traj.rotate.collocation_constraint.defects:v": [ - -0.0016376332341922898, - -0.0015882838231674338, - -0.0015379272813052295, - ], - "traj.rto.v[final]": [0.0], - "traj.rto.collocation_constraint.defects:r": [ - -0.18371640640749326, - 0.13977015615550417, - 0.4632567187185013, - ], - "traj.rto.collocation_constraint.defects:v": [ - -0.1054141896634857, - -0.04992513115327707, - -0.022180601898172673, - ], - "traj.v1_to_vr.v_over_v_stall[final]": [0.015444600279553376], - "traj.v1_to_vr.collocation_constraint.defects:r": [ - -1.439666952045824, - -1.4688363603096213, - -1.4980057685734194, - ], - "traj.v1_to_vr.collocation_constraint.defects:v": [ - -0.06913547541722186, - -0.06826113323636616, - -0.06735858646903126, - ], - }, - "design_vars": { - "traj.br_to_v1.t_duration": [3.5], - "traj.br_to_v1.states:r": [0.8333333333333333, 1.6666666666666665, 2.5], - "traj.br_to_v1.states:v": [0.3333333333333333, 0.6666666666666666, 1.0], - "traj.climb.t_initial": [75.0], - "traj.climb.t_duration": [15.0], - "traj.climb.states:r": [ - 1.524, - 1.55448, - 1.5849600000000001, - 1.61544, - 1.64592, - 1.6764000000000001, - ], - "traj.climb.states:h": [2.1336, 4.2672, 6.400800000000001, 8.5344, 10.668000000000001], - "traj.climb.states:v": [ - 0.8231111111111112, - 0.8334, - 0.8436888888888889, - 0.853977777777778, - 0.8642666666666667, - 0.8745555555555556, - ], - "traj.climb.states:gam": [ - 0.3490658503988659, - 0.6981317007977318, - 1.0471975511965979, - 1.3962634015954636, - 1.7453292519943295, - ], - "traj.climb.controls:alpha": [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5], - "traj.rotate.t_initial": [70.0], - "traj.rotate.t_duration": [5.0], - "traj.rotate.states:r": [1.75, 1.7666666666666668, 1.7833333333333332, 1.8], - "traj.rotate.states:v": [0.8, 0.8166666666666668, 0.8333333333333333, 0.85], - "traj.rotate.controls:alpha": [0.0, 0.0], - "traj.rto.t_initial": [35.0], - "traj.rto.t_duration": [35.0], - "traj.rto.states:r": [2.5, 3.3333333333333335, 4.166666666666667, 5.0], - "traj.rto.states:v": [1.1, 0.7333333333333334, 0.3666666666666667, 0.0], - "traj.v1_to_vr.t_initial": [35.0], - "traj.v1_to_vr.t_duration": [35.0], - "traj.v1_to_vr.states:r": [2.5, 1.7666666666666668, 1.0333333333333332, 0.3], - "traj.v1_to_vr.states:v": [1.0, 1.0333333333333334, 1.0666666666666667, 1.1], - }, - "objectives": {"traj.rto.states:r": [5000.0]}, -} - - -@use_tempdirs -class TestBalancedFieldLengthForDocs(unittest.TestCase): - def _make_problem(self, tx, optimizer='IPOPT'): - p = om.Problem() - - # Use IPOPT if available, with fallback to SLSQP - if optimizer is not None: - p.driver = om.pyOptSparseDriver(optimizer="IPOPT") - p.driver.options["optimizer"] = optimizer - p.driver.declare_coloring() - p.driver.options["print_results"] = False - if optimizer == "IPOPT": - p.driver.opt_settings["print_level"] = 0 - - # First Phase: Brake release to V1 - both engines operable - br_to_v1 = dm.Phase( - ode_class=BalancedFieldODEComp, - transcription=tx(num_segments=3), - ode_init_kwargs={"mode": "runway"}, - ) - br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0) - br_to_v1.add_state("r", fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0) - br_to_v1.add_state("v", fix_initial=True, lower=0, ref=100.0, defect_ref=100.0) - br_to_v1.add_parameter("alpha", val=0.0, opt=False, units="deg") - br_to_v1.add_timeseries_output("*") - - # Second Phase: Rejected takeoff at V1 - no engines operable - rto = dm.Phase( - ode_class=BalancedFieldODEComp, - transcription=tx(num_segments=3), - ode_init_kwargs={"mode": "runway"}, - ) - rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - rto.add_state("r", fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rto.add_state("v", fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - rto.add_parameter("alpha", val=0.0, opt=False, units="deg") - rto.add_timeseries_output("*") - - # Third Phase: V1 to Vr - single engine operable - v1_to_vr = dm.Phase( - ode_class=BalancedFieldODEComp, - transcription=tx(num_segments=3), - ode_init_kwargs={"mode": "runway"}, - ) - v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0) - v1_to_vr.add_state("r", fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - v1_to_vr.add_state("v", fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - v1_to_vr.add_parameter("alpha", val=0.0, opt=False, units="deg") - v1_to_vr.add_timeseries_output("*") - - # Fourth Phase: Rotate - single engine operable - rotate = dm.Phase( - ode_class=BalancedFieldODEComp, - transcription=tx(num_segments=3), - ode_init_kwargs={"mode": "runway"}, - ) - rotate.set_time_options(fix_initial=False, duration_bounds=(1.0, 5), duration_ref=1.0) - rotate.add_state("r", fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - rotate.add_state("v", fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - rotate.add_control( - "alpha", - order=1, - opt=True, - units="deg", - lower=0, - upper=10, - ref=10, - val=[0, 10], - control_type="polynomial", - ) - rotate.add_timeseries_output("*") - - # Fifth Phase: Climb to target speed and altitude at end of runway. - climb = dm.Phase( - ode_class=BalancedFieldODEComp, - transcription=tx(num_segments=5), - ode_init_kwargs={"mode": "climb"}, - ) - climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0) - climb.add_state("r", fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0) - climb.add_state("h", fix_initial=True, lower=0, ref=1.0, defect_ref=1.0) - climb.add_state("v", fix_initial=False, lower=0, ref=100.0, defect_ref=100.0) - climb.add_state("gam", fix_initial=True, lower=0, ref=0.05, defect_ref=0.05) - climb.add_control("alpha", opt=True, units="deg", lower=-10, upper=15, ref=10) - climb.add_timeseries_output("*") - - # Instantiate the trajectory and add phases - traj = dm.Trajectory() - p.model.add_subsystem("traj", traj) - traj.add_phase("br_to_v1", br_to_v1) - traj.add_phase("rto", rto) - traj.add_phase("v1_to_vr", v1_to_vr) - traj.add_phase("rotate", rotate) - traj.add_phase("climb", climb) - - # Add parameters common to multiple phases to the trajectory - traj.add_parameter( - "m", - val=174200.0, - opt=False, - units="lbm", - desc="aircraft mass", - targets={ - "br_to_v1": ["m"], - "v1_to_vr": ["m"], - "rto": ["m"], - "rotate": ["m"], - "climb": ["m"], - }, - ) - - traj.add_parameter( - "T_nominal", - val=27000 * 2, - opt=False, - units="lbf", - static_target=True, - desc="nominal aircraft thrust", - targets={"br_to_v1": ["T"]}, - ) - - traj.add_parameter( - "T_engine_out", - val=27000, - opt=False, - units="lbf", - static_target=True, - desc="thrust under a single engine", - targets={"v1_to_vr": ["T"], "rotate": ["T"], "climb": ["T"]}, - ) - - traj.add_parameter( - "T_shutdown", - val=0.0, - opt=False, - units="lbf", - static_target=True, - desc="thrust when engines are shut down for rejected takeoff", - targets={"rto": ["T"]}, - ) - - traj.add_parameter( - "mu_r_nominal", - val=0.03, - opt=False, - units=None, - static_target=True, - desc="nominal runway friction coefficient", - targets={"br_to_v1": ["mu_r"], "v1_to_vr": ["mu_r"], "rotate": ["mu_r"]}, - ) - - traj.add_parameter( - "mu_r_braking", - val=0.3, - opt=False, - units=None, - static_target=True, - desc="runway friction coefficient under braking", - targets={"rto": ["mu_r"]}, - ) - - traj.add_parameter( - "h_runway", - val=0.0, - opt=False, - units="ft", - desc="runway altitude", - targets={"br_to_v1": ["h"], "v1_to_vr": ["h"], "rto": ["h"], "rotate": ["h"]}, - ) - - traj.add_parameter( - "rho", - val=1.225, - opt=False, - units="kg/m**3", - static_target=True, - desc="atmospheric density", - targets={"br_to_v1": ["rho"], "v1_to_vr": ["rho"], "rto": ["rho"], "rotate": ["rho"]}, - ) - - traj.add_parameter( - "S", - val=124.7, - opt=False, - units="m**2", - static_target=True, - desc="aerodynamic reference area", - targets={ - "br_to_v1": ["S"], - "v1_to_vr": ["S"], - "rto": ["S"], - "rotate": ["S"], - "climb": ["S"], - }, - ) - - traj.add_parameter( - "CD0", - val=0.03, - opt=False, - units=None, - static_target=True, - desc="zero-lift drag coefficient", - targets={ - f"{phase}": ["CD0"] for phase in ["br_to_v1", "v1_to_vr", "rto", "rotate", "climb"] - }, - ) - - traj.add_parameter( - "AR", - val=9.45, - opt=False, - units=None, - static_target=True, - desc="wing aspect ratio", - targets={ - f"{phase}": ["AR"] for phase in ["br_to_v1", "v1_to_vr", "rto", "rotate", "climb"] - }, - ) - - traj.add_parameter( - "e", - val=801, - opt=False, - units=None, - static_target=True, - desc="Oswald span efficiency factor", - targets={ - f"{phase}": ["e"] for phase in ["br_to_v1", "v1_to_vr", "rto", "rotate", "climb"] - }, - ) - - traj.add_parameter( - "span", - val=35.7, - opt=False, - units="m", - static_target=True, - desc="wingspan", - targets={ - f"{phase}": ["span"] for phase in ["br_to_v1", "v1_to_vr", "rto", "rotate", "climb"] - }, - ) - - traj.add_parameter( - "h_w", - val=1.0, - opt=False, - units="m", - static_target=True, - desc="height of wing above CG", - targets={ - f"{phase}": ["h_w"] for phase in ["br_to_v1", "v1_to_vr", "rto", "rotate", "climb"] - }, - ) - - traj.add_parameter( - "CL0", - val=0.5, - opt=False, - units=None, - static_target=True, - desc="zero-alpha lift coefficient", - targets={ - f"{phase}": ["CL0"] for phase in ["br_to_v1", "v1_to_vr", "rto", "rotate", "climb"] - }, - ) - - traj.add_parameter( - "CL_max", - val=2.0, - opt=False, - units=None, - static_target=True, - desc="maximum lift coefficient for linear fit", - targets={ - f"{phase}": ["CL_max"] - for phase in ["br_to_v1", "v1_to_vr", "rto", "rotate", "climb"] - }, - ) - - traj.add_parameter( - "alpha_max", - val=10.0, - opt=False, - units="deg", - static_target=True, - desc="angle of attack at maximum lift", - targets={ - f"{phase}": ["alpha_max"] - for phase in ["br_to_v1", "v1_to_vr", "rto", "rotate", "climb"] - }, - ) - - # Standard "end of first phase to beginning of second phase" linkages - # Alpha changes from being a parameter in v1_to_vr to a polynomial control - # in rotate, to a dynamic control in `climb`. - traj.link_phases(["br_to_v1", "v1_to_vr"], vars=["time", "r", "v"]) - traj.link_phases(["v1_to_vr", "rotate"], vars=["time", "r", "v", "alpha"]) - traj.link_phases(["rotate", "climb"], vars=["time", "r", "v", "alpha"]) - traj.link_phases(["br_to_v1", "rto"], vars=["time", "r", "v"]) - - # Less common "final value of r must be the match at ends of two phases". - traj.add_linkage_constraint( - phase_a="rto", - var_a="r", - loc_a="final", - phase_b="climb", - var_b="r", - loc_b="final", - ref=1000, - ) - - # Define the constraints and objective for the optimal control problem - v1_to_vr.add_boundary_constraint("v_over_v_stall", loc="final", lower=1.2, ref=100) - - rto.add_boundary_constraint("v", loc="final", equals=0.0, ref=100, linear=True) - - rotate.add_boundary_constraint("F_r", loc="final", equals=0, ref=100000) - - climb.add_boundary_constraint("h", loc="final", equals=35, ref=35, units="ft", linear=True) - climb.add_boundary_constraint("gam", loc="final", equals=5, ref=5, units="deg", linear=True) - climb.add_path_constraint("gam", lower=0, upper=5, ref=5, units="deg") - climb.add_path_constraint("h", lower=0, upper=35, ref=1.0, units="ft") - climb.add_boundary_constraint("v_over_v_stall", loc="final", lower=1.25, ref=1.25) - - rto.add_objective("r", loc="final", ref=1.0) - - for phase_name, phase in traj._phases.items(): - if "T_nominal" in phase.parameter_options: - phase.add_timeseries_output("T_nominal", output_name="T") - if "T_engine_out" in phase.parameter_options: - phase.add_timeseries_output("T_engine_out", output_name="T") - if "T_shutdown" in phase.parameter_options: - phase.add_timeseries_output("T_shutdown", output_name="T") - phase.add_timeseries_output("alpha") - - # - # Setup the problem and set the initial guess - # - p.setup() - - br_to_v1.set_time_val(initial=0.0, duration=35.0) - br_to_v1.set_state_val("r", [0, 2500.0]) - br_to_v1.set_state_val("v", [0.0, 100.0]) - br_to_v1.set_parameter_val("alpha", 0.0, units="deg") - - v1_to_vr.set_time_val(initial=35.0, duration=35.0) - v1_to_vr.set_state_val("r", [2500, 300.0]) - v1_to_vr.set_state_val("v", [100, 110.0]) - v1_to_vr.set_parameter_val("alpha", 0.0, units="deg") - - rto.set_time_val(initial=35.0, duration=35.0) - rto.set_state_val("r", [2500, 5000.0]) - rto.set_state_val("v", [110, 0.0]) - rto.set_parameter_val("alpha", 0.0, units="deg") - - rotate.set_time_val(initial=70.0, duration=5.0) - rotate.set_state_val("r", [1750, 1800.0]) - rotate.set_state_val("v", [80, 85.0]) - rotate.set_control_val("alpha", 0.0, units="deg") - - climb.set_time_val(initial=75.0, duration=15.0) - climb.set_state_val("r", [5000, 5500.0], units="ft") - climb.set_state_val("v", [160, 170.0], units="kn") - climb.set_state_val("h", [0.0, 35.0], units="ft") - climb.set_state_val("gam", [0.0, 5.0], units="deg") - climb.set_control_val("alpha", 5.0, units="deg") - - return p - - @require_pyoptsparse(optimizer="IPOPT") - def test_balanced_field_length_for_docs(self): - for tx in (dm.Radau, dm.GaussLobatto): - p = self._make_problem(tx, optimizer='IPOPT') - - traj = p.model.traj - - result = dm.run_problem(p, run_driver=True, simulate=True) - - sol_db = p.get_outputs_dir() / "dymos_solution.db" - sim_db = traj.sim_prob.get_outputs_dir() / "dymos_simulation.db" - - sol = om.CaseReader(sol_db).get_case("final") - sim = om.CaseReader(sim_db).get_case("final") - - sol_r_f_climb = sol.get_val("traj.climb.timeseries.r")[-1, ...] - sol_r_f_rto = sol.get_val("traj.rto.timeseries.r")[-1, ...] - sim_r_f_climb = sim.get_val("traj.climb.timeseries.r")[-1, ...] - sim_r_f_rto = sim.get_val("traj.rto.timeseries.r")[-1, ...] - - self.assertTrue(result["success"]) - assert_near_equal(2114.387, sol_r_f_climb, tolerance=0.01) - assert_near_equal(2114.387, sol_r_f_rto, tolerance=0.01) - assert_near_equal(2114.387, sim_r_f_climb, tolerance=0.01) - assert_near_equal(2114.387, sim_r_f_rto, tolerance=0.01) - - def test_no_regression(self): - """Test that there are no regressions in the initial values of the driver vars.""" - import json - - # For now we only do this with GaussLobatto until we remove the legacy radau method. - for tx in (dm.GaussLobatto,): - p = self._make_problem(tx, optimizer=None) - dm.run_problem(p, run_driver=False, simulate=False) - driver_vars = p.list_driver_vars(out_stream=None) - - desvars = {} - constraints = {} - objs = {} - - for name, meta in driver_vars["constraints"]: - constraints[name] = meta["val"].tolist() - - for name, meta in driver_vars["design_vars"]: - desvars[name] = meta["val"].tolist() - - for name, meta in driver_vars["objectives"]: - objs[name] = meta["val"].tolist() - - vars = {"constraints": constraints, "design_vars": desvars, "objectives": objs} - - # If we change dymos we might have to regenerate the regression - # data. Dump the JSON and save it in the regression variables - # defined above. - # print(json.dumps(vars, indent=' ')) - - reg_data = regression_data[tx] - - errors = {} - - for var, val in vars["constraints"].items(): - try: - assert_near_equal( - np.asarray(val), reg_data["constraints"][var], tolerance=1.0e-9 - ) - except ValueError as e: - errors[var] = {"actual": val, "expected": reg_data["constraints"][var]} - - for var, val in vars["design_vars"].items(): - try: - assert_near_equal( - np.asarray(val), reg_data["design_vars"][var], tolerance=1.0e-9 - ) - except ValueError as e: - errors[var] = {"actual": val, "expected": reg_data["design_vars"][var]} - - for var, val in vars["objectives"].items(): - try: - assert_near_equal( - np.asarray(val), reg_data["objectives"][var], tolerance=1.0e-9 - ) - except ValueError as e: - errors[var] = {"actual": val, "expected": reg_data["objectives"][var]} - - if errors: - msg = f"Outputs for transcription {str(tx.__name__)} showed a regression.\n" - msg += json.dumps(errors, indent=" ") - self.fail(msg) - - -if __name__ == "__main__": - unittest.main() diff --git a/dymos/examples/brachistochrone/doc/brachistochrone_ode.py b/dymos/examples/brachistochrone/doc/brachistochrone_ode.py deleted file mode 100644 index 151d612005..0000000000 --- a/dymos/examples/brachistochrone/doc/brachistochrone_ode.py +++ /dev/null @@ -1,73 +0,0 @@ -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('g', default=9.80665, desc='gravitational acceleration in m/s**2') - - def setup(self): - nn = self.options['num_nodes'] - - # Inputs - self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s') - - self.add_input('theta', val=np.ones(nn), desc='angle of wire', units='rad') - - 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']) - - self.add_output('check', val=np.zeros(nn), desc='check solution: v/sin(theta) = constant', - units='m/s') - - # Setup partials - ar = np.arange(self.options['num_nodes'], dtype=int) - - self.declare_partials(of='vdot', wrt='theta', rows=ar, cols=ar) - - self.declare_partials(of='xdot', wrt='v', rows=ar, cols=ar) - self.declare_partials(of='xdot', wrt='theta', rows=ar, cols=ar) - - self.declare_partials(of='ydot', wrt='v', rows=ar, cols=ar) - self.declare_partials(of='ydot', wrt='theta', rows=ar, cols=ar) - - self.declare_partials(of='check', wrt='v', rows=ar, cols=ar) - self.declare_partials(of='check', wrt='theta', rows=ar, cols=ar) - - def compute(self, inputs, outputs): - theta = inputs['theta'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - g = self.options['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 = self.options['g'] - v = inputs['v'] - - 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 diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py deleted file mode 100644 index f358a2c07f..0000000000 --- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone.py +++ /dev/null @@ -1,403 +0,0 @@ -import os -import unittest - -try: - import matplotlib - import matplotlib.pyplot as plt - - matplotlib.use('Agg') - plt.style.use('ggplot') -except ImportError: - matplotlib = None - -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -@use_tempdirs -class TestBrachistochroneForDocs(unittest.TestCase): - - def tearDown(self): - for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']: - if os.path.exists(filename): - os.remove(filename) - - def test_brachistochrone_partials(self): - import numpy as np - import openmdao.api as om - from dymos.utils.testing_utils import assert_check_partials - from dymos.examples.brachistochrone.doc.brachistochrone_ode import BrachistochroneODE - - num_nodes = 5 - - p = om.Problem(model=om.Group()) - - ivc = p.model.add_subsystem('vars', om.IndepVarComp()) - ivc.add_output('v', shape=(num_nodes,), units='m/s') - ivc.add_output('theta', shape=(num_nodes,), units='deg') - - p.model.add_subsystem('ode', BrachistochroneODE(num_nodes=num_nodes)) - - p.model.connect('vars.v', 'ode.v') - p.model.connect('vars.theta', 'ode.theta') - - p.setup(force_alloc_complex=True) - - p.set_val('vars.v', 10*np.random.random(num_nodes)) - p.set_val('vars.theta', 10*np.random.uniform(1, 179, num_nodes)) - - p.run_model() - cpd = p.check_partials(method='cs', compact_print=True) - assert_check_partials(cpd) - - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_brachistochrone_for_docs_gauss_lobatto(self): - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - from dymos.examples.plotting import plot_results - from dymos.examples.brachistochrone import BrachistochroneODE - import matplotlib.pyplot as plt - - # - # 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)) - - 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) - - 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 - # - phase.set_time_val(initial=0.0, duration=2.0) - - phase.set_state_val('x', [0, 10]) - phase.set_state_val('y', [10, 5]) - phase.set_state_val('v', [0, 9.9]) - - phase.set_control_val('theta', [5, 100.5]) - - # - # 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) - - # Generate the explicitly simulated trajectory - exp_out = traj.simulate() - - plot_results([('traj.phase0.timeseries.x', 'traj.phase0.timeseries.y', - 'x (m)', 'y (m)'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.theta', - 'time (s)', 'theta (deg)')], - title='Brachistochrone Solution\nHigh-Order Gauss-Lobatto Method', - p_sol=p, p_sim=exp_out) - - plt.show() - - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_brachistochrone_for_docs_radau(self): - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - from dymos.examples.plotting import plot_results - from dymos.examples.brachistochrone import BrachistochroneODE - - # - # 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.Radau(num_segments=10))) - - # - # Set the variables - # - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - - 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) - - 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 - # - phase.set_time_val(initial=0.0, duration=2.0) - - phase.set_state_val('x', [0, 10]) - phase.set_state_val('y', [10, 5]) - phase.set_state_val('v', [0, 9.9]) - - phase.set_control_val('theta', [5, 100.5]) - - # - # 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) - - # Generate the explicitly simulated trajectory - exp_out = traj.simulate() - - plot_results([('traj.phase0.timeseries.x', 'traj.phase0.timeseries.y', - 'x (m)', 'y (m)'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.theta', - 'time (s)', 'theta (deg)')], - title='Brachistochrone Solution\nRadau Pseudospectral Method', - p_sol=p, p_sim=exp_out) - - plt.show() - - @require_pyoptsparse(optimizer='IPOPT') - def test_brachistochrone_for_docs_coloring_demo(self): - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - from dymos.examples.plotting import plot_results - from dymos.examples.brachistochrone import BrachistochroneODE - - # - # Initialize the Problem and the optimization driver - # - p = om.Problem(model=om.Group()) - p.driver = om.pyOptSparseDriver(optimizer='IPOPT') - p.driver.declare_coloring(tol=1.0E-12) - - # - # Create a trajectory and add a phase to it - # - traj = p.model.add_subsystem('traj', dm.Trajectory()) - - # - # In this case the phase has many segments to demonstrate the impact of coloring. - # - phase = traj.add_phase('phase0', - dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.Radau(num_segments=100))) - - # - # Set the variables - # - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - - 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) - - 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 - # - phase.set_time_val(initial=0.0, duration=2.0) - - phase.set_state_val('x', [0, 10]) - phase.set_state_val('y', [10, 5]) - phase.set_state_val('v', [0, 9.9]) - - phase.set_control_val('theta', [5, 100.5]) - - # - # 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) - - # Generate the explicitly simulated trajectory - exp_out = traj.simulate() - - plot_results([('traj.phase0.timeseries.x', 'traj.phase0.timeseries.y', - 'x (m)', 'y (m)'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.theta', - 'time (s)', 'theta (deg)')], - title='Brachistochrone Solution\nRadau Pseudospectral Method', - p_sol=p, p_sim=exp_out) - - plt.show() - - @require_pyoptsparse(optimizer='IPOPT') - def test_brachistochrone_for_docs_coloring_demo_solve_segments(self): - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - from dymos.examples.plotting import plot_results - from dymos.examples.brachistochrone import BrachistochroneODE - - # - # Initialize the Problem and the optimization driver - # - p = om.Problem(model=om.Group()) - p.driver = om.pyOptSparseDriver(optimizer='IPOPT') - p.driver.opt_settings['print_level'] = 0 - # p.driver.declare_coloring() - - # - # Create a trajectory and add a phase to it - # - traj = p.model.add_subsystem('traj', dm.Trajectory()) - - # - # In this case the phase has many segments to demonstrate the impact of coloring. - # - phase = traj.add_phase('phase0', - dm.Phase(ode_class=BrachistochroneODE, - transcription=dm.Radau(num_segments=100, - solve_segments='forward'))) - - # - # Set the variables - # - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - - phase.add_state('x', fix_initial=True) - - phase.add_state('y', fix_initial=True) - - phase.add_state('v', fix_initial=True) - - 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) - - # - # Replace state terminal bounds with nonlinear constraints - # - 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', loc='final', scaler=10) - - p.model.linear_solver = om.DirectSolver() - - # - # Setup the Problem - # - p.setup() - - # - # Set the initial values - # - phase.set_time_val(initial=0.0, duration=2.0) - - phase.set_state_val('x', [0, 10]) - phase.set_state_val('y', [10, 5]) - phase.set_state_val('v', [0, 9.9]) - - phase.set_control_val('theta', [5, 100.5]) - - # - # 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) - - # Generate the explicitly simulated trajectory - exp_out = traj.simulate() - - plot_results([('traj.phase0.timeseries.x', 'traj.phase0.timeseries.y', - 'x (m)', 'y (m)'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.theta', - 'time (s)', 'theta (deg)')], - title='Brachistochrone Solution\nRadau Pseudospectral Method', - p_sol=p, p_sim=exp_out) - - plt.show() - - -if __name__ == '__main__': # pragma: no cover - unittest.main() diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py deleted file mode 100644 index cd13c68966..0000000000 --- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py +++ /dev/null @@ -1,177 +0,0 @@ -import os -import unittest - -try: - import matplotlib - import matplotlib.pyplot as plt - - matplotlib.use('Agg') - plt.style.use('ggplot') -except ImportError: - matplotlib = None - - -from openmdao.utils.testing_utils import use_tempdirs - - -@use_tempdirs -class TestBrachistochroneStaticGravity(unittest.TestCase): - - def tearDown(self): - for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']: - if os.path.exists(filename): - os.remove(filename) - - def test_brachistochrone_partials(self): - import numpy as np - import openmdao.api as om - from dymos.utils.testing_utils import assert_check_partials - from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE - - num_nodes = 5 - - p = om.Problem(model=om.Group()) - - ivc = p.model.add_subsystem('vars', om.IndepVarComp()) - ivc.add_output('v', shape=(num_nodes,), units='m/s') - ivc.add_output('theta', shape=(num_nodes,), units='deg') - - p.model.add_subsystem('ode', BrachistochroneODE(num_nodes=num_nodes, static_gravity=True)) - - p.model.connect('vars.v', 'ode.v') - p.model.connect('vars.theta', 'ode.theta') - - p.setup(force_alloc_complex=True) - - p.set_val('vars.v', 10*np.random.random(num_nodes)) - p.set_val('vars.theta', 10*np.random.uniform(1, 179, num_nodes)) - - p.run_model() - cpd = p.check_partials(method='cs', compact_print=True) - assert_check_partials(cpd) - - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_brachistochrone_static_gravity(self): - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - import matplotlib.pyplot as plt - from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE - - # - # 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, - ode_init_kwargs={'static_gravity': True}, - transcription=dm.GaussLobatto(num_segments=10))) - - # - # Set the variables - # - phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - - phase.add_state('x', rate_source='xdot', - targets=None, - units='m', - fix_initial=True, fix_final=True, solve_segments=False) - - phase.add_state('y', rate_source='ydot', - targets=None, - units='m', - fix_initial=True, fix_final=True, solve_segments=False) - - phase.add_state('v', rate_source='vdot', - targets=['v'], - units='m/s', - fix_initial=True, fix_final=False, solve_segments=False) - - phase.add_control('theta', targets=['theta'], - continuity=True, rate_continuity=True, - units='deg', lower=0.01, upper=179.9) - - phase.add_parameter('g', targets=['g'], static_target=True, opt=False) - - # - # Minimize time at the end of the phase - # - phase.add_objective('time', loc='final', scaler=10) - - # - # Setup the Problem - # - p.setup() - - # - # Set the initial values - # The initial time is fixed, and we set that fixed value here. - # The optimizer is allowed to modify t_duration, but an initial guess is provided here. - # - phase.set_time_val(initial=0, duration=2) - - # Guesses for states are provided at all state_input nodes. - # We use the phase.interpolate method to linearly interpolate values onto the state input nodes. - # Since fix_initial=True for all states and fix_final=True for x and y, the initial or final - # values of the interpolation provided here will not be changed by the optimizer. - phase.set_state_val('x', [0, 10]) - phase.set_state_val('y', [10, 5]) - phase.set_state_val('v', [0, 9.9]) - - # Guesses for controls are provided at all control_input node. - # Here phase.interpolate is used to linearly interpolate values onto the control input nodes. - phase.set_control_val('theta', [5, 100.5]) - - # Set the value for gravitational acceleration. - phase.set_parameter_val('g', 9.80665) - - # - # Solve for the optimal trajectory - # - dm.run_problem(p, simulate=True) - - # Test the results - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, - tolerance=1.0E-3) - - # Load the explicitly simulated trajectory - sim_db = p.model.traj.sim_prob.get_outputs_dir() / 'dymos_simulation.db' - - exp_out = om.CaseReader(sim_db).get_case('final') - - # Extract the timeseries from the implicit solution and the explicit simulation - x = p.get_val('traj.phase0.timeseries.x') - y = p.get_val('traj.phase0.timeseries.y') - t = p.get_val('traj.phase0.timeseries.time') - theta = p.get_val('traj.phase0.timeseries.theta') - - x_exp = exp_out.get_val('traj.phase0.timeseries.x') - y_exp = exp_out.get_val('traj.phase0.timeseries.y') - t_exp = exp_out.get_val('traj.phase0.timeseries.time') - theta_exp = exp_out.get_val('traj.phase0.timeseries.theta') - - fig, axes = plt.subplots(nrows=2, ncols=1) - - axes[0].plot(x, y, 'o') - axes[0].plot(x_exp, y_exp, '-') - axes[0].set_xlabel('x (m)') - axes[0].set_ylabel('y (m)') - - axes[1].plot(t, theta, 'o') - axes[1].plot(t_exp, theta_exp, '-') - axes[1].set_xlabel('time (s)') - axes[1].set_ylabel(r'$\theta$ (deg)') - - plt.show() - - -if __name__ == '__main__': # pragma: no cover - unittest.main() diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_tandem_phases.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_tandem_phases.py deleted file mode 100644 index f774758668..0000000000 --- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_tandem_phases.py +++ /dev/null @@ -1,163 +0,0 @@ -import unittest - -try: - import matplotlib -except ImportError: - matplotlib = None - -import numpy as np -import openmdao.api as om -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -class BrachistochroneArclengthODE(om.ExplicitComponent): - - def initialize(self): - self.options.declare('num_nodes', types=int) - - def setup(self): - nn = self.options['num_nodes'] - - # Inputs - self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s') - self.add_input('theta', val=np.zeros(nn), desc='angle of wire', units='rad') - self.add_output('Sdot', val=np.zeros(nn), desc='rate of change of arclength', units='m/s') - - # Setup partials - arange = np.arange(nn) - - self.declare_partials(of='Sdot', wrt='v', rows=arange, cols=arange) - self.declare_partials(of='Sdot', wrt='theta', rows=arange, cols=arange) - - def compute(self, inputs, outputs): - theta = inputs['theta'] - v = inputs['v'] - outputs['Sdot'] = np.sqrt(1.0 + (1.0/np.tan(theta))**2) * v * np.sin(theta) - - def compute_partials(self, inputs, jacobian): - theta = inputs['theta'] - v = inputs['v'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - tan_theta = np.tan(theta) - cot_theta = 1.0 / tan_theta - csc_theta = 1.0 / sin_theta - - jacobian['Sdot', 'v'] = sin_theta * np.sqrt(1.0 + cot_theta**2) - jacobian['Sdot', 'theta'] = v * (cos_theta * (cot_theta**2 + 1) - cot_theta * csc_theta) / \ - (np.sqrt(1 + cot_theta**2)) - - -@use_tempdirs -class TestBrachistochroneTandemPhases(unittest.TestCase): - - @require_pyoptsparse(optimizer='SLSQP') - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_brachistochrone_tandem_phases(self): - from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE - - import numpy as np - import matplotlib.pyplot as plt - plt.switch_backend('Agg') - import openmdao.api as om - import dymos as dm - - from openmdao.utils.assert_utils import assert_near_equal - - p = om.Problem(model=om.Group()) - - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'SLSQP' - p.driver.declare_coloring() - - # The transcription of the first phase - tx0 = dm.GaussLobatto(num_segments=10, order=3, compressed=False) - - # The transcription for the second phase (and the secondary timeseries outputs from the first phase) - tx1 = dm.Radau(num_segments=20, order=9, compressed=False) - - # - # First Phase: Integrate the standard brachistochrone ODE - # - phase0 = dm.Phase(ode_class=BrachistochroneODE, transcription=tx0) - - p.model.add_subsystem('phase0', phase0) - - phase0.set_time_options(fix_initial=True, duration_bounds=(.5, 10)) - - phase0.add_state('x', fix_initial=True, fix_final=False) - - phase0.add_state('y', fix_initial=True, fix_final=False) - - phase0.add_state('v', fix_initial=True, fix_final=False) - - phase0.add_control('theta', continuity=True, rate_continuity=True, - units='deg', lower=0.01, upper=179.9) - - phase0.add_parameter('g', units='m/s**2', val=9.80665) - - phase0.add_boundary_constraint('x', loc='final', equals=10) - phase0.add_boundary_constraint('y', loc='final', equals=5) - - # Add alternative timeseries output to provide control inputs for the next phase - phase0.add_timeseries('timeseries2', transcription=tx1, subset='control_input') - - # - # Second Phase: Integration of ArcLength - # - phase1 = dm.Phase(ode_class=BrachistochroneArclengthODE, transcription=tx1) - - p.model.add_subsystem('phase1', phase1) - - phase1.set_time_options(fix_initial=True, input_duration=True) - - phase1.add_state('S', fix_initial=True, fix_final=False, - rate_source='Sdot', units='m') - - phase1.add_control('theta', opt=False, units='deg', targets='theta') - phase1.add_control('v', opt=False, units='m/s', targets='v') - - # - # Connect the two phases - # - p.model.connect('phase0.t_duration_val', 'phase1.t_duration') - - p.model.connect('phase0.timeseries2.theta', 'phase1.controls:theta') - p.model.connect('phase0.timeseries2.v', 'phase1.controls:v') - - # Minimize arclength at the end of the second phase - phase1.add_objective('S', loc='final', ref=1) - - p.model.linear_solver = om.DirectSolver() - p.setup(check=True) - - phase0.set_time_val(initial=0.0, duration=2.0) - - phase0.set_state_val('x', [0, 10]) - phase0.set_state_val('y', [10, 5]) - phase0.set_state_val('v', [0, 9.9]) - - phase0.set_control_val('theta', [5, 100.5]) - - phase0.set_parameter_val('g', 9.80665) - - phase1.set_state_val('S', 0.0) - - dm.run_problem(p) - - expected = np.sqrt((10-0)**2 + (10 - 5)**2) - assert_near_equal(p.get_val('phase1.timeseries.S')[-1], expected, tolerance=1.0E-3) - - fig, (ax0, ax1) = plt.subplots(2, 1) - fig.tight_layout() - ax0.plot(p.get_val('phase0.timeseries.x'), p.get_val('phase0.timeseries.y'), '.') - ax0.set_xlabel('x (m)') - ax0.set_ylabel('y (m)') - ax1.plot(p.get_val('phase1.timeseries.time'), p.get_val('phase1.timeseries.S'), '+') - ax1.set_xlabel('t (s)') - ax1.set_ylabel('S (m)') - plt.show() - - -if __name__ == '__main__': # pragma: no cover - unittest.main() diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_upstream_control.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_upstream_control.py deleted file mode 100644 index b56597081b..0000000000 --- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_upstream_control.py +++ /dev/null @@ -1,149 +0,0 @@ -import unittest - -try: - import matplotlib -except ImportError: - matplotlib = None - -from openmdao.utils.testing_utils import use_tempdirs - - -@use_tempdirs -class TestBrachistochroneUpstreamControl(unittest.TestCase): - - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_brachistochrone_upstream_control(self): - import numpy as np - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - - import matplotlib.pyplot as plt - plt.switch_backend('Agg') - - from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE - - # - # Define the OpenMDAO problem - # - p = om.Problem(model=om.Group()) - - # Instantiate the transcription so we can get the number of nodes from it while - # building the problem. - tx = dm.GaussLobatto(num_segments=10, order=3) - - # Add an indep var comp to provide the external control values - ivc = p.model.add_subsystem('control_ivc', om.IndepVarComp(), promotes_outputs=['*']) - - # Add the output to provide the values of theta at the control input nodes of the transcription. - ivc.add_output('theta', shape=(tx.grid_data.subset_num_nodes['control_input']), units='rad') - - # Add this external control as a design variable - p.model.add_design_var('theta', units='rad', lower=1.0E-5, upper=np.pi) - # Connect this to controls:theta in the appropriate phase. - # connect calls are cached, so we can do this before we actually add the trajectory to the problem. - p.model.connect('theta', 'traj.phase0.controls:theta') - - # - # Define a Trajectory object - # - traj = dm.Trajectory() - - p.model.add_subsystem('traj', subsys=traj) - - # - # Define a Dymos Phase object with GaussLobatto Transcription - # - phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=tx) - - traj.add_phase(name='phase0', phase=phase) - - # - # Set the time options - # Time has no targets in our ODE. - # We fix the initial time so that the it is not a design variable in the optimization. - # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. - # - phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') - - # - # Set the time options - # Initial values of positions and velocity are all fixed. - # The final value of position are fixed, but the final velocity is a free variable. - # The equations of motion are not functions of position, so 'x' and 'y' have no targets. - # The rate source points to the output in the ODE which provides the time derivative of the - # given state. - phase.add_state('x', fix_initial=True, fix_final=True, units='m', rate_source='xdot') - phase.add_state('y', fix_initial=True, fix_final=True, units='m', rate_source='ydot') - phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', - rate_source='vdot', targets=['v']) - - # Define theta as a control. - # Use opt=False to allow it to be connected to an external source. - # Arguments lower and upper are no longer valid for an input control. - phase.add_control(name='theta', targets=['theta'], opt=False) - - # Minimize final time. - phase.add_objective('time', loc='final') - - # Set the driver. - p.driver = om.ScipyOptimizeDriver() - - # Allow OpenMDAO to automatically determine our sparsity pattern. - # Doing so can significant speed up the execution of Dymos. - p.driver.declare_coloring() - - # Setup the problem - p.setup(check=True) - - # Now that the OpenMDAO problem is setup, we can set the values of the states and controls. - p.set_val('traj.phase0.states:x', phase.interp('x', [0, 10]), units='m') - - p.set_val('traj.phase0.states:y', phase.interp('y', [10, 5]), units='m') - - p.set_val('traj.phase0.states:v', phase.interp('v', [0, 5]), units='m/s') - - p.set_val('traj.phase0.controls:theta', phase.interp('theta', [90, 90]), units='deg') - - # Run the driver to solve the problem - p.run_driver() - - # Test the results - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, - tolerance=1.0E-3) - - # Check the validity of our results by using scipy.integrate.solve_ivp to - # integrate the solution. - sim_out = traj.simulate() - - # Plot the results - fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 4.5)) - - axes[0].plot(p.get_val('traj.phase0.timeseries.x'), - p.get_val('traj.phase0.timeseries.y'), - 'ro', label='solution') - - axes[0].plot(sim_out.get_val('traj.phase0.timeseries.x'), - sim_out.get_val('traj.phase0.timeseries.y'), - 'b-', label='simulation') - - axes[0].set_xlabel('x (m)') - axes[0].set_ylabel('y (m/s)') - axes[0].legend() - axes[0].grid() - - axes[1].plot(p.get_val('traj.phase0.timeseries.time'), - p.get_val('traj.phase0.timeseries.theta', units='deg'), - 'ro', label='solution') - - axes[1].plot(sim_out.get_val('traj.phase0.timeseries.time'), - sim_out.get_val('traj.phase0.timeseries.theta', units='deg'), - 'b-', label='simulation') - - axes[1].set_xlabel('time (s)') - axes[1].set_ylabel(r'$\theta$ (deg)') - axes[1].legend() - axes[1].grid() - - plt.show() diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_upstream_state.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_upstream_state.py deleted file mode 100644 index 86611b0fde..0000000000 --- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_upstream_state.py +++ /dev/null @@ -1,167 +0,0 @@ -import unittest - -try: - import matplotlib -except ImportError: - matplotlib = None - -from openmdao.utils.testing_utils import use_tempdirs - - -@use_tempdirs -class TestBrachistochroneUpstreamState(unittest.TestCase): - - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_brachistochrone_upstream_state(self): - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - - import matplotlib.pyplot as plt - plt.switch_backend('Agg') - - from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE - - # - # Define the OpenMDAO problem - # - p = om.Problem(model=om.Group()) - - # Instantiate the transcription so we can get the number of nodes from it while - # building the problem. - tx = dm.GaussLobatto(num_segments=10, order=3) - - # Add an indep var comp to provide the external control values - ivc = p.model.add_subsystem('states_ivc', om.IndepVarComp(), promotes_outputs=['*']) - - # Add the output to provide the values of theta at the control input nodes of the transcription. - ivc.add_output('x0', shape=(1,), units='m') - - # Connect x0 to the state error component so we can constrain the given value of x0 - # to be equal to the value chosen in the phase. - p.model.connect('x0', 'state_error_comp.x0_target') - p.model.connect('traj.phase0.timeseries.x', 'state_error_comp.x0_actual', src_indices=[0], flat_src_indices=True) - - # - # Define a Trajectory object - # - traj = dm.Trajectory() - - p.model.add_subsystem('traj', subsys=traj) - - p.model.add_subsystem('state_error_comp', - om.ExecComp('x0_error = x0_target - x0_actual', - x0_error={'units': 'm'}, - x0_target={'units': 'm'}, - x0_actual={'units': 'm'})) - - p.model.add_constraint('state_error_comp.x0_error', equals=0.0) - - # - # Define a Dymos Phase object with GaussLobatto Transcription - # - phase = dm.Phase(ode_class=BrachistochroneODE, - transcription=tx) - - traj.add_phase(name='phase0', phase=phase) - - # - # Set the time options - # Time has no targets in our ODE. - # We fix the initial time so that the it is not a design variable in the optimization. - # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. - # - phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') - - # - # Set the time options - # Initial values of positions and velocity are all fixed. - # The final value of position are fixed, but the final velocity is a free variable. - # The equations of motion are not functions of position, so 'x' and 'y' have no targets. - # The rate source points to the output in the ODE which provides the time derivative of the - # given state. - phase.add_state('x', fix_initial=False, fix_final=True, units='m', rate_source='xdot') - phase.add_state('y', fix_initial=True, fix_final=True, units='m', rate_source='ydot') - phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', - rate_source='vdot', targets=['v']) - - # Define theta as a control. - # Use opt=False to allow it to be connected to an external source. - # Arguments lower and upper are no longer valid for an input control. - phase.add_control(name='theta', units='rad', targets=['theta']) - - # Minimize final time. - phase.add_objective('time', loc='final') - - # Set the driver. - p.driver = om.ScipyOptimizeDriver() - - # Allow OpenMDAO to automatically determine our sparsity pattern. - # Doing so can significant speed up the execution of Dymos. - p.driver.declare_coloring() - - # Setup the problem - p.setup(check=True) - - # Now that the OpenMDAO problem is setup, we can set the values of the states. - p.set_val('x0', 0.0, units='m') - - # Here we're intentially setting the intiial x value to something other than zero, just - # to demonstrate that the optimizer brings it back in line with the value of x0 set above. - p.set_val('traj.phase0.states:x', - phase.interp('x', [1, 10]), - units='m') - - p.set_val('traj.phase0.states:y', - phase.interp('y', [10, 5]), - units='m') - - p.set_val('traj.phase0.states:v', - phase.interp('v', [0, 5]), - units='m/s') - - p.set_val('traj.phase0.controls:theta', - phase.interp('theta', [90, 90]), - units='deg') - - # Run the driver to solve the problem - dm.run_problem(p, make_plots=True) - - # Test the results - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, - tolerance=1.0E-3) - - # Check the validity of our results by using scipy.integrate.solve_ivp to - # integrate the solution. - sim_out = traj.simulate() - - # Plot the results - fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 4.5)) - - axes[0].plot(p.get_val('traj.phase0.timeseries.x'), - p.get_val('traj.phase0.timeseries.y'), - 'ro', label='solution') - - axes[0].plot(sim_out.get_val('traj.phase0.timeseries.x'), - sim_out.get_val('traj.phase0.timeseries.y'), - 'b-', label='simulation') - - axes[0].set_xlabel('x (m)') - axes[0].set_ylabel('y (m/s)') - axes[0].legend() - axes[0].grid() - - axes[1].plot(p.get_val('traj.phase0.timeseries.time'), - p.get_val('traj.phase0.timeseries.theta', units='deg'), - 'ro', label='solution') - - axes[1].plot(sim_out.get_val('traj.phase0.timeseries.time'), - sim_out.get_val('traj.phase0.timeseries.theta', units='deg'), - 'b-', label='simulation') - - axes[1].set_xlabel('time (s)') - axes[1].set_ylabel(r'$\theta$ (deg)') - axes[1].legend() - axes[1].grid() - - plt.show() diff --git a/dymos/examples/cannonball/doc/__init__.py b/dymos/examples/cannonball/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py deleted file mode 100644 index 81ebd9008c..0000000000 --- a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py +++ /dev/null @@ -1,256 +0,0 @@ -import unittest - -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -@use_tempdirs -class TestTwoPhaseCannonballForDocs(unittest.TestCase): - - @require_pyoptsparse(optimizer='SLSQP') - def test_two_phase_cannonball_for_docs(self): - import numpy as np - - import openmdao.api as om - from openmdao.components.interp_util.interp import InterpND - from openmdao.utils.assert_utils import assert_near_equal - - import dymos as dm - from dymos.models.atmosphere.atmos_1976 import USatm1976Data - - ############################################# - # Component for the design part of the model - ############################################# - class CannonballSizeComp(om.ExplicitComponent): - """ - Compute the area and mass of a cannonball with a given radius and density. - - Notes - ----- - This component is not vectorized with 'num_nodes' as is the usual way - with Dymos, but is instead intended to compute a scalar mass and reference - area from scalar radius and density inputs. This component does not reside - in the ODE but instead its outputs are connected to the trajectory via - input design parameters. - """ - def setup(self): - self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m') - self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3') - - self.add_output(name='mass', shape=(1,), desc='cannonball mass', units='kg') - self.add_output(name='S', shape=(1,), desc='aerodynamic reference area', units='m**2') - - self.declare_partials(of='mass', wrt='dens') - self.declare_partials(of='mass', wrt='radius') - - self.declare_partials(of='S', wrt='radius') - - def compute(self, inputs, outputs): - radius = inputs['radius'] - dens = inputs['dens'] - - outputs['mass'] = (4/3.) * dens * np.pi * radius ** 3 - outputs['S'] = np.pi * radius ** 2 - - def compute_partials(self, inputs, partials): - radius = inputs['radius'] - dens = inputs['dens'] - - partials['mass', 'dens'] = (4/3.) * np.pi * radius ** 3 - partials['mass', 'radius'] = 4. * dens * np.pi * radius ** 2 - - partials['S', 'radius'] = 2 * np.pi * radius - - ############################################# - # Build the ODE class - ############################################# - class CannonballODE(om.ExplicitComponent): - """ - Cannonball ODE assuming flat earth and accounting for air resistance - """ - - def initialize(self): - self.options.declare('num_nodes', types=int) - - def setup(self): - nn = self.options['num_nodes'] - - # static parameters - self.add_input('m', units='kg') - self.add_input('S', units='m**2') - # 0.5 good assumption for a sphere - self.add_input('CD', 0.5) - - # time varying inputs - self.add_input('h', units='m', shape=nn) - self.add_input('v', units='m/s', shape=nn) - self.add_input('gam', units='rad', shape=nn) - - # state rates - self.add_output('v_dot', shape=nn, units='m/s**2', tags=['dymos.state_rate_source:v']) - self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam']) - self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h']) - self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r']) - self.add_output('ke', shape=nn, units='J') - - # Ask OpenMDAO to compute the partial derivatives using finite-difference - # with a partial coloring algorithm for improved performance, and use - # a graph coloring algorithm to automatically detect the sparsity pattern. - self.declare_coloring(wrt='*', method='fd') - - 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') - - def compute(self, inputs, outputs): - - gam = inputs['gam'] - v = inputs['v'] - h = inputs['h'] - m = inputs['m'] - S = inputs['S'] - CD = inputs['CD'] - - GRAVITY = 9.80665 # m/s**2 - - rho = self.rho_interp.interpolate(inputs['h']) - - q = 0.5*rho*inputs['v']**2 - qS = q * S - D = qS * CD - cgam = np.cos(gam) - sgam = np.sin(gam) - outputs['v_dot'] = - D/m-GRAVITY*sgam - outputs['gam_dot'] = -(GRAVITY/v)*cgam - outputs['h_dot'] = v*sgam - outputs['r_dot'] = v*cgam - outputs['ke'] = 0.5*m*v**2 - - ############################################# - # Setup the Dymos problem - ############################################# - - p = om.Problem(model=om.Group()) - - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'SLSQP' - p.driver.declare_coloring() - - p.model.add_subsystem('size_comp', CannonballSizeComp(), - promotes_inputs=['radius', 'dens']) - p.model.set_input_defaults('dens', val=7.87, units='g/cm**3') - p.model.add_design_var('radius', lower=0.01, upper=0.10, - ref0=0.01, ref=0.10, units='m') - - traj = p.model.add_subsystem('traj', dm.Trajectory()) - - transcription = dm.Radau(num_segments=5, order=3, compressed=True) - ascent = dm.Phase(ode_class=CannonballODE, transcription=transcription) - - ascent = traj.add_phase('ascent', ascent) - - # All initial states except flight path angle are fixed - # Final flight path angle is fixed (we will set it to zero - # so that the phase ends at apogee). - # The output of the ODE which provides the rate source for each state - # is obtained from the tags used on those outputs in the ODE. - # The units of the states are automatically inferred by multiplying the units - # of those rates by the time units. - ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100), - duration_ref=100, units='s') - ascent.set_state_options('r', fix_initial=True, fix_final=False) - ascent.set_state_options('h', fix_initial=True, fix_final=False) - ascent.set_state_options('gam', fix_initial=False, fix_final=True) - ascent.set_state_options('v', fix_initial=False, fix_final=False) - - ascent.add_parameter('S', units='m**2', static_target=True) - ascent.add_parameter('m', units='kg', static_target=True) - - # Limit the muzzle energy - ascent.add_boundary_constraint('ke', loc='initial', - upper=400000, lower=0, ref=100000) - - # Second Phase (descent) - transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True) - descent = dm.Phase(ode_class=CannonballODE, transcription=transcription) - - traj.add_phase('descent', descent) - - # All initial states and time are free, since - # they will be linked to the final states of ascent. - # Final altitude is fixed, because we will set - # it to zero so that the phase ends at ground impact) - descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100), - duration_ref=100, units='s') - descent.add_state('r') - descent.add_state('h', fix_initial=False, fix_final=True) - descent.add_state('gam', fix_initial=False, fix_final=False) - descent.add_state('v', fix_initial=False, fix_final=False) - - descent.add_parameter('S', units='m**2', static_target=True) - descent.add_parameter('m', units='kg', static_target=True) - - descent.add_objective('r', loc='final', scaler=-1.0) - - # Add internally-managed design parameters to the trajectory. - traj.add_parameter('CD', - targets={'ascent': ['CD'], 'descent': ['CD']}, - val=0.5, units=None, opt=False, static_target=True) - - # Add externally-provided design parameters to the trajectory. - # In this case, we connect 'm' to pre-existing input parameters - # named 'm' in each phase. - traj.add_parameter('m', units='kg', val=1.0, - targets={'ascent': 'm', 'descent': 'm'}, static_target=True) - - # In this case, by omitting targets, we're connecting these - # parameters to parameters with the same name in each phase. - traj.add_parameter('S', units='m**2', val=0.005, static_target=True) - - # Link Phases (link time and all state variables) - traj.link_phases(phases=['ascent', 'descent'], vars=['*']) - - # Issue Connections - p.model.connect('size_comp.mass', 'traj.parameters:m') - p.model.connect('size_comp.S', 'traj.parameters:S') - - # A linear solver at the top level can improve performance. - p.model.linear_solver = om.DirectSolver() - - # Finish Problem Setup - p.setup() - - ############################################# - # Set constants and initial guesses - ############################################# - p.set_val('radius', 0.05, units='m') - p.set_val('dens', 7.87, units='g/cm**3') - - traj.set_parameter_val('CD', 0.5) - - ascent.set_time_val(initial=0.0, duration=10.0) - - ascent.set_state_val('r', [0, 100]) - ascent.set_state_val('h', [0, 100]) - ascent.set_state_val('v', [200, 150]) - ascent.set_state_val('gam', [25, 0], units='deg') - - descent.set_time_val(initial=10.0, duration=10.0) - - descent.set_state_val('r', [100, 200]) - descent.set_state_val('h', [100, 0]) - descent.set_state_val('v', [150, 200]) - descent.set_state_val('gam', [0, -45], units='deg') - - ##################################################### - # Run the optimization and final explicit simulation - ##################################################### - dm.run_problem(p, simulate=True) - - assert_near_equal(p.get_val('traj.descent.states:r')[-1], - 3183.25, tolerance=1.0E-2) - - -if __name__ == '__main__': # pragma: no cover - unittest.main() diff --git a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py b/dymos/examples/cannonball/test/test_cannonball_implicit_duration.py similarity index 100% rename from dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py rename to dymos/examples/cannonball/test/test_cannonball_implicit_duration.py diff --git a/dymos/examples/double_integrator/doc/__init__.py b/dymos/examples/double_integrator/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/double_integrator/doc/test_doc_double_integrator.py b/dymos/examples/double_integrator/doc/test_doc_double_integrator.py deleted file mode 100644 index 8e3dbbe264..0000000000 --- a/dymos/examples/double_integrator/doc/test_doc_double_integrator.py +++ /dev/null @@ -1,111 +0,0 @@ -import os -import unittest - -try: - import matplotlib - import matplotlib.pyplot as plt - - matplotlib.use('Agg') - plt.style.use('ggplot') -except ImportError: - matplotlib = None - -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -@use_tempdirs -class TestDoubleIntegratorForDocs(unittest.TestCase): - - @classmethod - def tearDownClass(cls): - for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out']: - if os.path.exists(filename): - os.remove(filename) - - @require_pyoptsparse(optimizer='SLSQP') - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_double_integrator_for_docs(self): - import matplotlib.pyplot as plt - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - from dymos.examples.plotting import plot_results - from dymos.examples.double_integrator.double_integrator_ode import DoubleIntegratorODE - - # Initialize the problem and assign the driver - p = om.Problem(model=om.Group()) - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'SLSQP' - p.driver.declare_coloring() - - # Setup the trajectory and its phase - traj = p.model.add_subsystem('traj', dm.Trajectory()) - - transcription = dm.Radau(num_segments=30, order=3, compressed=False) - - phase = traj.add_phase('phase0', - dm.Phase(ode_class=DoubleIntegratorODE, transcription=transcription)) - - # - # Set the options for our variables. - # - phase.set_time_options(fix_initial=True, fix_duration=True, units='s') - phase.add_state('v', fix_initial=True, fix_final=True, rate_source='u', units='m/s') - phase.add_state('x', fix_initial=True, rate_source='v', units='m') - - phase.add_control('u', units='m/s**2', scaler=0.01, continuity=False, rate_continuity=False, - rate2_continuity=False, shape=(1, ), lower=-1.0, upper=1.0) - - # - # Maximize distance travelled. - # - phase.add_objective('x', loc='final', scaler=-1) - - p.model.linear_solver = om.DirectSolver() - - # - # Setup the problem and set our initial values. - # - p.setup(check=True) - - phase.set_time_val(initial=0.0, duration=1.0) - phase.set_state_val('x', [0, 0.25]) - phase.set_state_val('v', [0, 0]) - phase.set_control_val('u', [1, -1]) - - # - # Solve the problem. - # - dm.run_problem(p) - - # - # Verify that the results are correct. - # - x = p.get_val('traj.phase0.timeseries.x') - v = p.get_val('traj.phase0.timeseries.v') - - assert_near_equal(x[0], 0.0, tolerance=1.0E-4) - assert_near_equal(x[-1], 0.25, tolerance=1.0E-4) - - assert_near_equal(v[0], 0.0, tolerance=1.0E-4) - assert_near_equal(v[-1], 0.0, tolerance=1.0E-4) - - # - # Simulate the explicit solution and plot the results. - # - exp_out = traj.simulate() - - plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.x', - 'time (s)', 'x $(m)$'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.v', - 'time (s)', 'v $(m/s)$'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.u', - 'time (s)', 'u $(m/s^2)$')], - title='Double Integrator Solution\nRadau Pseudospectral Method', - p_sol=p, p_sim=exp_out) - - plt.show() - - -if __name__ == "__main__": - unittest.main() diff --git a/dymos/examples/finite_burn_orbit_raise/doc/__init__.py b/dymos/examples/finite_burn_orbit_raise/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/finite_burn_orbit_raise/doc/test_doc_finite_burn_orbit_raise.py b/dymos/examples/finite_burn_orbit_raise/doc/test_doc_finite_burn_orbit_raise.py deleted file mode 100644 index b5e68a5620..0000000000 --- a/dymos/examples/finite_burn_orbit_raise/doc/test_doc_finite_burn_orbit_raise.py +++ /dev/null @@ -1,241 +0,0 @@ -import unittest - -try: - import matplotlib.pyplot as plt - - plt.switch_backend('Agg') - plt.style.use('ggplot') -except ImportError: - plt = None - -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -@use_tempdirs -class TestFiniteBurnOrbitRaise(unittest.TestCase): - - @require_pyoptsparse(optimizer='IPOPT') - @unittest.skipIf(plt is None, "This test requires matplotlib") - def test_finite_burn_orbit_raise(self): - import numpy as np - import matplotlib.pyplot as plt - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - - import dymos as dm - from dymos.examples.finite_burn_orbit_raise.finite_burn_eom import FiniteBurnODE - - p = om.Problem(model=om.Group()) - - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'IPOPT' - p.driver.declare_coloring() - - traj = dm.Trajectory() - - traj.add_parameter('c', opt=False, val=1.5, units='DU/TU', - targets={'burn1': ['c'], 'coast': ['c'], 'burn2': ['c']}) - - # First Phase (burn) - - burn1 = dm.Phase(ode_class=FiniteBurnODE, - transcription=dm.Radau(num_segments=5, order=3, compressed=False)) - - burn1 = traj.add_phase('burn1', burn1) - - 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=100.0, - rate_source='r_dot', units='DU') - burn1.add_state('theta', fix_initial=True, fix_final=False, defect_scaler=100.0, - rate_source='theta_dot', units='rad') - burn1.add_state('vr', fix_initial=True, fix_final=False, defect_scaler=100.0, - rate_source='vr_dot', units='DU/TU') - burn1.add_state('vt', fix_initial=True, fix_final=False, defect_scaler=100.0, - rate_source='vt_dot', units='DU/TU') - burn1.add_state('accel', fix_initial=True, fix_final=False, - 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, rate_continuity_scaler=0.001, rate2_continuity_scaler=0.001, - lower=-30, upper=30) - # Second Phase (Coast) - coast = dm.Phase(ode_class=FiniteBurnODE, - transcription=dm.Radau(num_segments=5, order=3, compressed=False)) - - 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=100.0, - rate_source='r_dot', targets=['r'], units='DU') - coast.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=100.0, - rate_source='theta_dot', targets=['theta'], units='rad') - coast.add_state('vr', fix_initial=False, fix_final=False, defect_scaler=100.0, - rate_source='vr_dot', targets=['vr'], units='DU/TU') - coast.add_state('vt', fix_initial=False, fix_final=False, defect_scaler=100.0, - rate_source='vt_dot', targets=['vt'], units='DU/TU') - coast.add_state('accel', fix_initial=True, fix_final=True, - rate_source='at_dot', targets=['accel'], units='DU/TU**2') - coast.add_state('deltav', fix_initial=False, fix_final=False, - rate_source='deltav_dot', units='DU/TU') - - coast.add_parameter('u1', opt=False, val=0.0, units='deg', targets=['u1']) - - # Third Phase (burn) - burn2 = dm.Phase(ode_class=FiniteBurnODE, - transcription=dm.Radau(num_segments=5, order=3, compressed=False)) - - traj.add_phase('coast', coast) - traj.add_phase('burn2', burn2) - - 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=100.0, - rate_source='r_dot', units='DU') - burn2.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=100.0, - rate_source='theta_dot', units='rad') - burn2.add_state('vr', fix_initial=False, fix_final=True, defect_scaler=1000.0, - rate_source='vr_dot', units='DU/TU') - burn2.add_state('vt', fix_initial=False, fix_final=True, defect_scaler=1000.0, - rate_source='vt_dot', units='DU/TU') - burn2.add_state('accel', fix_initial=False, fix_final=False, defect_scaler=1.0, - rate_source='at_dot', units='DU/TU**2') - burn2.add_state('deltav', fix_initial=False, fix_final=False, defect_scaler=1.0, - rate_source='deltav_dot', units='DU/TU') - - burn2.add_objective('deltav', loc='final', scaler=100.0) - - burn2.add_control('u1', rate_continuity=True, rate2_continuity=True, units='deg', - scaler=0.01, lower=-90, upper=90) - - burn1.add_timeseries_output('pos_x') - coast.add_timeseries_output('pos_x') - burn2.add_timeseries_output('pos_x') - - burn1.add_timeseries_output('pos_y') - coast.add_timeseries_output('pos_y') - burn2.add_timeseries_output('pos_y') - - # Link Phases - traj.link_phases(phases=['burn1', 'coast', 'burn2'], - vars=['time', 'r', 'theta', 'vr', 'vt', 'deltav']) - - traj.link_phases(phases=['burn1', 'burn2'], vars=['accel']) - - p.model.add_subsystem('traj', subsys=traj) - - # Finish Problem Setup - - # Needed to move the direct solver down into the phases for use with MPI. - # - After moving down, used fewer iterations (about 30 less) - - p.driver.add_recorder(om.SqliteRecorder('two_burn_orbit_raise_example.db')) - - p.setup(check=True, mode='fwd') - - # 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.set_time_val(initial=0.0, duration=2.25) - burn1.set_state_val('r', [1, 1.5]) - burn1.set_state_val('theta', [0, 1.7]) - burn1.set_state_val('vr', [0, 0]) - burn1.set_state_val('vt', [1, 1]) - burn1.set_state_val('accel', [0.1, 0.0]) - burn1.set_state_val('deltav', [0, 0.1]) - burn1.set_control_val('u1', [-3.5, 13.0]) - - coast.set_time_val(initial=2.25, duration=3.0) - coast.set_state_val('r', [1.3, 1.5]) - coast.set_state_val('theta', [2.1767, 1.7]) - coast.set_state_val('vr', [0.3285, 0]) - coast.set_state_val('vt', [0.97, 1]) - coast.set_state_val('accel', [0, 0]) - - burn2.set_time_val(initial=5.25, duration=1.75) - burn2.set_state_val('r', [1, 3]) - burn2.set_state_val('theta', [0, 4]) - burn2.set_state_val('vr', [0, 0]) - burn2.set_state_val('vt', [1, np.sqrt(1 / 3)]) - burn2.set_state_val('accel', [0.1, 0.0]) - burn2.set_state_val('deltav', [0.1, 0.2]) - burn2.set_control_val('u1', [0, 0]) - - dm.run_problem(p) - - assert_near_equal(p.get_val('traj.burn2.states:deltav')[-1], 0.3995, - tolerance=2.0E-3) - - # - # Plot results - # - traj = p.model.traj - exp_out = traj.simulate() - - fig = plt.figure(figsize=(8, 4)) - fig.suptitle('Two Burn Orbit Raise Solution') - ax_u1 = plt.subplot2grid((2, 2), (0, 0)) - ax_deltav = plt.subplot2grid((2, 2), (1, 0)) - ax_xy = plt.subplot2grid((2, 2), (0, 1), rowspan=2) - - span = np.linspace(0, 2 * np.pi, 100) - ax_xy.plot(np.cos(span), np.sin(span), 'k--', lw=1) - ax_xy.plot(3 * np.cos(span), 3 * np.sin(span), 'k--', lw=1) - ax_xy.set_xlim(-4.5, 4.5) - ax_xy.set_ylim(-4.5, 4.5) - - ax_xy.set_xlabel(r'x ($R_e$)') - ax_xy.set_ylabel(r'y ($R_e$)') - - ax_u1.set_xlabel(r'time ($TU$)') - ax_u1.set_ylabel(r'$u_1$ ($deg$)') - ax_u1.grid(True) - - ax_deltav.set_xlabel('time ($TU$)') - ax_deltav.set_ylabel(r'${\Delta}v$ ($DU/TU$)') - ax_deltav.grid(True) - - t_sol = dict((phs, p.get_val('traj.{0}.timeseries.time'.format(phs))) - for phs in ['burn1', 'coast', 'burn2']) - x_sol = dict((phs, p.get_val('traj.{0}.timeseries.pos_x'.format(phs))) - for phs in ['burn1', 'coast', 'burn2']) - y_sol = dict((phs, p.get_val('traj.{0}.timeseries.pos_y'.format(phs))) - for phs in ['burn1', 'coast', 'burn2']) - dv_sol = dict((phs, p.get_val('traj.{0}.timeseries.deltav'.format(phs))) - for phs in ['burn1', 'coast', 'burn2']) - u1_sol = dict((phs, p.get_val('traj.{0}.timeseries.u1'.format(phs), units='deg')) - for phs in ['burn1', 'burn2']) - - t_exp = dict((phs, exp_out.get_val('traj.{0}.timeseries.time'.format(phs))) - for phs in ['burn1', 'coast', 'burn2']) - x_exp = dict((phs, exp_out.get_val('traj.{0}.timeseries.pos_x'.format(phs))) - for phs in ['burn1', 'coast', 'burn2']) - y_exp = dict((phs, exp_out.get_val('traj.{0}.timeseries.pos_y'.format(phs))) - for phs in ['burn1', 'coast', 'burn2']) - dv_exp = dict((phs, exp_out.get_val('traj.{0}.timeseries.deltav'.format(phs))) - for phs in ['burn1', 'coast', 'burn2']) - u1_exp = dict((phs, exp_out.get_val('traj.{0}.timeseries.u1'.format(phs), - units='deg')) - for phs in ['burn1', 'burn2']) - - for phs in ['burn1', 'coast', 'burn2']: - try: - ax_u1.plot(t_exp[phs], u1_exp[phs], '-', marker=None, color='C0') - ax_u1.plot(t_sol[phs], u1_sol[phs], 'o', mfc='C1', mec='C1', ms=3) - except KeyError: - pass - - ax_deltav.plot(t_exp[phs], dv_exp[phs], '-', marker=None, color='C0') - ax_deltav.plot(t_sol[phs], dv_sol[phs], 'o', mfc='C1', mec='C1', ms=3) - - ax_xy.plot(x_exp[phs], y_exp[phs], '-', marker=None, color='C0', label='explicit') - ax_xy.plot(x_sol[phs], y_sol[phs], 'o', mfc='C1', mec='C1', ms=3, label='implicit') - - plt.show() - - -if __name__ == '__main__': # pragma: no cover - unittest.main() diff --git a/dymos/examples/hyper_sensitive/doc/__init__.py b/dymos/examples/hyper_sensitive/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/hyper_sensitive/doc/test_doc_hyper_sensitive.py b/dymos/examples/hyper_sensitive/doc/test_doc_hyper_sensitive.py deleted file mode 100644 index 71e5731f75..0000000000 --- a/dymos/examples/hyper_sensitive/doc/test_doc_hyper_sensitive.py +++ /dev/null @@ -1,121 +0,0 @@ -from __future__ import print_function, division, absolute_import - -import os -import unittest -import numpy as np - -try: - import matplotlib - import matplotlib.pyplot as plt - - matplotlib.use('Agg') - plt.style.use('ggplot') -except ImportError: - matplotlib = None - -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -tf = 10.0 - - -def solution(): - sqrt_two = np.sqrt(2) - val = sqrt_two * tf - c1 = (1.5 * np.exp(-val) - 1) / (np.exp(-val) - np.exp(val)) - c2 = (1 - 1.5 * np.exp(val)) / (np.exp(-val) - np.exp(val)) - - ui = c1 * (1 + sqrt_two) + c2 * (1 - sqrt_two) - uf = c1 * (1 + sqrt_two) * np.exp(val) + c2 * (1 - sqrt_two) * np.exp(-val) - J = 0.5 * (c1 ** 2 * (1 + sqrt_two) * np.exp(2 * val) + c2 ** 2 * (1 - sqrt_two) * np.exp(-2 * val) - - (1 + sqrt_two) * c1 ** 2 - (1 - sqrt_two) * c2 ** 2) - return ui, uf, J - - -@use_tempdirs -class TestHyperSensitive(unittest.TestCase): - - @classmethod - def tearDownClass(cls): - for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out']: - if os.path.exists(filename): - os.remove(filename) - - @require_pyoptsparse(optimizer='SLSQP') - def test_hyper_sensitive_for_docs(self): - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - from dymos.examples.plotting import plot_results - from dymos.examples.hyper_sensitive.hyper_sensitive_ode import HyperSensitiveODE - - # Initialize the problem and assign the driver - p = om.Problem(model=om.Group()) - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'SLSQP' - p.driver.declare_coloring() - - # Setup the trajectory and its phase - traj = p.model.add_subsystem('traj', dm.Trajectory()) - - transcription = dm.Radau(num_segments=30, order=3, compressed=False) - - phase = traj.add_phase('phase0', - dm.Phase(ode_class=HyperSensitiveODE, transcription=transcription)) - - phase.set_time_options(fix_initial=True, fix_duration=True) - phase.add_state('x', fix_initial=True, fix_final=False, rate_source='x_dot', targets=['x']) - phase.add_state('xL', fix_initial=True, fix_final=False, rate_source='L', targets=['xL']) - phase.add_control('u', opt=True, targets=['u']) - - phase.add_boundary_constraint('x', loc='final', equals=1) - - phase.add_objective('xL', loc='final') - - p.setup(check=True) - - phase.set_time_val(initial=0, duration=tf) - phase.set_state_val('x', [1.5, 1]) - phase.set_state_val('xL', [0, 1]) - phase.set_control_val('u', [-0.6, 2.4]) - - # - # Solve the problem. - # - dm.run_problem(p) - - # - # Verify that the results are correct. - # - ui, uf, J = solution() - - assert_near_equal(p.get_val('traj.phase0.timeseries.u')[0], - ui, - tolerance=1.5e-2) - - assert_near_equal(p.get_val('traj.phase0.timeseries.u')[-1], - uf, - tolerance=1.5e-2) - - assert_near_equal(p.get_val('traj.phase0.timeseries.xL')[-1], - J, - tolerance=1e-2) - - # - # Simulate the explicit solution and plot the results. - # - exp_out = traj.simulate() - - if matplotlib is not None: - plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.x', - 'time (s)', 'x $(m)$'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.u', - 'time (s)', 'u $(m/s^2)$')], - title='Hyper Sensitive Problem Solution\nRadau Pseudospectral Method', - p_sol=p, p_sim=exp_out) - - plt.show() - - -if __name__ == "__main__": - unittest.main() diff --git a/dymos/examples/length_constrained_brachistochrone/doc/__init__.py b/dymos/examples/length_constrained_brachistochrone/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/balanced_field/doc/__init__.py b/dymos/examples/length_constrained_brachistochrone/test/__init__.py similarity index 100% rename from dymos/examples/balanced_field/doc/__init__.py rename to dymos/examples/length_constrained_brachistochrone/test/__init__.py diff --git a/dymos/examples/length_constrained_brachistochrone/doc/test_doc_length_constrained_brachistochrone.py b/dymos/examples/length_constrained_brachistochrone/test/test_length_constrained_brachistochrone.py similarity index 100% rename from dymos/examples/length_constrained_brachistochrone/doc/test_doc_length_constrained_brachistochrone.py rename to dymos/examples/length_constrained_brachistochrone/test/test_length_constrained_brachistochrone.py diff --git a/dymos/examples/robertson_problem/doc/__init__.py b/dymos/examples/robertson_problem/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/brachistochrone/doc/__init__.py b/dymos/examples/robertson_problem/test/__init__.py similarity index 100% rename from dymos/examples/brachistochrone/doc/__init__.py rename to dymos/examples/robertson_problem/test/__init__.py diff --git a/dymos/examples/robertson_problem/doc/test_doc_robertson_problem.py b/dymos/examples/robertson_problem/test/test_robertson_problem.py similarity index 97% rename from dymos/examples/robertson_problem/doc/test_doc_robertson_problem.py rename to dymos/examples/robertson_problem/test/test_robertson_problem.py index d4898ae3f3..040b8ab15b 100644 --- a/dymos/examples/robertson_problem/doc/test_doc_robertson_problem.py +++ b/dymos/examples/robertson_problem/test/test_robertson_problem.py @@ -15,7 +15,7 @@ @use_tempdirs -class TestRobertsonProblemForDocs(unittest.TestCase): +class TestRobertsonProblem(unittest.TestCase): def robertson_problem(self, t_final=1.0): @@ -62,7 +62,7 @@ def robertson_problem(self, t_final=1.0): @unittest.skipIf(matplotlib is None, "This test requires matplotlib") @unittest.skipIf(sys.platform == 'win32', "LSODA accuracy issues with explicit shooting on Windows") - def test_robertson_problem_for_docs(self): + def test_robertson_problem(self): import openmdao.api as om from dymos.utils.testing_utils import assert_check_partials diff --git a/dymos/examples/robot_arm/doc/__init__.py b/dymos/examples/robot_arm/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py deleted file mode 100644 index 788cb5310d..0000000000 --- a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py +++ /dev/null @@ -1,114 +0,0 @@ -import os -import unittest - -import numpy as np -import openmdao.api as om -import dymos as dm - -from openmdao.utils.assert_utils import assert_near_equal -from openmdao.utils.general_utils import printoptions -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - -from dymos.examples.robot_arm.robot_arm_ode import RobotArmODE - - -@use_tempdirs -class TestRobotArm(unittest.TestCase): - - def tearDown(self): - for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out']: - if os.path.exists(filename): - os.remove(filename) - - @require_pyoptsparse(optimizer='SLSQP') - def make_problem(self, tx, optimizer='SLSQP', numseg=30): - p = om.Problem(model=om.Group()) - p.driver = om.pyOptSparseDriver() - p.driver.declare_coloring() - p.driver.options['optimizer'] = optimizer - if optimizer == 'SNOPT': - p.driver.opt_settings['iSumm'] = 6 - p.driver.opt_settings['Verify level'] = 3 - elif optimizer == 'IPOPT': - p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based' - p.driver.opt_settings['max_iter'] = 500 - p.driver.opt_settings['print_level'] = 0 - traj = p.model.add_subsystem('traj', dm.Trajectory()) - - phase = traj.add_phase('phase', dm.Phase(ode_class=RobotArmODE, transcription=tx)) - phase.set_time_options(fix_initial=True, fix_duration=False) - - phase.add_state('x0', fix_initial=True, fix_final=True, rate_source='x0_dot', units='m') - phase.add_state('x1', fix_initial=True, fix_final=True, rate_source='x1_dot', units='rad') - phase.add_state('x2', fix_initial=True, fix_final=True, rate_source='x2_dot', units='rad') - phase.add_state('x3', fix_initial=True, fix_final=True, rate_source='x3_dot', units='m/s') - phase.add_state('x4', fix_initial=True, fix_final=True, rate_source='x4_dot', units='rad/s') - phase.add_state('x5', fix_initial=True, fix_final=True, rate_source='x5_dot', units='rad/s') - - phase.add_control('u0', opt=True, lower=-1, upper=1, scaler=0.1, units='m**2/s**2', - continuity=False, rate_continuity=False) - phase.add_control('u1', opt=True, lower=-1, upper=1, scaler=0.1, units='m**3*rad/s**2', - continuity=False, rate_continuity=False) - phase.add_control('u2', opt=True, lower=-1, upper=1, scaler=0.1, units='m**3*rad/s**2', - continuity=False, rate_continuity=False) - - phase.add_path_constraint('u0', lower=-1, upper=1, scaler=0.1) - phase.add_path_constraint('u1', lower=-1, upper=1, scaler=0.1) - phase.add_path_constraint('u2', lower=-1, upper=1, scaler=0.1) - - phase.add_objective('time', ref=0.1) - - phase.set_refine_options(refine=True, tol=1e-5, smoothness_factor=1.2) - - p.setup(check=True, force_alloc_complex=False, mode='auto') - - phase.set_time_val(initial=0, duration=10) - phase.set_state_val('x0', [4.5, 4.5]) - phase.set_state_val('x1', [0.0, 2 * np.pi / 3]) - phase.set_state_val('x2', [np.pi / 4, np.pi / 4]) - phase.set_state_val('x3', [0.0, 0.0]) - phase.set_state_val('x4', [0.0, 0.0]) - phase.set_state_val('x5', [0.0, 0.0]) - - return p - - def test_partials(self): - p = self.make_problem(tx=dm.Radau(num_segments=30, order=3), optimizer='SLSQP') - p.run_model() - with printoptions(linewidth=1024, edgeitems=100): - p.check_partials(method='fd', compact_print=True, out_stream=None) - - @require_pyoptsparse(optimizer='IPOPT') - def test_robot_arm_radau(self): - p = self.make_problem(tx=dm.Radau(num_segments=12, order=3), optimizer='IPOPT', numseg=12) - dm.run_problem(p, refine_iteration_limit=5) - - t = p.get_val('traj.phase.timeseries.time') - assert_near_equal(t[-1], 9.14138, tolerance=1e-3) - - @require_pyoptsparse(optimizer='IPOPT') - def test_robot_arm_gl(self): - p = self.make_problem(tx=dm.GaussLobatto(num_segments=20, order=3), optimizer='IPOPT') - dm.run_problem(p, refine_iteration_limit=3) - - t = p.get_val('traj.phase.timeseries.time') - assert_near_equal(t[-1], 9.14138, tolerance=1e-3) - - @require_pyoptsparse(optimizer='IPOPT') - def test_robot_arm_birkhoff_lgl(self): - tx = dm.Birkhoff(num_nodes=30, grid_type='lgl') - p = self.make_problem(tx=tx, optimizer='IPOPT') - dm.run_problem(p) - - t = p.get_val('traj.phase.timeseries.time') - assert_near_equal(t[-1], 9.14138, tolerance=1e-3) - - @require_pyoptsparse(optimizer='IPOPT') - def test_robot_arm_birkhoff_cgl(self): - tx = dm.Birkhoff(num_nodes=30, grid_type='cgl') - p = self.make_problem(tx=tx, optimizer='IPOPT') - dm.run_problem(p) - - t = p.get_val('traj.phase.timeseries.time') - - assert_near_equal(t[-1], 9.14138, tolerance=1e-3) diff --git a/dymos/examples/shuttle_reentry/doc/__init__.py b/dymos/examples/shuttle_reentry/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/shuttle_reentry/doc/test_doc_reentry.py b/dymos/examples/shuttle_reentry/doc/test_doc_reentry.py deleted file mode 100644 index 96e52add2f..0000000000 --- a/dymos/examples/shuttle_reentry/doc/test_doc_reentry.py +++ /dev/null @@ -1,110 +0,0 @@ -import os -import unittest - -try: - import matplotlib.pyplot as plt - plt.switch_backend('Agg') - plt.style.use('ggplot') -except ImportError: - plt = None - -import numpy as np - -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -@use_tempdirs -class TestReentryForDocs(unittest.TestCase): - - def tearDown(self): - for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']: - if os.path.exists(filename): - os.remove(filename) - - @require_pyoptsparse(optimizer='SLSQP') - @unittest.skipIf(plt is None, "This test requires matplotlib") - def test_reentry(self): - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - from dymos.examples.shuttle_reentry.shuttle_ode import ShuttleODE - from dymos.examples.plotting import plot_results - - # Instantiate the problem, add the driver, and allow it to use coloring - p = om.Problem(model=om.Group()) - p.driver = om.pyOptSparseDriver() - p.driver.declare_coloring() - p.driver.options['optimizer'] = 'SLSQP' - - # Instantiate the trajectory and add a phase to it - traj = p.model.add_subsystem('traj', dm.Trajectory()) - phase0 = traj.add_phase('phase0', - dm.Phase(ode_class=ShuttleODE, - transcription=dm.Radau(num_segments=15, order=3))) - - phase0.set_time_options(fix_initial=True, units='s', duration_ref=200) - phase0.add_state('h', fix_initial=True, fix_final=True, units='ft', rate_source='hdot', - lower=0, ref0=75000, ref=300000, defect_ref=1000) - phase0.add_state('gamma', fix_initial=True, fix_final=True, units='rad', - rate_source='gammadot', - lower=-89. * np.pi / 180, upper=89. * np.pi / 180) - phase0.add_state('phi', fix_initial=True, fix_final=False, units='rad', - rate_source='phidot', lower=0, upper=89. * np.pi / 180) - phase0.add_state('psi', fix_initial=True, fix_final=False, units='rad', - rate_source='psidot', lower=0, upper=90. * np.pi / 180) - phase0.add_state('theta', fix_initial=True, fix_final=False, units='rad', - rate_source='thetadot', - lower=-89. * np.pi / 180, upper=89. * np.pi / 180) - phase0.add_state('v', fix_initial=True, fix_final=True, units='ft/s', - rate_source='vdot', lower=0, ref0=2500, ref=25000) - phase0.add_control('alpha', units='rad', opt=True, lower=-np.pi / 2, upper=np.pi / 2, ) - phase0.add_control('beta', units='rad', opt=True, lower=-89 * np.pi / 180, upper=1 * np.pi / 180, ) - - # The original implementation by Betts includes a heating rate path constraint. - # This will work with the SNOPT optimizer but SLSQP has difficulty converging the solution. - # phase0.add_path_constraint('q', lower=0, upper=70, ref=70) - phase0.add_timeseries_output('q', shape=(1,)) - - phase0.add_objective('theta', loc='final', ref=-0.01) - - p.setup(check=True) - - phase0.set_time_val(initial=0, duration=2000, units='s') - phase0.set_state_val('h', [260000, 80000], units='ft') - phase0.set_state_val('gamma', [-1, -5], units='deg') - phase0.set_state_val('phi', [0, 75], units='deg') - phase0.set_state_val('psi', [90, 10], units='deg') - phase0.set_state_val('theta', [0, 25], units='deg') - phase0.set_state_val('v', [25600, 2500], units='ft/s') - phase0.set_control_val('alpha', 17.4, units='deg') - phase0.set_control_val('beta', [-75, 0], units='deg') - - # Run the driver - dm.run_problem(p) - - # Check the validity of the solution - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 2008.59, - tolerance=1e-3) - assert_near_equal(p.get_val('traj.phase0.timeseries.theta', units='deg')[-1], - 34.1412, tolerance=1e-3) - - # Run the simulation to check if the model is physically valid - sim_out = traj.simulate() - - # Plot the results - - plot_results([('traj.phase0.timeseries.time', 'traj.phase0.timeseries.alpha', - 'time (s)', 'alpha (rad)'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.beta', - 'time (s)', 'beta (rad)'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.theta', - 'time (s)', 'theta (rad)'), - ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.q', - 'time (s)', 'q (btu/ft/ft/s')], title='Reentry Solution', p_sol=p, - p_sim=sim_out) - - plt.show() - - -if __name__ == '__main__': # pragma: no cover - unittest.main() diff --git a/dymos/examples/ssto/doc/__init__.py b/dymos/examples/ssto/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/ssto/doc/test_doc_ssto_earth.py b/dymos/examples/ssto/doc/test_doc_ssto_earth.py deleted file mode 100644 index 3ade491c93..0000000000 --- a/dymos/examples/ssto/doc/test_doc_ssto_earth.py +++ /dev/null @@ -1,142 +0,0 @@ -import unittest - -try: - import matplotlib -except ImportError: - matplotlib = None - -from openmdao.utils.assert_utils import assert_near_equal -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -@use_tempdirs -class TestDocSSTOEarth(unittest.TestCase): - - @require_pyoptsparse(optimizer='SLSQP') - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_doc_ssto_earth(self): - import matplotlib.pyplot as plt - 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, targets=['theta']) - phase.add_parameter('thrust', units='N', opt=False, val=2100000.0, targets=['thrust']) - - # - # 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) - - phase.set_time_val(initial=0.0, duration=150.0) - phase.set_state_val('x', [0, 1.15E5]) - phase.set_state_val('y', [0, 1.85E5]) - phase.set_state_val('vy', [1.0E-6, 0]) - phase.set_state_val('m', [117000, 1163]) - phase.set_control_val('theta', [1.5, -0.76]) - phase.set_parameter_val('thrust', 2.1, units='MN') - - # - # Solve the Problem - # - dm.run_problem(p) - - 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) - # - # Get the explicitly simulated results - # - exp_out = traj.simulate() - - # - # Plot the results - # - fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 8)) - - axes[0].plot(p.get_val('traj.phase0.timeseries.x'), - p.get_val('traj.phase0.timeseries.y'), - marker='o', - ms=4, - linestyle='None', - label='solution') - - axes[0].plot(exp_out.get_val('traj.phase0.timeseries.x'), - exp_out.get_val('traj.phase0.timeseries.y'), - marker=None, - linestyle='-', - label='simulation') - - axes[0].set_xlabel('range (m)') - axes[0].set_ylabel('altitude (m)') - axes[0].set_aspect('equal') - - axes[1].plot(p.get_val('traj.phase0.timeseries.time'), - p.get_val('traj.phase0.timeseries.theta'), - marker='o', - ms=4, - linestyle='None') - - axes[1].plot(exp_out.get_val('traj.phase0.timeseries.time'), - exp_out.get_val('traj.phase0.timeseries.theta'), - linestyle='-', - marker=None) - - axes[1].set_xlabel('time (s)') - axes[1].set_ylabel('theta (deg)') - - plt.suptitle('Single Stage to Orbit Solution Using Linear Tangent Guidance') - fig.legend(loc='lower center', ncol=2) - - plt.show() - - -if __name__ == "__main__": - unittest.main() diff --git a/dymos/examples/ssto/doc/test_doc_ssto_linear_tangent_guidance.py b/dymos/examples/ssto/doc/test_doc_ssto_linear_tangent_guidance.py deleted file mode 100644 index 12308adad2..0000000000 --- a/dymos/examples/ssto/doc/test_doc_ssto_linear_tangent_guidance.py +++ /dev/null @@ -1,273 +0,0 @@ -import unittest - -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -@use_tempdirs -class TestDocSSTOLinearTangentGuidance(unittest.TestCase): - - @require_pyoptsparse(optimizer='SLSQP') - def test_doc_ssto_linear_tangent_guidance(self): - import numpy as np - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - - g = 1.61544 # lunar gravity, m/s**2 - - class LaunchVehicle2DEOM(om.ExplicitComponent): - """ - Simple 2D Cartesian Equations of Motion for a launch vehicle subject to thrust and drag. - """ - def initialize(self): - self.options.declare('num_nodes', types=int) - - def setup(self): - nn = self.options['num_nodes'] - - # Inputs - self.add_input('vx', - val=np.zeros(nn), - desc='x velocity', - units='m/s') - - self.add_input('vy', - val=np.zeros(nn), - desc='y velocity', - units='m/s') - - self.add_input('m', - val=np.zeros(nn), - desc='mass', - units='kg') - - self.add_input('theta', - val=np.zeros(nn), - desc='pitch angle', - units='rad') - - self.add_input('thrust', - val=2100000 * np.ones(nn), - desc='thrust', - units='N') - - self.add_input('Isp', - val=265.2 * np.ones(nn), - desc='specific impulse', - units='s') - - # Outputs - 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('vxdot', - val=np.zeros(nn), - desc='x acceleration magnitude', - units='m/s**2') - - self.add_output('vydot', - val=np.zeros(nn), - desc='y acceleration magnitude', - units='m/s**2') - - self.add_output('mdot', - val=np.zeros(nn), - desc='mass rate of change', - units='kg/s') - - # Setup partials - ar = np.arange(self.options['num_nodes']) - - self.declare_partials(of='xdot', wrt='vx', rows=ar, cols=ar, val=1.0) - self.declare_partials(of='ydot', wrt='vy', rows=ar, cols=ar, val=1.0) - - self.declare_partials(of='vxdot', wrt='vx', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='m', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='theta', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='thrust', rows=ar, cols=ar) - - self.declare_partials(of='vydot', wrt='m', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='theta', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='vy', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='thrust', rows=ar, cols=ar) - - self.declare_partials(of='mdot', wrt='thrust', rows=ar, cols=ar) - self.declare_partials(of='mdot', wrt='Isp', rows=ar, cols=ar) - - def compute(self, inputs, outputs): - theta = inputs['theta'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - vx = inputs['vx'] - vy = inputs['vy'] - m = inputs['m'] - F_T = inputs['thrust'] - Isp = inputs['Isp'] - - outputs['xdot'] = vx - outputs['ydot'] = vy - outputs['vxdot'] = F_T * cos_theta / m - outputs['vydot'] = F_T * sin_theta / m - g - outputs['mdot'] = -F_T / (g * Isp) - - def compute_partials(self, inputs, jacobian): - theta = inputs['theta'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - m = inputs['m'] - F_T = inputs['thrust'] - Isp = inputs['Isp'] - - # jacobian['vxdot', 'vx'] = -CDA * rho * vx / m - jacobian['vxdot', 'm'] = -(F_T * cos_theta) / m ** 2 - jacobian['vxdot', 'theta'] = -(F_T / m) * sin_theta - jacobian['vxdot', 'thrust'] = cos_theta / m - - # jacobian['vydot', 'vy'] = -CDA * rho * vy / m - jacobian['vydot', 'm'] = -(F_T * sin_theta) / m ** 2 - jacobian['vydot', 'theta'] = (F_T / m) * cos_theta - jacobian['vydot', 'thrust'] = sin_theta / m - - jacobian['mdot', 'thrust'] = -1.0 / (g * Isp) - jacobian['mdot', 'Isp'] = F_T / (g * Isp ** 2) - - class LinearTangentGuidanceComp(om.ExplicitComponent): - """ Compute pitch angle from static controls governing linear expression for - pitch angle tangent as function of time. - """ - - def initialize(self): - self.options.declare('num_nodes', types=int) - - def setup(self): - nn = self.options['num_nodes'] - - self.add_input('a_ctrl', - val=np.zeros(nn), - desc='linear tangent slope', - units='1/s') - - self.add_input('b_ctrl', - val=np.zeros(nn), - desc='tangent of theta at t=0', - units=None) - - self.add_input('time_phase', - val=np.zeros(nn), - desc='time', - units='s') - - self.add_output('theta', - val=np.zeros(nn), - desc='pitch angle', - units='rad') - - # Setup partials - arange = np.arange(self.options['num_nodes']) - - self.declare_partials(of='theta', wrt='a_ctrl', rows=arange, cols=arange, val=1.0) - self.declare_partials(of='theta', wrt='b_ctrl', rows=arange, cols=arange, val=1.0) - self.declare_partials(of='theta', wrt='time_phase', rows=arange, cols=arange, val=1.0) - - def compute(self, inputs, outputs): - a = inputs['a_ctrl'] - b = inputs['b_ctrl'] - t = inputs['time_phase'] - outputs['theta'] = np.arctan(a * t + b) - - def compute_partials(self, inputs, jacobian): - a = inputs['a_ctrl'] - b = inputs['b_ctrl'] - t = inputs['time_phase'] - - x = a * t + b - denom = x ** 2 + 1.0 - - jacobian['theta', 'a_ctrl'] = t / denom - jacobian['theta', 'b_ctrl'] = 1.0 / denom - jacobian['theta', 'time_phase'] = a / denom - - class LaunchVehicleLinearTangentODE(om.Group): - """ - The LaunchVehicleLinearTangentODE for this case consists of a guidance component and - the EOM. Guidance is simply an OpenMDAO ExecComp which computes the arctangent of the - tan_theta variable. - """ - def initialize(self): - self.options.declare('num_nodes', types=int, - desc='Number of nodes to be evaluated in the RHS') - - def setup(self): - nn = self.options['num_nodes'] - self.add_subsystem('guidance', LinearTangentGuidanceComp(num_nodes=nn)) - self.add_subsystem('eom', LaunchVehicle2DEOM(num_nodes=nn)) - self.connect('guidance.theta', 'eom.theta') - - # - # Setup and solve the optimal control problem - # - p = om.Problem(model=om.Group()) - - p.driver = om.pyOptSparseDriver() - p.driver.declare_coloring() - - traj = dm.Trajectory() - p.model.add_subsystem('traj', traj) - - phase = dm.Phase(ode_class=LaunchVehicleLinearTangentODE, - transcription=dm.GaussLobatto(num_segments=10, order=5, compressed=True)) - - traj.add_phase('phase0', phase) - - phase.set_time_options(fix_initial=True, duration_bounds=(10, 1000), - targets=['guidance.time_phase']) - - phase.add_state('x', fix_initial=True, lower=0, rate_source='eom.xdot', units='m') - phase.add_state('y', fix_initial=True, lower=0, rate_source='eom.ydot', units='m') - phase.add_state('vx', fix_initial=True, lower=0, rate_source='eom.vxdot', targets=['eom.vx'], units='m/s') - phase.add_state('vy', fix_initial=True, rate_source='eom.vydot', targets=['eom.vy'], units='m/s') - phase.add_state('m', fix_initial=True, rate_source='eom.mdot', targets=['eom.m'], units='kg') - - phase.add_boundary_constraint('y', loc='final', equals=1.85E5, linear=True) - phase.add_boundary_constraint('vx', loc='final', equals=1627.0) - phase.add_boundary_constraint('vy', loc='final', equals=0) - - phase.add_parameter('a_ctrl', units='1/s', opt=True, targets=['guidance.a_ctrl']) - phase.add_parameter('b_ctrl', units=None, opt=True, targets=['guidance.b_ctrl']) - phase.add_parameter('thrust', units='N', opt=False, val=3.0 * 50000.0 * 1.61544, targets=['eom.thrust']) - phase.add_parameter('Isp', units='s', opt=False, val=1.0E6, targets=['eom.Isp']) - - phase.add_objective('time', index=-1, scaler=0.01) - - p.model.linear_solver = om.DirectSolver() - - phase.add_timeseries_output('guidance.theta', units='deg') - - p.setup(check=True) - - phase.set_time_val(initial=0.0, duration=500.0) - phase.set_state_val('x', [0, 350000.0]) - phase.set_state_val('y', [0, 185000.0]) - phase.set_state_val('vx', [0, 1627.0]) - phase.set_state_val('vy', [1.0E-6, 0.0]) - phase.set_state_val('m', 50000) - phase.set_parameter_val('a_ctrl', -0.01) - phase.set_parameter_val('b_ctrl', 3.0) - - dm.run_problem(p, simulate=True) - - # - # Check the results. - # - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 481, tolerance=0.01) - - -if __name__ == "__main__": - unittest.main() diff --git a/dymos/examples/ssto/doc/test_doc_ssto_polynomial_control.py b/dymos/examples/ssto/doc/test_doc_ssto_polynomial_control.py deleted file mode 100644 index d106662a18..0000000000 --- a/dymos/examples/ssto/doc/test_doc_ssto_polynomial_control.py +++ /dev/null @@ -1,317 +0,0 @@ -import unittest - -try: - import matplotlib - import matplotlib.pyplot as plt - - matplotlib.use('Agg') - plt.style.use('ggplot') -except ImportError: - matplotlib = None - -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse - - -@use_tempdirs -class TestDocSSTOPolynomialControl(unittest.TestCase): - - @require_pyoptsparse(optimizer='SLSQP') - @unittest.skipIf(matplotlib is None, "This test requires matplotlib") - def test_doc_ssto_polynomial_control(self): - import numpy as np - import matplotlib.pyplot as plt - import openmdao.api as om - from openmdao.utils.assert_utils import assert_near_equal - import dymos as dm - - g = 1.61544 # lunar gravity, m/s**2 - - class LaunchVehicle2DEOM(om.ExplicitComponent): - """ - Simple 2D Cartesian Equations of Motion for a launch vehicle subject to thrust and drag. - """ - def initialize(self): - self.options.declare('num_nodes', types=int) - - def setup(self): - nn = self.options['num_nodes'] - - # Inputs - self.add_input('vx', - val=np.zeros(nn), - desc='x velocity', - units='m/s') - - self.add_input('vy', - val=np.zeros(nn), - desc='y velocity', - units='m/s') - - self.add_input('m', - val=np.zeros(nn), - desc='mass', - units='kg') - - self.add_input('theta', - val=np.zeros(nn), - desc='pitch angle', - units='rad') - - self.add_input('thrust', - val=2100000 * np.ones(nn), - desc='thrust', - units='N') - - self.add_input('Isp', - val=265.2 * np.ones(nn), - desc='specific impulse', - units='s') - - # Outputs - 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('vxdot', - val=np.zeros(nn), - desc='x acceleration magnitude', - units='m/s**2') - - self.add_output('vydot', - val=np.zeros(nn), - desc='y acceleration magnitude', - units='m/s**2') - - self.add_output('mdot', - val=np.zeros(nn), - desc='mass rate of change', - units='kg/s') - - # Setup partials - ar = np.arange(self.options['num_nodes']) - - self.declare_partials(of='xdot', wrt='vx', rows=ar, cols=ar, val=1.0) - self.declare_partials(of='ydot', wrt='vy', rows=ar, cols=ar, val=1.0) - - self.declare_partials(of='vxdot', wrt='vx', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='m', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='theta', rows=ar, cols=ar) - self.declare_partials(of='vxdot', wrt='thrust', rows=ar, cols=ar) - - self.declare_partials(of='vydot', wrt='m', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='theta', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='vy', rows=ar, cols=ar) - self.declare_partials(of='vydot', wrt='thrust', rows=ar, cols=ar) - - self.declare_partials(of='mdot', wrt='thrust', rows=ar, cols=ar) - self.declare_partials(of='mdot', wrt='Isp', rows=ar, cols=ar) - - def compute(self, inputs, outputs): - theta = inputs['theta'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - vx = inputs['vx'] - vy = inputs['vy'] - m = inputs['m'] - F_T = inputs['thrust'] - Isp = inputs['Isp'] - - outputs['xdot'] = vx - outputs['ydot'] = vy - outputs['vxdot'] = F_T * cos_theta / m - outputs['vydot'] = F_T * sin_theta / m - g - outputs['mdot'] = -F_T / (g * Isp) - - def compute_partials(self, inputs, jacobian): - theta = inputs['theta'] - cos_theta = np.cos(theta) - sin_theta = np.sin(theta) - m = inputs['m'] - F_T = inputs['thrust'] - Isp = inputs['Isp'] - - # jacobian['vxdot', 'vx'] = -CDA * rho * vx / m - jacobian['vxdot', 'm'] = -(F_T * cos_theta) / m ** 2 - jacobian['vxdot', 'theta'] = -(F_T / m) * sin_theta - jacobian['vxdot', 'thrust'] = cos_theta / m - - # jacobian['vydot', 'vy'] = -CDA * rho * vy / m - jacobian['vydot', 'm'] = -(F_T * sin_theta) / m ** 2 - jacobian['vydot', 'theta'] = (F_T / m) * cos_theta - jacobian['vydot', 'thrust'] = sin_theta / m - - jacobian['mdot', 'thrust'] = -1.0 / (g * Isp) - jacobian['mdot', 'Isp'] = F_T / (g * Isp ** 2) - - class LaunchVehicleLinearTangentODE(om.Group): - """ - The LaunchVehicleLinearTangentODE for this case consists of a guidance component and - the EOM. Guidance is simply an OpenMDAO ExecComp which computes the arctangent of the - tan_theta variable. - """ - - def initialize(self): - self.options.declare('num_nodes', types=int, - desc='Number of nodes to be evaluated in the RHS') - - def setup(self): - nn = self.options['num_nodes'] - - self.add_subsystem('guidance', om.ExecComp('theta=arctan(tan_theta)', - theta={'val': np.ones(nn), - 'units': 'rad'}, - tan_theta={'val': np.ones(nn)})) - - self.add_subsystem('eom', LaunchVehicle2DEOM(num_nodes=nn)) - - self.connect('guidance.theta', 'eom.theta') - - # - # Setup and solve the optimal control problem - # - p = om.Problem(model=om.Group()) - - traj = p.model.add_subsystem('traj', dm.Trajectory()) - - phase = dm.Phase(ode_class=LaunchVehicleLinearTangentODE, - transcription=dm.Radau(num_segments=20, order=3, compressed=False)) - traj.add_phase('phase0', phase) - - phase.set_time_options(fix_initial=True, duration_bounds=(10, 1000), units='s') - - # - # Set the state options. We include rate_source, units, and targets here since the ODE - # is not decorated with their default values. - # - phase.add_state('x', fix_initial=True, lower=0, rate_source='eom.xdot', units='m') - phase.add_state('y', fix_initial=True, lower=0, rate_source='eom.ydot', units='m') - phase.add_state('vx', fix_initial=True, lower=0, rate_source='eom.vxdot', - units='m/s', targets=['eom.vx']) - phase.add_state('vy', fix_initial=True, rate_source='eom.vydot', - units='m/s', targets=['eom.vy']) - phase.add_state('m', fix_initial=True, rate_source='eom.mdot', - units='kg', targets=['eom.m']) - - # - # The tangent of theta is modeled as a linear polynomial over the duration of the phase. - # - phase.add_control('tan_theta', order=1, units=None, opt=True, - targets=['guidance.tan_theta'], control_type='polynomial') - - # - # Parameters values for thrust and specific impulse are design parameters. They are - # provided by an IndepVarComp in the phase, but with opt=False their values are not - # design variables in the optimization problem. - # - phase.add_parameter('thrust', units='N', opt=False, val=3.0 * 50000.0 * 1.61544, - targets=['eom.thrust']) - phase.add_parameter('Isp', units='s', opt=False, val=1.0E6, targets=['eom.Isp']) - - # - # Set the boundary constraints. These are all states which could also be handled - # by setting fix_final=True and including the correct final value in the initial guess. - # - phase.add_boundary_constraint('y', loc='final', equals=1.85E5, linear=True) - phase.add_boundary_constraint('vx', loc='final', equals=1627.0) - phase.add_boundary_constraint('vy', loc='final', equals=0) - - phase.add_objective('time', index=-1, scaler=0.01) - - # - # Add theta as a timeseries output since it's not included by default. - # - phase.add_timeseries_output('guidance.theta', units='deg') - - # - # Set the optimizer - # - p.driver = om.pyOptSparseDriver() - p.driver.options['optimizer'] = 'SLSQP' - p.driver.declare_coloring() - - # - # We don't strictly need to define a linear solver here since our problem is entirely - # feed-forward with no iterative loops. It's good practice to add one, however, since - # failing to do so can cause incorrect derivatives if iterative processes are ever - # introduced to the system. - # - p.model.linear_solver = om.DirectSolver() - - p.setup(check=True) - - # - # Assign initial guesses for the independent variables in the problem. - # - phase.set_time_val(initial=0.0, duration=500.0) - phase.set_state_val('x', [0, 350000.0]) - phase.set_state_val('y', [0, 185000.0]) - phase.set_state_val('vx', [0, 1627.0]) - phase.set_state_val('vy', [1.0E-6, 0.0]) - phase.set_state_val('m', 50000) - phase.set_control_val('tan_theta', [[0.5 * np.pi], [0.0]]) - - # - # Solve the problem. - # - dm.run_problem(p) - - # - # Check the results. - # - assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 481, tolerance=0.01) - - # - # Get the explitly simulated results - # - exp_out = traj.simulate() - - # - # Plot the results - # - fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 8)) - - axes[0].plot(p.get_val('traj.phase0.timeseries.x'), - p.get_val('traj.phase0.timeseries.y'), - marker='o', - ms=4, - linestyle='None', - label='solution') - - axes[0].plot(exp_out.get_val('traj.phase0.timeseries.x'), - exp_out.get_val('traj.phase0.timeseries.y'), - marker=None, - linestyle='-', - label='simulation') - - axes[0].set_xlabel('range (m)') - axes[0].set_ylabel('altitude (m)') - axes[0].set_aspect('equal') - - axes[1].plot(p.get_val('traj.phase0.timeseries.time'), - p.get_val('traj.phase0.timeseries.theta'), - marker='o', - ms=4, - linestyle='None') - - axes[1].plot(exp_out.get_val('traj.phase0.timeseries.time'), - exp_out.get_val('traj.phase0.timeseries.theta'), - linestyle='-', - marker=None) - - axes[1].set_xlabel('time (s)') - axes[1].set_ylabel('theta (deg)') - - plt.suptitle('Single Stage to Orbit Solution Using Polynomial Controls') - fig.legend(loc='lower center', ncol=2) - - plt.show() - - -if __name__ == "__main__": - unittest.main() diff --git a/dymos/examples/vanderpol/doc/__init__.py b/dymos/examples/vanderpol/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py b/dymos/examples/vanderpol/doc/test_doc_vanderpol.py deleted file mode 100644 index f0b527468a..0000000000 --- a/dymos/examples/vanderpol/doc/test_doc_vanderpol.py +++ /dev/null @@ -1,132 +0,0 @@ -import os -import unittest - -from openmdao.utils.testing_utils import use_tempdirs -from openmdao.utils.mpi import MPI - - -@use_tempdirs -class TestVanderpolForDocs(unittest.TestCase): - def tearDown(self): - for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']: - if os.path.exists(filename): - os.remove(filename) - - def test_vanderpol_for_docs_simulation(self): - import dymos as dm - from dymos.examples.vanderpol.vanderpol_dymos import vanderpol - - # Create the Dymos problem instance - p = vanderpol(transcription='gauss-lobatto', num_segments=75) - - dm.run_problem(p, run_driver=False, simulate=True, make_plots=True) - - def test_vanderpol_for_docs_optimize(self): - import dymos as dm - from dymos.examples.vanderpol.vanderpol_dymos import vanderpol - - # Create the Dymos problem instance - p = vanderpol(transcription='gauss-lobatto', num_segments=75, - transcription_order=3, compressed=True, optimizer='SLSQP') - - dm.run_problem(p, simulate=True, make_plots=True) - - def test_vanderpol_for_docs_optimize_refine(self): - import dymos as dm - from dymos.examples.vanderpol.vanderpol_dymos import vanderpol - from openmdao.utils.assert_utils import assert_near_equal - - # Create the Dymos problem instance - p = vanderpol(transcription='gauss-lobatto', num_segments=15, - transcription_order=3, compressed=True, optimizer='SLSQP') - - # Enable grid refinement and find optimal control solution to stop oscillation - p.model.traj.phases.phase0.set_refine_options(refine=True) - - dm.run_problem(p, refine_iteration_limit=10, simulate=True, make_plots=True) - - assert_near_equal(p.get_val('traj.phase0.states:x0')[-1, ...], 0.0) - assert_near_equal(p.get_val('traj.phase0.states:x1')[-1, ...], 0.0) - assert_near_equal(p.get_val('traj.phase0.states:J')[-1, ...], 5.2808, tolerance=1.0E-3) - assert_near_equal(p.get_val('traj.phase0.controls:u')[-1, ...], 0.0, tolerance=1.0E-3) - - -@unittest.skipUnless(MPI, "MPI is required.") -@use_tempdirs -class TestVanderpolDelayMPI(unittest.TestCase): - N_PROCS = 2 - - def test_vanderpol_delay_mpi(self): - import openmdao.api as om - import dymos as dm - from dymos.examples.vanderpol.vanderpol_ode import VanderpolODE - from openmdao.utils.assert_utils import assert_near_equal - - DELAY = 0.005 - - p = om.Problem(model=om.Group()) - p.driver = om.ScipyOptimizeDriver() - - p.driver.options['optimizer'] = 'SLSQP' - p.driver.declare_coloring() - - # define a Trajectory object and add to model - traj = dm.Trajectory() - p.model.add_subsystem('traj', subsys=traj) - - t = dm.Radau(num_segments=30, order=3) - - # define a Phase as specified above and add to Phase - phase = dm.Phase(ode_class=VanderpolODE, transcription=t, - ode_init_kwargs={'delay': DELAY, 'distrib': True}) - traj.add_phase(name='phase0', phase=phase) - - t_final = 15 - phase.set_time_options(fix_initial=True, fix_duration=True, duration_val=t_final, units='s') - - # set the State time options - phase.add_state('x0', fix_initial=False, fix_final=False, - rate_source='x0dot', - units='V/s', - targets='x0') # target required because x0 is an input - phase.add_state('x1', fix_initial=False, fix_final=False, - rate_source='x1dot', - units='V', - targets='x1') # target required because x1 is an input - phase.add_state('J', fix_initial=False, fix_final=False, - rate_source='Jdot', - units=None) - - # define the control - phase.add_control(name='u', units=None, lower=-0.75, upper=1.0, continuity=True, - rate_continuity=True, - targets='u') # target required because u is an input - - # add constraints - phase.add_boundary_constraint('x0', loc='initial', equals=1.0) - phase.add_boundary_constraint('x1', loc='initial', equals=1.0) - phase.add_boundary_constraint('J', loc='initial', equals=0.0) - - phase.add_boundary_constraint('x0', loc='final', equals=0.0) - phase.add_boundary_constraint('x1', loc='final', equals=0.0) - - # define objective to minimize - phase.add_objective('J', loc='final') - - # setup the problem - p.setup(check=True) - - phase.set_time_val(0.0, t_final) - - # add a linearly interpolated initial guess for the state and control curves - phase.set_state_val('x0', [1, 0]) - phase.set_state_val('x0', [1, 0]) - phase.set_state_val('x0', [0, 1]) - phase.set_control_val('u', -0.75) - - dm.run_problem(p, run_driver=True, simulate=False) - - assert_near_equal(p.get_val('traj.phase0.states:x0')[-1, ...], 0.0) - assert_near_equal(p.get_val('traj.phase0.states:x1')[-1, ...], 0.0) - assert_near_equal(p.get_val('traj.phase0.states:J')[-1, ...], 5.2808, tolerance=0.0025) - assert_near_equal(p.get_val('traj.phase0.controls:u')[-1, ...], 0.0, tolerance=1.0E-3) From 806f4bc743195d1830903ced91d2880a7a0cc6be Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 11 May 2026 14:29:25 -0400 Subject: [PATCH 2/2] cleanup references to doc folders in the docs. --- .../contributing/contributing.ipynb | 20 +++++-------- .../brachistochrone_upstream_controls.ipynb | 2 +- .../brachistochrone_upstream_states.ipynb | 30 +++++++++---------- .../min_time_climb/min_time_climb.ipynb | 8 ++--- docs/dymos_book/faq/tandem_phases.ipynb | 10 +++---- .../dymos_book/faq/use_partial_coloring.ipynb | 10 +++---- .../scripts/fitter_animation_1_seg.py | 2 +- dymos/examples/min_time_climb/doc/__init__.py | 0 .../{doc => test}/aero_partial_coloring.py | 0 .../dynamic_pressure_comp_partial_coloring.py | 0 .../min_time_climb_ode_partial_coloring.py | 0 ...est_doc_min_time_climb_partial_coloring.py | 2 +- dymos/transcriptions/test/test_set_vals.py | 2 +- 13 files changed, 40 insertions(+), 46 deletions(-) delete mode 100644 dymos/examples/min_time_climb/doc/__init__.py rename dymos/examples/min_time_climb/{doc => test}/aero_partial_coloring.py (100%) rename dymos/examples/min_time_climb/{doc => test}/dynamic_pressure_comp_partial_coloring.py (100%) rename dymos/examples/min_time_climb/{doc => test}/min_time_climb_ode_partial_coloring.py (100%) diff --git a/docs/dymos_book/contributing/contributing.ipynb b/docs/dymos_book/contributing/contributing.ipynb index 76bed17bc6..fe0fdcc05e 100644 --- a/docs/dymos_book/contributing/contributing.ipynb +++ b/docs/dymos_book/contributing/contributing.ipynb @@ -78,13 +78,7 @@ "Adding a new example is a great way to contribute to Dymos.\n", "It's a great introduction to the Dymos development process, and examples provide a great way for users to learn to apply Dymos in new applications.\n", "Submit new examples via [the Dymos issues on Github](https://github.com/OpenMDAO/dymos/issues).\n", - "A new example should do the following:\n", - "\n", - "- Include a new directory under the `dymos/examples` directory.\n", - "- A unittest should be included in a `doc` subfolder within the example directory.\n", - "- The unittest method should be self-contained (it should include all imports necessary to run the example).\n", - "- If you want to include output and/or plots from the example in the documentation (highly recommended), decorate the test with the `@dymos.utils.doc_utils.save_for_docs` decorator. This will save the text and plot outputs from the test for inclusion in the Dymos documentation.\n", - "- A new markdown file should be added under `mkdocs/docs/examples/` within the Dymos repository.\n", + "A new example should include a new documentation file with the code necessary to implement the example.\n", "\n", "The Dymos docs are built on [JupyterBook](https://jupyterbook.org/intro.html) which allows users to run any page of the documentation by opening it in colab as a [Jupyter Notebook](https://jupyter.org). For those wanting to contribute, they are able to contribute by writing their own Jupyter Notebooks. Below are some important ways on how to build notebooks for Dymos.\n", "\n", @@ -137,7 +131,7 @@ }, "outputs": [], "source": [ - "om.display_source(\"dymos.examples.brachistochrone.doc.brachistochrone_ode\")" + "om.display_source(\"dymos.examples.brachistochrone.brachistochrone_ode\")" ] }, { @@ -148,7 +142,7 @@ "source": [ "import numpy as np\n", "import openmdao.api as om\n", - "from dymos.examples.brachistochrone.doc.brachistochrone_ode import BrachistochroneODE\n", + "from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE\n", "\n", "num_nodes = 5\n", "\n", @@ -165,8 +159,8 @@ "\n", "p.setup(force_alloc_complex=True)\n", "\n", - "p.set_val('vars.v', 10*np.random.random(num_nodes))\n", - "p.set_val('vars.theta', 10*np.random.uniform(1, 179, num_nodes))\n", + "p.set_val('vars.v', 10 * np.random.random(num_nodes))\n", + "p.set_val('vars.theta', 10 * np.random.uniform(1, 179, num_nodes))\n", "\n", "p.run_model()\n", "cpd = p.check_partials(method='cs', compact_print=True)" @@ -318,7 +312,7 @@ } }, "kernelspec": { - "display_name": "Python 3", + "display_name": "dev-local", "language": "python", "name": "python3" }, @@ -332,7 +326,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.3" + "version": "3.13.13" } }, "nbformat": 4, diff --git a/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_controls.ipynb b/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_controls.ipynb index a1e6870eaf..3de9816abf 100644 --- a/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_controls.ipynb +++ b/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_controls.ipynb @@ -73,7 +73,7 @@ }, "outputs": [], "source": [ - "om.display_source(\"dymos.examples.brachistochrone.doc.brachistochrone_ode\")" + "om.display_source(\"dymos.examples.brachistochrone.brachistochrone_ode\")" ] }, { diff --git a/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_states.ipynb b/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_states.ipynb index faf62aab77..ca7813ffe2 100644 --- a/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_states.ipynb +++ b/docs/dymos_book/examples/brachistochrone/brachistochrone_upstream_states.ipynb @@ -72,23 +72,23 @@ }, "outputs": [], "source": [ - "om.display_source(\"dymos.examples.brachistochrone.doc.brachistochrone_ode\")" + "om.display_source(\"dymos.examples.brachistochrone.brachistochrone_ode\")" ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "active-ipynb", - "remove-input", - "remove-output" - ] - }, - "outputs": [], - "source": [ - "%matplotlib inline" + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "active-ipynb", + "remove-input", + "remove-output" ] + }, + "outputs": [], + "source": [ + "%matplotlib inline" + ] }, { "cell_type": "code", @@ -272,7 +272,7 @@ } }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "dev-local", "language": "python", "name": "python3" }, @@ -286,7 +286,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.0" + "version": "3.13.13" }, "orphan": true }, diff --git a/docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb b/docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb index 4f22396eed..e7624e3751 100644 --- a/docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb +++ b/docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb @@ -101,7 +101,7 @@ }, "outputs": [], "source": [ - "om.display_source(\"dymos.examples.min_time_climb.doc.aero_partial_coloring\")" + "om.display_source(\"dymos.examples.min_time_climb.test.aero_partial_coloring\")" ] }, { @@ -114,7 +114,7 @@ "from dymos.models.eom import FlightPathEOM2D\n", "from dymos.examples.min_time_climb.prop import PropGroup\n", "from dymos.models.atmosphere import USatm1976Comp\n", - "from dymos.examples.min_time_climb.doc.aero_partial_coloring import AeroGroup\n", + "from dymos.examples.min_time_climb.test.aero_partial_coloring import AeroGroup\n", "\n", "\n", "class MinTimeClimbODE(om.Group):\n", @@ -348,7 +348,7 @@ } }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "dev-local", "language": "python", "name": "python3" }, @@ -362,7 +362,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.0" + "version": "3.13.13" } }, "nbformat": 4, diff --git a/docs/dymos_book/faq/tandem_phases.ipynb b/docs/dymos_book/faq/tandem_phases.ipynb index 17e7ba1570..f267de5736 100644 --- a/docs/dymos_book/faq/tandem_phases.ipynb +++ b/docs/dymos_book/faq/tandem_phases.ipynb @@ -86,7 +86,7 @@ }, "outputs": [], "source": [ - "om.display_source('dymos.examples.brachistochrone.doc.test_doc_brachistochrone_tandem_phases.BrachistochroneArclengthODE')" + "om.display_source('dymos.examples.brachistochrone.test.test_tandem_phases.BrachistochroneArclengthODE')" ] }, { @@ -138,7 +138,7 @@ " def compute(self, inputs, outputs):\n", " theta = inputs['theta']\n", " v = inputs['v']\n", - " outputs['Sdot'] = np.sqrt(1.0 + (1.0/np.tan(theta))**2) * v * np.sin(theta)\n", + " outputs['Sdot'] = np.sqrt(1.0 + (1.0 / np.tan(theta))**2) * v * np.sin(theta)\n", "\n", " def compute_partials(self, inputs, jacobian):\n", " theta = inputs['theta']\n", @@ -303,7 +303,7 @@ "source": [ "from openmdao.utils.assert_utils import assert_near_equal\n", "\n", - "expected = np.sqrt((10-0)**2 + (10 - 5)**2)\n", + "expected = np.sqrt((10 - 0) ** 2 + (10 - 5)**2)\n", "assert_near_equal(p.get_val('phase1.timeseries.S')[-1], expected, tolerance=1.0E-3)" ] } @@ -319,7 +319,7 @@ } }, "kernelspec": { - "display_name": "py312", + "display_name": "dev-local", "language": "python", "name": "python3" }, @@ -333,7 +333,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.8" + "version": "3.13.13" } }, "nbformat": 4, diff --git a/docs/dymos_book/faq/use_partial_coloring.ipynb b/docs/dymos_book/faq/use_partial_coloring.ipynb index 2860192b6a..d3eb1c904a 100644 --- a/docs/dymos_book/faq/use_partial_coloring.ipynb +++ b/docs/dymos_book/faq/use_partial_coloring.ipynb @@ -71,7 +71,7 @@ }, "outputs": [], "source": [ - "om.display_source(\"dymos.examples.min_time_climb.doc.dynamic_pressure_comp_partial_coloring.DynamicPressureCompFD\")" + "om.display_source(\"dymos.examples.min_time_climb.test.dynamic_pressure_comp_partial_coloring.DynamicPressureCompFD\")" ] }, { @@ -211,7 +211,7 @@ }, "outputs": [], "source": [ - "om.display_source(\"dymos.examples.min_time_climb.doc.min_time_climb_ode_partial_coloring.MinTimeClimbODE\")" + "om.display_source(\"dymos.examples.min_time_climb.test.min_time_climb_ode_partial_coloring.MinTimeClimbODE\")" ] }, { @@ -242,7 +242,7 @@ "import openmdao.api as om\n", "\n", "import dymos as dm\n", - "from dymos.examples.min_time_climb.doc.min_time_climb_ode_partial_coloring import MinTimeClimbODE\n", + "from dymos.examples.min_time_climb.test.min_time_climb_ode_partial_coloring import MinTimeClimbODE\n", "\n", "\n", "for fd in (False, True):\n", @@ -441,7 +441,7 @@ } }, "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "dev-local", "language": "python", "name": "python3" }, @@ -455,7 +455,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.13.13" } }, "nbformat": 4, diff --git a/docs/dymos_book/getting_started/scripts/fitter_animation_1_seg.py b/docs/dymos_book/getting_started/scripts/fitter_animation_1_seg.py index 0098dd6776..7c5c074a92 100644 --- a/docs/dymos_book/getting_started/scripts/fitter_animation_1_seg.py +++ b/docs/dymos_book/getting_started/scripts/fitter_animation_1_seg.py @@ -8,7 +8,7 @@ from matplotlib.animation import FuncAnimation # plt.switch_backend('Agg') # disable plotting to the screen -from dymos.examples.oscillator.doc.oscillator_ode import OscillatorODE +from dymos.examples.oscillator.oscillator_ode import OscillatorODE NUM_SEG = 1 ORDER = 5 diff --git a/dymos/examples/min_time_climb/doc/__init__.py b/dymos/examples/min_time_climb/doc/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/dymos/examples/min_time_climb/doc/aero_partial_coloring.py b/dymos/examples/min_time_climb/test/aero_partial_coloring.py similarity index 100% rename from dymos/examples/min_time_climb/doc/aero_partial_coloring.py rename to dymos/examples/min_time_climb/test/aero_partial_coloring.py diff --git a/dymos/examples/min_time_climb/doc/dynamic_pressure_comp_partial_coloring.py b/dymos/examples/min_time_climb/test/dynamic_pressure_comp_partial_coloring.py similarity index 100% rename from dymos/examples/min_time_climb/doc/dynamic_pressure_comp_partial_coloring.py rename to dymos/examples/min_time_climb/test/dynamic_pressure_comp_partial_coloring.py diff --git a/dymos/examples/min_time_climb/doc/min_time_climb_ode_partial_coloring.py b/dymos/examples/min_time_climb/test/min_time_climb_ode_partial_coloring.py similarity index 100% rename from dymos/examples/min_time_climb/doc/min_time_climb_ode_partial_coloring.py rename to dymos/examples/min_time_climb/test/min_time_climb_ode_partial_coloring.py diff --git a/dymos/examples/min_time_climb/test/test_doc_min_time_climb_partial_coloring.py b/dymos/examples/min_time_climb/test/test_doc_min_time_climb_partial_coloring.py index bc2684e64f..6ef41ec0d4 100644 --- a/dymos/examples/min_time_climb/test/test_doc_min_time_climb_partial_coloring.py +++ b/dymos/examples/min_time_climb/test/test_doc_min_time_climb_partial_coloring.py @@ -21,7 +21,7 @@ def test_min_time_climb_partial_coloring(self): from openmdao.utils.assert_utils import assert_near_equal import dymos as dm - from dymos.examples.min_time_climb.doc.min_time_climb_ode_partial_coloring import MinTimeClimbODE + from dymos.examples.min_time_climb.test.min_time_climb_ode_partial_coloring import MinTimeClimbODE # # Instantiate the problem and configure the optimization driver diff --git a/dymos/transcriptions/test/test_set_vals.py b/dymos/transcriptions/test/test_set_vals.py index a09ec7b419..009e80970f 100644 --- a/dymos/transcriptions/test/test_set_vals.py +++ b/dymos/transcriptions/test/test_set_vals.py @@ -8,7 +8,7 @@ import numpy as np import openmdao.api as om import dymos as dm -from dymos.examples.brachistochrone.doc.brachistochrone_ode import BrachistochroneODE +from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE @use_tempdirs