Dymos/OpenMDAO phase linkage solution reading error

Question

I am running and investigating the Dymos/OpenMDAO example provided here: SSTO Earth Launch

It runs fine, with different boundary conditions and parameters.

My goal is to make the trajectory in two phases, so that I can simulate a two stage rocket, where the first stage will detach and produce a mass "jump".

I am building an step by step model. So the first is to set up two phases, same boundary conditions (initial and final) and link them using:

traj.link_phases(phases=['phase0', 'phase1'], vars=['*'], connected=True)

It works as expected, same trajectory, but with two phases.

I am preparing the model to drop the weight of the burnt first stage. But now only changing the linkages of all variables to be explicit (control variables excluded):

traj.link_phases(phases=['phase0', 'phase1'], vars=['time','x','y','vx','vy','m'])

The case works, but when loading the solution for plotting i got the error:

Exception has occurred: KeyError
'traj.linkages.phase0:m_final|phase1:m_initial'
KeyError: 'traj.linkages.phase0:m_final|phase1:m_initial'

During handling of the above exception, another exception occurred:

  File "/dymos/ssto/main_phases_works.py", line 113, in <module>
    sol = om.CaseReader('dymos_solution.db').get_case('final')
KeyError: 'traj.linkages.phase0:m_final|phase1:m_initial'

Apparently, it is not reading the linkages information for plotting? Anyone know why?

Thank you jorge

This may be a bug. Could you please post an issue at github.com/OpenMDAO/dymos with the sample code if possible so that we can resolve it? — Rob Falck, Commented Jun 25 at 12:11
Sure I will, let the grab the original files and post them there — jordi nosaurio, Commented Jun 26 at 13:03
I know what is going on, the right comparison is with: traj.link_phases(phases=['phase0', 'phase1'], vars=['time','x','y','vx','vy','m'], connected=True) And works, with connected = True. What I do not know is why it does not work without it, since it should create constrains for the optimizer. — jordi nosaurio, Commented Jun 26 at 19:46
It's likely just a logic error when loading the case for plotting. — Rob Falck, Commented Jun 26 at 20:35
I cannot get it converged when using connected = False. even for just the mass. I am going to post the code — jordi nosaurio, Commented Jun 27 at 21:21

jordi nosaurio · Accepted Answer · 2024-06-27 22:01:06Z

The modified SSTO code is for just one phase (just for proper comparison)

import matplotlib.pyplot as plt
import openmdao.api as om
import dymos as dm
from dymos.examples.plotting import plot_results

p = om.Problem(model=om.Group())
p.driver = om.pyOptSparseDriver()

p.driver.options['optimizer'] = 'IPOPT'

p.driver.opt_settings['derivative_test'] = 'first-order'
p.driver.opt_settings['mu_strategy'] = 'adaptive'
p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
p.driver.opt_settings['mu_init'] = 0.01
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
p.driver.opt_settings['max_iter']= 1000
p.driver.opt_settings['tol']= 1e-2
p.driver.declare_coloring()

from launch_vehicle_ode import LaunchVehicleODE

traj = dm.Trajectory()

phase0 = dm.Phase(ode_class=LaunchVehicleODE, ode_init_kwargs={'CD': 1.0, 'S': 0.2},
                 transcription=dm.Radau(num_segments=10, order=3, compressed=False))

traj.add_phase('phase0', phase0)

p.model.add_subsystem('traj', traj)
#
phase0.set_time_options(    fix_initial=True, duration_bounds=(1, 500))
phase0.add_state('x',       fix_initial=True, rate_source='xdot', fix_final=True, lower = 0.0)
phase0.add_state('y',       fix_initial=True,rate_source='ydot',  fix_final=True, lower = 0.0)
phase0.add_state('vx',      fix_initial=True,rate_source='vxdot', fix_final=True, lower = 0.0)
phase0.add_state('vy',      fix_initial=True, rate_source='vydot',fix_final=True, lower = 0.0)
phase0.add_state('m',       fix_final=True,rate_source='mdot', lower = 500.0)
phase0.add_control('theta', fix_initial=True, units='deg', lower=-90.0, upper=90.0, targets=['theta'])
phase0.add_parameter('thrust', val = 25000.0, opt=False, targets=['thrust'])
#
phase0.add_objective('time', loc='final')

