Define an AnalyticPhase for which the solution to the ODE is known analytically. Rather than taking an ODE, take a system that provides a solution for any time.

This would be useful to have on some of the project's we're currently working, so I'm bumping the priority up.

For the initial cut, we need a new Phase called AnalyticPhase that inherits fromPhase and overrides some of its methods.

Handling states

In an analytic phase, states are really just outputs of the form f(x0, t, p).
That is analytic phase will, at least initially, fundamentally just solve the IVP analytically, providing states just as any other output.
We'll maintain an official list of states for the sake of keeping things consistent with other phases in dymos when it comes to linking them together and stuff like that.

We will override the set_state_options method so that it doesn't have a rate_src, but instead has val_src - the output in the ODE that provides the value of the state at a given time.

Handling controls

Initially there will be no support for controls in AnalyticPhase.
Attempting to add a control or polynomial control will raise an exception.
I'm not certain that there's any example where an ODE with dynamic controls has an analytic solution.

TImeseries outputs

Timeseries outputs will still be provided in the form phasename.timeseries.* in order to be consistent with other phases, to support constraints and linkages.

Transcriptions

I don't think it makes sense that AnalyticPhases have transcriptions. In this case we're not really converting an infinite-dimensional optimal control problem to a discrete optimization problem.

Still, the mechanisms in dymos dictate that the transcription helps to setup/configure the phase, so we'll have an AnalyticTranscription that will just automatically be associated with the AnalyticPhase.

The next question is, where do the nodes go?
For this initial cut I think it makes sense that AnalyticPhase has a num_nodes option and provides output nodes by using the formula np.linspace(t_initial, t_initial + t_duration, num_nodes).

Constraints

We should still support path and boundary constraints as we do with the other Phase/transcription types.
From a control viewpoint, this phase is similar to ExplicitShooting, in that it fundamentally is governed by values of the initial states, parameters, initial time and duration.

Possible example cases

1. Propagating a 2-body orbit using f and g series. Example here

2. Propagate some simple ODE, such as the one here

SSTO Linear Tangent docs need to be updated - currently don't contain the correct code.
Multi-phase cannonball contains a math display error that needs to be corrected.

Trajectory objects are Groups which serve to organize a set of one or more phases.

Each Trajectory group consists of two primary subsystems:

• phases
  A ParallelGroup within the Trajectory that can hold each phase, and potentially evaluate them in parallel
• linkages
  A PhaseLinkageComp that manages all of the linkage constraints between phases

Trajectory objects as currently envisioned have the following primary methods:

• add_phase(phase_name, phase_object, **kwargs)

  Phases are added to a trajectory with an add_phase method. Functionally, it's very similar to the add_subsystem method of ParallelGroup. Phases are cached in an internal dictionary _phases, and then actually added to the phases subgroup during setup.

• link_phase(phases, vars=['*'], locs=('++', '--'))
  Trajectory objects provide a link_phases method that provides an ordering of phases in the trajectory. The goal of link_phases is to make common phase linkages (sequential) trivial, and complicated phase linkages (parallel, multiple branches, etc) tractable.

  Argument phases provides a series of two or more phases that follow one another in sequence. Argument vars provides the name of the time, state, controls, or control rates to be linked. Providing the value '*' in the list here will imply that time and all common state variables between the phases are to be linked. Providing a tuple-pair as one of the elements in vars allows variables of different names to be connected across phases.

  Finally locs provides the location of the connection for each pair of phases in the given sequence, using the --, -+, +-, ++ convention for states and controls.

  Ultimately link_phases should probably provide the same functionality as OTIS's global constraints, where any variables can be connected between two phases.

• simulate(times, parallel=True)
  Simulate each phase in the trajectory, returning the results of each simulation in a dictionary. Use python's multiprocessing module to simulate each phase simultaneously when possible, to improve performance.

• get_values(var, phases=None, nodes='all', flat=False)
  Return the values of a given variable at the given nodes through the given sequence of phases in the trajectory. Omitting phases returns the variable in all phases. If not present, return nan for each node in the phase. Might have to think about how we treat variables which may have different shapes in different phases, since plotting a 3-vector from various phases, and then just a 1-vector of nans from a phase where the variable is not present. If flat, rather than returning a dictionary, return a continuous array across phases.

For instance, in the double integrator example, setting opt_duration=False causes the following.

NameError: Output 't_duration' does not exist for connection in 'phase0' from 't_duration' to 'continuity_comp.t_duration'.

RuntimeError: Output units of 'kn/s' for 'climb.control_rate_comp.control_rates:TAS_rate' are incompatible with input units of 'kn/s/s' for 'climb.path_constraints.all_values:TAS_rate'.

slightly more convenient then wrapping it in a list/tuple if getting only one phase.

