Skip to content

The Brachistochrone

Things you'll learn through this example

  • How to define a basic Dymos ODE system.
  • How to test the partials of your ODE system.
  • Adding a Trajectory object with a single Phase to an OpenMDAO Problem.
  • Imposing boundary conditions on states with simple bounds via fix_initial and fix_final.
  • Using the Phase.interpolate` method to set a linear guess for state and control values across the Phase.
  • Checking the validity of the result through explicit simulation via the Trajectory.simulate method.

The brachistochrone is one of the most well-known optimal control problems. It was originally posed as a challenge by Johann Bernoulli.

The brachistochrone problem

Given two points A and B in a vertical plane, find the path AMB down which a movable point M must by virtue of its weight fall from A to B in the shortest possible time.

  • Johann Bernoulli, Acta Eruditorum, June 1696

We seek to find the optimal shape of a wire between two points (A and B) such that a bead sliding without friction along the wire moves from point A to point B in minimum time.

Brachistochrone free body diagram

State variables

In this implementation, three state variables are used to define the configuration of the system at any given instant in time.

  • x: The horizontal position of the particle at an instant in time.
  • y: The vertical position of the particle at an instant in time.
  • v: The speed of the particle at an instant in time.

System dynamics

From the free-body diagram above, the evolution of the state variables is given by the following ordinary differential equations (ODE).

\begin{align} \frac{d x}{d t} &= v \sin(\theta) \\ \frac{d y}{d t} &= -v \cos(\theta) \\ \frac{d v}{d t} &= g \cos(\theta) \end{align}

Control variables

This system has a single control variable.

  • \theta: The angle between the gravity vector and the tangent to the curve at the current instant in time.

The initial and final conditions

In this case, starting point A is given as (0, 10). The point moving along the curve will begin there with zero initial velocity.

The initial conditions are:

\begin{align} x_0 &= 0 \\ y_0 &= 10 \\ v_0 &= 0 \end{align}

The end point B is given as (10, 5). The point will end there, but the velocity at that point is not constrained.

The final conditions are:

\begin{align} x_f &= 10 \\ y_f &= 5 \\ v_f &= \mathrm{free} \end{align}

Defining the ODE as an OpenMDAO System

In Dymos, the ODE is an OpenMDAO System (a Component, or a Group of components). The following ExplicitComponent computes the state rates for the brachistochrone problem.

More detail on the workings of an ExplicitComponent can be found in the OpenMDAO documentation. In summary:

  • initialize: Called at setup, and used to define options for the component. ALL Dymos ODE components should have the property num_nodes, which defines the number of points at which the outputs are simultaneously computed.
  • setup: Used to add inputs and outputs to the component, and declare which outputs (and indices of outputs) are dependent on each of the inputs.
  • compute: Used to compute the outputs, given the inputs.
  • compute_partials: Used to compute the derivatives of the outputs w.r.t. each of the inputs analytically. This method may be omitted if finite difference or complex-step approximations are used, though analytic is recommended.
 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

Things to note about the ODE system

  • There is no input for the position states (x and y). The dynamics aren't functions of these states, so they aren't needed as inputs.
  • While g is an input to the system, since it will never change throughout the trajectory, it can be an option on the system. This way we don't have to define any partials w.r.t. g.
  • The output check is an auxiliary output, not a rate of the state variables. In this case, optimal control theory tells us that check should be constant throughout the trajectory, so it's a useful output from the ODE.

Testing the ODE

Now that the ODE system is defined, it is strongly recommended to test the analytic partials before using it in optimization. If the partials are incorrect, then the optimization will almost certainly fail. Fortunately, OpenMDAO makes testing derivatives easy with the check_partials method. The assert_check_partials method in openmdao.utils.assert_utils can be used in test frameworks to verify the correctness of the partial derivatives in a model.

The following is a test method which creates a new OpenMDAO problem whose model contains the ODE class. The problem is setup with the force_alloc_complex=True argument to enable complex-step approximation of the derivatives. Complex step typically produces derivative approximations with an error on the order of 1.0E-16, as opposed to ~1.0E-6 for forward finite difference approximations.

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)
-----------------------------------
Component: BrachistochroneODE 'ode'
-----------------------------------
'<output>' wrt '<variable>' | calc mag.  | check mag. | a(cal-chk) | r(cal-chk)
-------------------------------------------------------------------------------

'check'    wrt 'theta'      | 7.8782e+01 | 7.8782e+01 | 1.1102e-16 | 1.4092e-18
'check'    wrt 'v'          | 6.4046e+00 | 6.4046e+00 | 6.2804e-16 | 9.8060e-17
'vdot'     wrt 'theta'      | 1.3160e+01 | 1.3160e+01 | 2.0351e-15 | 1.5464e-16
'vdot'     wrt 'v'          | 0.0000e+00 | 0.0000e+00 | 0.0000e+00 | nan
'xdot'     wrt 'theta'      | 6.2973e+00 | 6.2973e+00 | 2.4825e-16 | 3.9422e-17
'xdot'     wrt 'v'          | 1.3419e+00 | 1.3419e+00 | 6.2063e-17 | 4.6249e-17
'ydot'     wrt 'theta'      | 4.4741e+00 | 4.4741e+00 | 2.5438e-16 | 5.6857e-17
'ydot'     wrt 'v'          | 1.7886e+00 | 1.7886e+00 | 1.1102e-16 | 6.2071e-17

##################################################################
Sub Jacobian with Largest Relative Error: BrachistochroneODE 'ode'
##################################################################
'<output>' wrt '<variable>' | calc mag.  | check mag. | a(cal-chk) | r(cal-chk)
-------------------------------------------------------------------------------
'vdot'     wrt 'theta'      | 1.3160e+01 | 1.3160e+01 | 2.0351e-15 | 1.5464e-16

Solving the Problem

The following script fully defines the brachistochrone problem with Dymos and solves it. In this section we'll walk through each step.

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

p['traj.phase0.states:x'] = phase.interpolate(ys=[0, 10], nodes='state_input')
p['traj.phase0.states:y'] = phase.interpolate(ys=[10, 5], nodes='state_input')
p['traj.phase0.states:v'] = phase.interpolate(ys=[0, 9.9], nodes='state_input')
p['traj.phase0.controls:theta'] = phase.interpolate(ys=[5, 100.5], nodes='control_input')

#
# Solve for the optimal trajectory
#
dm.run_problem(p)

# Test the results


# Generate the explicitly simulated trajectory
exp_out = traj.simulate()

plot_results([('traj.phase0.timeseries.states:x', 'traj.phase0.timeseries.states:y',
               'x (m)', 'y (m)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.controls:theta',
               'time (s)', 'theta (deg)')],
             title='Brachistochrone Solution\nHigh-Order Gauss-Lobatto Method',
             p_sol=p, p_sim=exp_out)

plt.show()
Full total jacobian was computed 3 times, taking 0.073863 seconds.
Total jacobian shape: (40, 51) 


Jacobian shape: (40, 51)  (18.92% nonzero)
FWD solves: 12   REV solves: 0
Total colors vs. total size: 12 vs 51  (76.5% improvement)

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

Simulating trajectory traj
Done simulating trajectory traj