p.model.linear_solver = om.DirectSolver()
#
p.setup(check=True)

p.set_val('traj.phase0.t_initial', 0.0)
p.set_val('traj.phase0.t_duration', 150.0)
p.set_val('traj.phase0.states:x', phase0.interp('x', [0, 80.0E3]))
p.set_val('traj.phase0.states:y', phase0.interp('y', [0, 50.0E3]))
p.set_val('traj.phase0.states:vx', phase0.interp('vx', [0.0, 2000.0]))
p.set_val('traj.phase0.states:vy', phase0.interp('vy', [10.0, 0.0]))
p.set_val('traj.phase0.states:m', phase0.interp('m', [1500.0, 500.0]))
p.set_val('traj.phase0.controls:theta', phase0.interp('theta', [89.0, 0.0]))
p.set_val('traj.phase0.parameters:thrust', 25000, units='N')
#
phase0.set_refine_options(tol=1.0E-2)

dm.run_problem(p, simulate=True, run_driver=True, refine_iteration_limit=1)

#POST
sol = om.CaseReader('dymos_solution.db').get_case('final')
sim = om.CaseReader('dymos_simulation.db').get_case('final')

plot_results([('traj.phase0.timeseries.x', 'traj.phase0.timeseries.y',
               'range (m)', 'altitude (m)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.y',
               'time (s)', 'altitude (m)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.vx',
               'time (s)', 'xspeed (m/s)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.vy',
               'time (s)', 'yspeed (m/s)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.theta',
               'time (s)', 'theta (deg)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.m',
               'time (s)', 'm (kg)'),
               ],
             title='Supersonic Minimum Time-to-Climb Solution',
             p_sol=sol, p_sim=sim)

plt.show()

The code for two phases is just a split it two and link them:

import matplotlib.pyplot as plt
import openmdao.api as om
import dymos as dm
from dymos.examples.plotting import plot_results

p = om.Problem(model=om.Group())
p.driver = om.pyOptSparseDriver()

p.driver.options['optimizer'] = 'IPOPT'

p.driver.opt_settings['derivative_test'] = 'first-order'
p.driver.opt_settings['mu_strategy'] = 'adaptive'
p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
p.driver.opt_settings['mu_init'] = 0.01
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
p.driver.opt_settings['max_iter']= 1000
p.driver.opt_settings['tol']= 1e-2
p.driver.declare_coloring()

from launch_vehicle_ode import LaunchVehicleODE

traj = dm.Trajectory()

phase0 = dm.Phase(ode_class=LaunchVehicleODE, ode_init_kwargs={'CD': 1.0, 'S': 0.2},
                 transcription=dm.Radau(num_segments=10, order=3, compressed=False))
phase1 = dm.Phase(ode_class=LaunchVehicleODE, ode_init_kwargs={'CD': 1.0, 'S': 0.2},
                 transcription=dm.Radau(num_segments=10, order=3, compressed=False))

traj.add_phase('phase0', phase0)
traj.add_phase('phase1', phase1)

p.model.add_subsystem('traj', traj)

phase0.set_time_options(fix_initial=True, duration_bounds=(1, 500))
phase0.add_state('x',  fix_initial=True, rate_source='xdot', lower = 0.0)
phase0.add_state('y',  fix_initial=True,rate_source='ydot', lower = 0.0)
phase0.add_state('vx', fix_initial=True,rate_source='vxdot', lower = 0.0)
phase0.add_state('vy', fix_initial=True, rate_source='vydot', lower = 0.0)
phase0.add_state('m',  rate_source='mdot', lower=1100.0)
phase0.add_control('theta', fix_initial=True, units='deg', lower=-90.0, upper=90.0, targets=['theta'])
phase0.add_parameter('thrust', val = 25000.0, opt=False, targets=['thrust'])