Move direction out of phase options and put on state options as time_direction.
Add 'solve_subsystems': True to the default continuity solver options.
Add a test case that uses phase (optimizer constrained) linkages.

set_values is the complement to get_values. Given a continuous time-series of data, interpolate as necessary and insert into the appropriate state/control values.

set_values should remove the need for the user directly calling phase.interpolate, and therefore has to support a similar option set (xs, ys, interpolation kind). Signature should be

phase.set_values(name, values, units=None, interp=None, interp_xs=None)

name is the dymos variable name, where internally we check to see if it is a time, state, control, or ODE variable.

values are the values to be assigned.

units are the units in which we are assigning the variable. Internally the appropriate unit conversion will occur.

interp is the type of interpolation to perform on the given values. Default should be None, which will raise an error if the given values aren't of the correct length.

interp_xs are the time/independent variable values to use, equivalent to the xs argument to phase.interpolate. If interp is None or 'linear' then xs are not necessary, otherwise an error should be raised in its absence.

Sometimes it is necessary for ODE to know the current elapsed time since a phase started (i.e. an event).

1. Add a new argument for time_phase targets to set_time_options
2. Add time_phase output to time_comp to be connected to any time_phase targets given
3. Allow time_phase to be used in the same contexts as time (boundary constraints, objectives, etc).

Test that ODE_class is a class and that it has a class property called "ode_options" that is an instance of dymos.ODEOptions

Compressed transcription seems to be giving a substantial performance advantage over uncompressed. It might even be worth looking into testing if pyoptsparse is ordering linear constraints in the SNOPT-recommended way.

import numpy as np

from openmdao.api import ExplicitComponent, Group, Problem
from dymos import declare_time, declare_state, declare_parameter, Phase, ODEOptions


n_traj = 4


class MyComp(ExplicitComponent):

    def initialize(self):
        self.options.declare('n_traj', types=int)
        self.options.declare('num_nodes', types=int)

    def setup(self):
        nn = self.options['num_nodes']
        n_traj = self.options['n_traj']
        self.add_input('time', val=np.zeros(nn))
        self.add_input('alpha', shape=np.zeros((n_traj, 2)).shape)

        self.add_output('y', val=np.zeros(nn))

    def compute(self, inputs, outputs):
        pass

    def compute_partials(self, inputs, partials):
        pass

class MyODE(Group):
    ode_options = ODEOptions()
    ode_options.declare_time(units='s', targets = ['comp.time'])
    ode_options.declare_state(name='F', rate_source='comp.y')
    ode_options.declare_parameter(name='alpha', shape=(4,2), targets = 'comp.alpha')

    def initialize(self):   
        self.options.declare('num_nodes', types=int)
        self.options.declare('n_traj', default=2, types=int)

    def setup(self):
        nn = self.options['num_nodes']
        n_traj = self.options['n_traj']

        self.add_subsystem(name='comp',
                           subsys=MyComp(num_nodes=nn,
                                         n_traj=n_traj))

p = Problem(model=Group())

phase = Phase(transcription='gauss-lobatto',
              ode_class=MyODE,
              num_segments=25,
              transcription_order=3,
              compressed=True)

p.model.add_subsystem('phase0', phase)

phase.add_input_parameter('alpha', val=np.ones((n_traj, 2)), units='m')

p.setup()

Results in error:

ValueError: The source indices [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] do not specify a valid shape for the connection 'phase0.input_params.input_parameters:alpha_out' to 'phase0.rhs_disc.comp.alpha'. The target shape is (2, 2) but indices are (50,).

Sometimes it's convenient to define the rate of a variable as a control, and then integrate it to keep track of the control value. Currently, doing so requires a 'pass-through' component in the ODE that just echoes the control value as an output.

This change would allow the user to declare_state with a rate_source of None in the ODE, which could then be overridden with a rate_source option to set_state_options.

For some problems, like the double integrator, the problem could potentially be solved without an ODE at all.

This enhancement allows trajectories to own design parameters. The design parameters can then be connected to each phase within the trajectory.

Syntax:
trajectory.add_design_parameter(name, val, units=None, targets=None, opt=True, input=False, lower, upper, adder, scaler, ref0, ref)

All options are the same for the phase add_design_parameter method except for targets. If given, targets is a dictionary that maps each phase name to the ODE parameter name to which it is being connected. This allows a single design parameter at the trajectory level to be connected to variables which may have different names in different phases. If targets is None, assume that the given name is also the name of the target parameter in each phase.

Options fixed_initial and fixed_final are not sufficient to determine the direction of integration of each variable. Change the current state option direction (used only by RungeKuttaPhase at the moment) to time_direction.

Options connected_initial and connected_final specify that either the initial state value or final state value is to be connected to another output.

