Skip to content

How do I connect a scalar input to the ODE?

By default, we recommend that users treat all ODE input variables as if they are potentially dynamic. This allows the user to use the input as either a dynamic control, or as a static design or input parameter. By default, parameters will "fan" the value out to all nodes. This allows the partials to be defined in a consistent fashion (generally a diagonal matrix for a scalar input and output) regardless of whether the input is static or dynamic.

But there are some cases in which the user may know that a variable will never have the potential to change throughout the trajectory. In these cases, we can reduce a bit of the data transfer OpenMDAO needs to perform by defining the input as a scalar in the ODE, rather than sizing it based on the number of nodes.

The Brachistochrone with a static input.

The local gravity g in the brachistochrone problem makes a good candidate for a static input parameter. The brachistochrone generally won't be in an environment where the local acceleration of gravity is varying by any significant amount.

In the slightly modified brachistochrone example below, we add a new option to the BrachistochroneODE static_gravity that allows us to decide whether gravity is a vectorized input or a scalar input to the ODE.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
import numpy as np
import openmdao.api as om


class BrachistochroneODE(om.ExplicitComponent):

    def initialize(self):
        self.options.declare('num_nodes', types=int)
        self.options.declare('static_gravity', types=(bool,), default=False,
                             desc='If True, treat gravity as a static (scalar) input, rather than '
                                  'having different values at each node.')

    def setup(self):
        nn = self.options['num_nodes']
        g_default_val = 9.80665 if self.options['static_gravity'] else 9.80665 * np.ones(nn)

        # Inputs
        self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s')

        self.add_input('g', val=g_default_val, desc='grav. acceleration', units='m/s/s')

        self.add_input('theta', val=np.ones(nn), desc='angle of wire', units='rad')

        self.add_output('xdot', val=np.zeros(nn), desc='velocity component in x', units='m/s',
                        tags=['state_rate_source:x', 'state_units:m'])

        self.add_output('ydot', val=np.zeros(nn), desc='velocity component in y', units='m/s',
                        tags=['state_rate_source:y', 'state_units:m'])

        self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2',
                        tags=['state_rate_source:v', 'state_units:m/s'])

        self.add_output('check', val=np.zeros(nn), desc='check solution: v/sin(theta) = constant',
                        units='m/s')

        # Setup partials
        arange = np.arange(self.options['num_nodes'])
        self.declare_partials(of='vdot', wrt='theta', rows=arange, cols=arange)

        self.declare_partials(of='xdot', wrt='v', rows=arange, cols=arange)
        self.declare_partials(of='xdot', wrt='theta', rows=arange, cols=arange)

        self.declare_partials(of='ydot', wrt='v', rows=arange, cols=arange)
        self.declare_partials(of='ydot', wrt='theta', rows=arange, cols=arange)

        self.declare_partials(of='check', wrt='v', rows=arange, cols=arange)
        self.declare_partials(of='check', wrt='theta', rows=arange, cols=arange)

        if self.options['static_gravity']:
            c = np.zeros(self.options['num_nodes'])
            self.declare_partials(of='vdot', wrt='g', rows=arange, cols=c)
        else:
            self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange)

    def compute(self, inputs, outputs):
        theta = inputs['theta']
        cos_theta = np.cos(theta)
        sin_theta = np.sin(theta)
        g = inputs['g']
        v = inputs['v']

        outputs['vdot'] = g * cos_theta
        outputs['xdot'] = v * sin_theta
        outputs['ydot'] = -v * cos_theta
        outputs['check'] = v / sin_theta

    def compute_partials(self, inputs, partials):
        theta = inputs['theta']
        cos_theta = np.cos(theta)
        sin_theta = np.sin(theta)
        g = inputs['g']
        v = inputs['v']

        partials['vdot', 'g'] = cos_theta
        partials['vdot', 'theta'] = -g * sin_theta

        partials['xdot', 'v'] = sin_theta
        partials['xdot', 'theta'] = v * cos_theta

        partials['ydot', 'v'] = -cos_theta
        partials['ydot', 'theta'] = v * sin_theta

        partials['check', 'v'] = 1 / sin_theta
        partials['check', 'theta'] = -v * cos_theta / sin_theta ** 2

In the corresponding run script, we pass {'static_gravity': True} as one of the ode_init_kwargs to the Phase, and declare g as a static design variable using the dynamic=False argument.

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(initial_bounds=(0, 0), 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'], dynamic=False, 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.
#
p.set_val('traj.phase0.t_initial', 0.0)
p.set_val('traj.phase0.t_duration', 2.0)

# 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.
p.set_val('traj.phase0.states:x', phase.interpolate(ys=[0, 10], nodes='state_input'))
p.set_val('traj.phase0.states:y', phase.interpolate(ys=[10, 5], nodes='state_input'))
p.set_val('traj.phase0.states:v', phase.interpolate(ys=[0, 9.9], nodes='state_input'))

# Guesses for controls are provided at all control_input node.
# Here phase.interpolate is used to linearly interpolate values onto the control input nodes.
p.set_val('traj.phase0.controls:theta', phase.interpolate(ys=[5, 100.5], nodes='control_input'))

# Set the value for gravitational acceleration.
p.set_val('traj.phase0.parameters:g', 9.80665)

#
# Solve for the optimal trajectory
#
dm.run_problem(p, simulate=True)

# Test the results


# Generate the explicitly simulated trajectory
exp_out = om.CaseReader('dymos_simulation.db').get_case('final')

# Extract the timeseries from the implicit solution and the explicit simulation
x = p.get_val('traj.phase0.timeseries.states:x')
y = p.get_val('traj.phase0.timeseries.states:y')
t = p.get_val('traj.phase0.timeseries.time')
theta = p.get_val('traj.phase0.timeseries.controls:theta')

x_exp = exp_out.get_val('traj.phase0.timeseries.states:x')
y_exp = exp_out.get_val('traj.phase0.timeseries.states:y')
t_exp = exp_out.get_val('traj.phase0.timeseries.time')
theta_exp = exp_out.get_val('traj.phase0.timeseries.controls: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()
Full total jacobian was computed 3 times, taking 0.090680 seconds.
Total jacobian shape: (40, 51) 


Jacobian shape: (40, 51)  (13.19% nonzero)
FWD solves: 8   REV solves: 0
Total colors vs. total size: 8 vs 51  (84.3% improvement)

Sparsity computed using tolerance: 1e-25
Time to compute sparsity: 0.090680 sec.
Time to compute coloring: 0.001454 sec.
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 18.016167304638337
            Iterations: 24
            Function evaluations: 24
            Gradient evaluations: 24
Optimization Complete
-----------------------------------

Simulating trajectory traj
Done simulating trajectory traj