phase1.set_time_options(duration_bounds=(1, 500))
phase1.add_state('x', rate_source='xdot',  fix_final=True, lower = 0.0)
phase1.add_state('y', rate_source='ydot',  fix_final=True, lower = 0.0)
phase1.add_state('vx',rate_source='vxdot', fix_final=True, lower = 0.0)
phase1.add_state('vy',rate_source='vydot', fix_final=True, lower = 0.0)
phase1.add_state('m', rate_source='mdot',   fix_final=True, lower = 500.0)
phase1.add_control('theta', units='deg', lower=-90.0, upper=90.0, targets=['theta'])
phase1.add_parameter('thrust', val = 25000.0, opt=False, targets=['thrust'])
#
# Set the options for our constraints and objective
#

phase1.add_objective('time', loc='final')

p.model.linear_solver = om.DirectSolver()

p.setup(check=True)

p.set_val('traj.phase0.t_initial', 0.0)
p.set_val('traj.phase0.t_duration', 150.0)
p.set_val('traj.phase0.states:x', phase0.interp('x',   [0, 50.0E3]))
p.set_val('traj.phase0.states:y', phase0.interp('y',   [0, 45.0E3]))
p.set_val('traj.phase0.states:vx', phase0.interp('vx', [0.0, 1000.0]))
p.set_val('traj.phase0.states:vy', phase0.interp('vy', [10.0, 0.0]))
p.set_val('traj.phase0.states:m', phase0.interp('m',   [1500.0, 1100.0]))
p.set_val('traj.phase0.controls:theta', phase0.interp('theta', [89.0, 0.0]))
p.set_val('traj.phase0.parameters:thrust', 25000, units='N')

p.set_val('traj.phase1.t_initial', 100.0)
p.set_val('traj.phase1.t_duration', 150.0)
p.set_val('traj.phase1.states:x', phase1.interp('x',   [50.0E3, 80.0E3]))
p.set_val('traj.phase1.states:y', phase1.interp('y',   [45.0E3, 50.0E3]))
p.set_val('traj.phase1.states:vx', phase1.interp('vx', [1000.0, 2000.0]))
p.set_val('traj.phase1.states:vy', phase1.interp('vy', [0.0, 0.0]))
p.set_val('traj.phase1.states:m', phase1.interp('m',   [1100.0, 500.0]))
p.set_val('traj.phase1.controls:theta', phase1.interp('theta', [0.0, 0.0]))
p.set_val('traj.phase1.parameters:thrust', 25000, units='N')

traj.link_phases(phases=['phase0', 'phase1'], vars=['*'], connected=True)

phase0.set_refine_options(tol=1.0E-2)
phase1.set_refine_options(tol=1.0E-2)
dm.run_problem(p, simulate=False, run_driver=True, refine_iteration_limit=1, 
               simulate_kwargs={'times_per_seg' : 100, 'rtol': 1.0E-2})

p.driver.scaling_report()

sol = om.CaseReader('dymos_solution.db').get_case('final')

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))

time_imp = {'phase0': p.get_val('traj.phase0.timeseries.time'),
            'phase1': p.get_val('traj.phase1.timeseries.time')}
r_imp = {'phase0': p.get_val('traj.phase0.timeseries.x'),
         'phase1': p.get_val('traj.phase1.timeseries.x')}
h_imp = {'phase0': p.get_val('traj.phase0.timeseries.y'),
         'phase1': p.get_val('traj.phase1.timeseries.y')}

axes.plot(r_imp['phase0'], h_imp['phase0'], 'bo')
axes.plot(r_imp['phase1'], h_imp['phase1'], 'ro')

axes.set_xlabel('range (m)')
axes.set_ylabel('altitude (m)')
axes.grid(True)

fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(10, 6))
states = ['x', 'y', 'vx', 'vy', 'm','theta']
for i, state in enumerate(states):
    x_imp = {'phase0': sol.get_val(f'traj.phase0.timeseries.{state}'),
             'phase1': sol.get_val(f'traj.phase1.timeseries.{state}')}
    axes[i].set_ylabel(state)
    axes[i].grid(True)
    axes[i].plot(time_imp['phase0'], x_imp['phase0'], 'bo')
    axes[i].plot(time_imp['phase1'], x_imp['phase1'], 'ro')

plt.show()

Linking the phases with connected=True is ok and I got convergence. Nevertheless my purpose is a "jump" in mass to simulate the drop of the burned first stage in this two stage rocket. So I am trying to link them as (also changing the mass initial conditions):

traj.link_phases(phases=['phase0', 'phase1'], vars=['time','x','y','vx','vy'], connected=True)
traj.add_linkage_constraint('phase0', 'phase1', 'm', 'm', loc_a='final', loc_b='initial',
                            mult_a=1.0, mult_b=-1.0, equals = 200.0, units='kg', connected=True)

Although it seems to interpret the right mass "jump" I cannot get it converged. I would like to remove the "connected" setting in the mass (it is now fully continuous as expected) but either I got no convergence and/or negative times in the second phase, or the plotting error appears.

For example I cannot converge this try:

import matplotlib.pyplot as plt
import openmdao.api as om
import dymos as dm
from dymos.examples.plotting import plot_results

p = om.Problem(model=om.Group())
p.driver = om.pyOptSparseDriver()

p.driver.options['optimizer'] = 'IPOPT'

p.driver.opt_settings['derivative_test'] = 'first-order'
p.driver.opt_settings['mu_strategy'] = 'adaptive'
p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
p.driver.opt_settings['mu_init'] = 0.01
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
p.driver.opt_settings['max_iter']= 1000
p.driver.declare_coloring()

from launch_vehicle_ode import LaunchVehicleODE

traj = dm.Trajectory()

phase0 = dm.Phase(ode_class=LaunchVehicleODE, ode_init_kwargs={'CD': 1.0, 'S': 0.2},
                 transcription=dm.Radau(num_segments=10, order=3, compressed=False))
phase1 = dm.Phase(ode_class=LaunchVehicleODE, ode_init_kwargs={'CD': 1.0, 'S': 0.2},
                 transcription=dm.Radau(num_segments=10, order=3, compressed=False))

traj.add_phase('phase0', phase0)
traj.add_phase('phase1', phase1)

p.model.add_subsystem('traj', traj)

phase0.set_time_options(fix_initial=True, duration_bounds=(1, 500))
phase0.add_state('x',  fix_initial=True, rate_source='xdot', lower = 0.0, ref=1.0E5, defect_ref=1.0E4)
phase0.add_state('y',  fix_initial=True,rate_source='ydot', lower = 0.0, ref=1.0E5, defect_ref=1.0E4)
phase0.add_state('vx', fix_initial=True,rate_source='vxdot', lower = 0.0, ref=1.0E3, defect_ref=1.0E3)
phase0.add_state('vy', fix_initial=True, rate_source='vydot', lower = 0.0, ref=1.0E3, defect_ref=1.0E3)
phase0.add_state('m',  rate_source='mdot', lower=1100.0, ref=1.0E3, defect_ref=1.0E2)
phase0.add_control('theta', fix_initial=True, units='deg', lower=-90.0, upper=90.0, targets=['theta'])
phase0.add_parameter('thrust', val = 25000.0, opt=False, targets=['thrust'])

phase1.set_time_options(duration_bounds=(1, 500))
phase1.add_state('x', rate_source='xdot',  fix_final=True, lower = 0.0, ref=1.0E5, defect_ref=1.0E4)
phase1.add_state('y', rate_source='ydot',  fix_final=True, lower = 0.0, ref=1.0E5, defect_ref=1.0E4)
phase1.add_state('vx',rate_source='vxdot', fix_final=True, lower = 0.0, ref=1.0E3, defect_ref=1.0E3)
phase1.add_state('vy',rate_source='vydot', fix_final=True, lower = 0.0, ref=1.0E3, defect_ref=1.0E3)
phase1.add_state('m', rate_source='mdot',   fix_final=True, upper = 900.0, ref=1.0E3, defect_ref=1.0E2)
phase1.add_control('theta', units='deg', lower=-90.0, upper=90.0, targets=['theta'])
phase1.add_parameter('thrust', val = 25000.0, opt=False, targets=['thrust'])