If options fixed_initial or fixed_final are True, then the optimizer owns those values and the initial/final values are not design variables. If they are fixed then they cannot be connected. fixed_final cannot be used in forward shooting (time_direction = 1) and option fixed_initial cannot be used with backward shooting time_direction = -1

If options fixed_initial or fixed_final are False, and the respective connect_initial or connect_final is False, then the initial/final value is governed by the optimizer as a design variable.

There are various combinations of these options that are invalid that need to be checked.

Favor the use of "design parameter" over static controls. Design parameters can exist in phases and in trajectory objects.

Unbounded initial or duration times have an upper bound of -INF_BOUND, which should be INF_BOUND.

it currently just gives a scalar t, which does not behave the same as the collocation phase.

Phases should provide a single interface for trajectory design and input parameters. This should be accomplished by using "promote as" to the appropriate variable name.

Dymos needs some indexing fixes to allow states, controls, and design/input parameters to be vectors. Also add a new test case that uses vectorized states and controls to verify that the implementation works.

Allow the user to specify phase-tau points at which the timeseries outputs should be provided. The timeseries output will then use interpolation to provide results at the desired points.

Use the Akima-spline based model that allows the sensitivity of atmospheric parameters w.r.t. altitude to be outputs.

*Restructure docs into the OpenMDAO setup,
*fix any Sphinx warnings,
*clean up pep 8 issues in docs files,
*add basic auto-documentation for dymos classes,
*fix auto-publishing to github-pages.

For example, if you do:
phase.set_time_options(fix_initial=True, fix_duration=True, initial=0.0, duration=2.5*3600.0)
the values for initial and duration are later replaced with their defaults, which is 1.0.

Additional API change: discuss with Rob the idea of changing initial and duration to initial_val and duration_val to be more consistent.

This change makes phase design parameters consistent with times in that they support external connections via the input=True argument to add_design_parameter.

As with time, once a design parameter is set to be input from an external source, the user is playing by big kid rules. Options on the design parameter dealing with optimization are invalid.

phase.simulate should take times as a number or an iterable.

