gistfile1.txt

import numpy as np

import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal

import dymos as dm
from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE


class ParamComp(om.ExplicitComponent):

    def setup(self):

        self.add_input('aircraft:wing:sweep_distribution', np.array([0.3, 0.5, 0.9]), units='m')

    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
        pass


class ODEGroup(om.Group):

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

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

        self.add_subsystem('p1', ParamComp(), promotes=['*'])
        self.add_subsystem('brach', BrachistochroneODE(num_nodes=nn), promotes=['*'])


p = om.Problem()
p.driver = om.ScipyOptimizeDriver()

traj = p.model.add_subsystem('traj', dm.Trajectory(), promotes=['*'])
phase = traj.add_phase('phase0',
                       dm.Phase(ode_class=ODEGroup,
                                transcription=dm.GaussLobatto(num_segments=10)))

#p.model.add_subsystem('p1', ParamComp(), promotes=['*'])

# Add a vector parameter.
traj.add_parameter('aircraft:wing:sweep_distribution',
                   opt=False,
                   units='m',
                   #val=np.array([0.4, 0.6, 0.95]),
                   shape=(3, ),
                   static_target=True,
                   targets={'phase0': ['aircraft:wing:sweep_distribution']})

# Need to set defaults for this variable too.
#p.model.set_input_defaults('aircraft:wing:sweep_distribution',
#                           val=np.array([0.4, 0.6, 0.95]),
#                           units='m')

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)

phase.add_control('theta', units='deg',
                  rate_continuity=False, lower=0.01, upper=179.9)

phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)

# Minimize time at the end of the phase
phase.add_objective('time', loc='final', scaler=10)

phase.add_path_constraint('theta_rate', lower=0, upper=100)

p.model.linear_solver = om.DirectSolver()

p.setup()

p['traj.phase0.t_initial'] = 0.0
p['traj.phase0.t_duration'] = 2.0

p['traj.phase0.states:x'] = phase.interp('x', ys=[0, 10])
p['traj.phase0.states:y'] = phase.interp('y', ys=[10, 5])
p['traj.phase0.states:v'] = phase.interp('v', ys=[0, 9.9])
p['traj.phase0.controls:theta'] = phase.interp('theta', ys=[0.9, 101.5])

# Solve for the optimal trajectory
p.run_driver()

# Test the results
assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)

print('done')