phase1.add_objective('time', loc='final')

p.model.linear_solver = om.DirectSolver()

p.setup(check=True)

phase0.set_time_val(initial=0.0, duration=150.0)
phase0.set_state_val('x',  [0, 50.0E3])
phase0.set_state_val('y',  [0, 45.0E3])
phase0.set_state_val('vx', [0.0, 1000.0])
phase0.set_state_val('vy', [10.0, 0.0])
phase0.set_state_val('m',  [1500.0, 1100.0])
phase0.set_control_val('theta', [89.0, 0.0],units='deg')
phase0.set_parameter_val('thrust', 25000.0, units='N')

phase1.set_time_val(initial=0.0, duration=150.0)
phase1.set_state_val('x',  [50.0E3, 80.0E3])
phase1.set_state_val('y',  [45.0E3, 50.0E3])
phase1.set_state_val('vx', [1000.0, 2000.0])
phase1.set_state_val('vy', [0.0, 0.0])
phase1.set_state_val('m',  [900.0, 400.0])
phase1.set_control_val('theta', [0.0, 0.0],units='deg')
phase1.set_parameter_val('thrust', 25000.0, units='N')

traj.link_phases(phases=['phase0', 'phase1'], vars=['time','x','y','vx','vy'], connected=True)
traj.add_linkage_constraint('phase0', 'phase1', 'm', 'm', loc_a='final', loc_b='initial',
                            mult_a=1.0, mult_b=-1.0, equals = 200.0, units='kg', ref=1.0E3, connected=True)

dm.run_problem(p, simulate=False, run_driver=True, refine_iteration_limit=5,
               simulate_kwargs={'times_per_seg' : 100, 'rtol': 1.0E-3})

p.driver.scaling_report()

sol = om.CaseReader('dymos_solution.db').get_case('final')

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))

time_imp = {'phase0': p.get_val('traj.phase0.timeseries.time'),
            'phase1': p.get_val('traj.phase1.timeseries.time')}
r_imp = {'phase0': p.get_val('traj.phase0.timeseries.x'),
         'phase1': p.get_val('traj.phase1.timeseries.x')}
h_imp = {'phase0': p.get_val('traj.phase0.timeseries.y'),
         'phase1': p.get_val('traj.phase1.timeseries.y')}

axes.plot(r_imp['phase0'], h_imp['phase0'], 'bo')
axes.plot(r_imp['phase1'], h_imp['phase1'], 'ro')

axes.set_xlabel('range (m)')
axes.set_ylabel('altitude (m)')
axes.grid(True)

fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(10, 6))
states = ['x', 'y', 'vx', 'vy', 'm','theta']
for i, state in enumerate(states):
    x_imp = {'phase0': sol.get_val(f'traj.phase0.timeseries.{state}'),
             'phase1': sol.get_val(f'traj.phase1.timeseries.{state}')}
    axes[i].set_ylabel(state)
    axes[i].grid(True)

    axes[i].plot(time_imp['phase0'], x_imp['phase0'], 'bo')
    axes[i].plot(time_imp['phase1'], x_imp['phase1'], 'ro')

plt.show()

I am puzzled about why appears to be such a difficult problem. Any clues?

thank you

Collectives™ on Stack Overflow

Dymos/OpenMDAO phase linkage solution reading error

1 Answer 1

Not the answer you're looking for? Browse other questions tagged
optimization
controls
openmdao
or ask your own question.

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Not the answer you're looking for? Browse other questions tagged optimizationcontrolsopenmdao or ask your own question.

Related

Not the answer you're looking for? Browse other questions tagged
optimization
controls
openmdao
or ask your own question.