• If a number, return that many evenly distributed points across the phase (np.linspace(t_initial, t_initial+t_duration, N).
• If an iterable, require it to be a monotonic increasing or decreasing sequence of times within the phase at which output is requested>
• It should be a warning to request an extrapolation of time

The commercial aircraft example is very nice, but I believe that the brachistochrone model was accidentally inserted instead of the aircraft model.

https://openmdao.github.io/dymos/examples/commercial_aircraft.html#the-ode-system-brachistochrone-ode-py

Encompasses all work performed leading up to the 1.0.0 major release.

In declare_parameter and declare_state, check that name is a valid openmdao variable name (no dots, etc). If it's invalid, raise a sensible ValueError.

WARNING: The following components have multiple inputs connected to the same source, which can introduce unnecessary data transfer overhead:
phase0.timeseries has inputs ['all_values:time', 'all_values:time_phase'] connected to phase0.time.time

This warning is present in the example from issue #138

Now that the post-processing get_values methods have been replaced by timeseries outputs, a trajectory should have a timeseries output so users can access timeseries across the entire trajectory.

This isn't completely trivial since a state/control variable may not exists in all phases across a trajectory, or more change from being a state in one phase to a control or design parameter in another.

Add custom control profiles at the phase level.

phase.add_control(..., control_class=X)

The control class will declare it's on design variables, promoted to 'controls:' in the phase.
The control class will provide control values at the control discretization nodes, which will then be connected through the standard control interpolation functions.

Add a new method (find_var(name, nodes)) to OpenMDAO Trajectories/Phases the path and indices of a variable/subset combination. This information is the same as that used by get_values to return data, but leaves accessing the actual data to the user.

Use the same time comp in GLM, Gauss-Lobatto, and Radau-PS phases

Part 1:

Phases get a new component timeseries_output_comp which aggregates all output variables.

By default, the following timeseries outputs are provided:

• time
• time_phase
• states
• controls
• design parameters
• input parameters
• constraints
• objective

User may add additional timeseries outputs using
phase.add_timeseries_output(ode_rel_path)

Outputs are accessible via path_to_phase.outputs:output_name

Part 2:

Remove the complexity of PhaseSimulationResults and TrajSimulationResults by having simulate only record the specified timeseries outputs.

When using add_input_parameter to pass static information from outside of Dymos into a phase, there is a shape mismatch error.

from openmdao.api import Group, ExecComp, Problem
from openmdao.api import DirectSolver, pyOptSparseDriver

from dymos import declare_time, declare_state, declare_parameter
from dymos import Phase
from dymos.utils.lgl import lgl
from dymos.models.eom import FlightPathEOM2D
         
import numpy as np
       
@declare_time(units='s')
@declare_state('v', rate_source='eom.v_dot', units='m/s')
@declare_state('h', rate_source='eom.h_dot', units='m')
@declare_parameter('m', targets='sum.m', units='kg', shape=(2, 2))

class TrajectoryODE(Group):
    def initialize(self):
        self.options.declare('num_nodes', types=int)    
    
    def setup(self):
        nn = self.options['num_nodes']
        
        self.add_subsystem('sum', ExecComp('m_tot = sum(m)', 
                                           m={'value':np.zeros(nn)}, 
                                           m_tot={'value':np.zeros(nn)}))
        
        self.add_subsystem('eom', FlightPathEOM2D(num_nodes=nn))
        
        self.connect('sum.m_tot','eom.m')
        
def BallDrop(optimizer='SLSQP', num_segments=1, transcription_order=5):
    
    p = Problem(model=Group())
    
    p.driver = pyOptSparseDriver()
    p.driver.options['optimizer'] = optimizer
    p.driver.options['dynamic_simul_derivs'] = True

    seg_ends, _ = lgl(num_segments + 1)

    phase = Phase(transcription='radau-ps',
                  ode_class=TrajectoryODE,
                  num_segments=num_segments, transcription_order=transcription_order,
                  segment_ends=seg_ends,
                  )    
    
    p.model.add_subsystem('phase0', phase)
    
    phase.set_time_options(initial_bounds=(0.0, 100.0), duration_bounds=(0.,100.))
    
    phase.set_state_options('h', fix_initial=True, fix_final=True, lower=0.0, units='m')
    phase.set_state_options('v', fix_initial=True, fix_final=False, units='m/s')    

    phase.add_input_parameter('m', val=[[1,1],[2,2]], units='kg')
    
    p.model.linear_solver = DirectSolver()
    
    p.setup(check=True, force_alloc_complex=True)

    p['phase0.t_initial'] = 0.0
    p['phase0.t_duration'] = 100.0

    p['phase0.states:h'] = phase.interpolate(ys=[20, 0], nodes='state_input')
    p['phase0.states:v'] = phase.interpolate(ys=[0, -5], nodes='state_input')

    
    p.run_model()

    print('Time:')
    print(p.get_val('phase0.timeseries.time'))
    print('')
    print('Height:')
    print(p.get_val('phase0.timeseries.states:h'))
    print('')
    print('Velocity:')
    print(p.get_val('phase0.timeseries.states:v'))

if __name__ == '__main__':
    p = BallDrop()

Traceback (most recent call last):
  File "mre_target_shape.py", line 84, in <module>
    p = BallDrop()
  File "mre_target_shape.py", line 62, in BallDrop
    p.setup(check=True, force_alloc_complex=True)
  File "/mnt/c/OMDAO/OpenMDAO.git/openmdao/core/problem.py", line 755, in setup
    derivatives)
  File "/mnt/c/OMDAO/OpenMDAO.git/openmdao/core/system.py", line 697, in _setup
    self._setup_connections(recurse=recurse)
  File "/mnt/c/OMDAO/OpenMDAO.git/openmdao/core/group.py", line 913, in _setup_connections
    subsys._setup_connections(recurse)
  File "/mnt/c/OMDAO/OpenMDAO.git/openmdao/core/group.py", line 1084, in _setup_connections
    len(out_shape), source_dimensions))
ValueError: The source indices [[[0 1]  [2 3]] [[0 1]  [2 3]] [[0 1]  [2 3]] [[0 1]  [2 3]] [[0 1]  [2 3]] [[0 1]  [2 3]]] do not specify a valid shape for the connection 'phase0.input_params.input_parameters:m_out' to 'phase0.rhs_all.sum.m'. The source has 3 dimensions but the indices expect 2.

Adding phase-wide interpolated controls will mean we have some controls whose input values are defined at some set of nodes in phase tau space. Adding this functionality will do the appropriate interpolation for the user.

Adds explicit phases which will initially use fixed-step RK4 integration across each segment of a phase.

By declaring an ODE parameter with dynamic = False, it can only be connected to design parameters or input parameters.

some ODE's

put GLM phase alongside gauss-lobatto and radau-ps phases in the file structure.
get rid of the ozone subfolder

This will make it easier for the user to find default values for arguments.
It will change Dymos such that calling set_state_options on the same state will override previous values.

• set_state_options
• set_time_options
• add_path_constraint
• add_boundary_constraint

When specifying phase linkages, add a new option to specify whether they are to be handled by a solver or by the optimizer:

traj = Trajectory(...)

...

traj.link_phases(phases=['a', 'b', 'c'], vars=['x', 'y', 'v'], connected=True|False, connect_start=True|False, connect_end=True|False)

I'm open to different option names in phase_linkages, like solver|optimizer?

Design parameters should be able to be used as objectives. Add a new test that extremizes a design parameter to verify operation.