Brachistochrone

Background

The brachistochrone asks for the curve along which a bead, released from rest and moving without friction under uniform gravity, reaches a lower point in the least time. Its solution is a cycloid rather than a straight line, making it a standard test for optimal-control and calculus-of-variations methods.

Problem formulation

Here \(y\) is positive downward and the path angle \(\theta\) is measured from the downward vertical. The minimum-time problem is

\[ \begin{aligned} \min_{x,y,v,\theta,t_f}\quad & t_f=\int_0^{t_f}1\,\mathrm{d}t \\ \text{subject to}\quad & \dot{x}=v\sin\theta, \\ & \dot{y}=v\cos\theta, \\ & \dot{v}=g\cos\theta, \qquad g=9.81\ \mathrm{m/s^2}, \\ & [x(0),y(0),v(0)]=[0,0,0], \\ & [x(t_f),y(t_f)]=[2,2], \\ & v\ge0, \qquad 0\le\theta\le\frac{\pi}{2}. \end{aligned} \]

The final speed and final time are free.

Variables and units

Symbol	Meaning	Unit
\(t\)	Time	s
\(x\)	Horizontal displacement	m
\(y\)	Downward displacement	m
\(v\)	Speed along the curve	m/s
\(\theta\)	Path angle from the downward vertical	rad
\(g\)	Gravitational acceleration	\(9.81\ \mathrm{m/s^2}\)

Modeling choices

Using downward displacement keeps gravitational acceleration positive in the model. The speed constraint excludes a nonphysical reversal of the path parameter, and the angle bounds select monotone motion to the right and downward. The initial guess follows a straight descending path with speed \(v=\sqrt{2gy}\), which is consistent with conservation of mechanical energy.

Run the example

From the repository root, run:

python -m examples.brachistochrone

To save the figure without opening a window:

python -m examples.brachistochrone --save brachistochrone.png --no-show

Key implementation

import numpy as np
import sympy as sp

from pockit.lobatto import System, linear_guess

system = System(0)
phase = system.new_phase(["x", "y", "speed"], ["path_angle"])
_, _, speed = phase.x
(path_angle,) = phase.u

phase.set_dynamics(
    [
        speed * sp.sin(path_angle),
        speed * sp.cos(path_angle),
        9.81 * sp.cos(path_angle),
    ]
)
phase.set_integral([1.0])
phase.set_phase_constraint(
    [speed, path_angle], [0.0, 0.0], [np.inf, np.pi / 2.0]
)
phase.set_boundary_condition([0.0, 0.0, 0.0], [2.0, 2.0, None], 0.0, None)
phase.set_discretization(10, 8)
system.set_phase([phase])
system.set_objective(phase.I[0])

guess = linear_guess(phase, 0.0)
guess.t_f = 1.0

The complete example solves with Ipopt, constructs the analytical cycloid from the endpoint geometry, and checks the numerical travel time against it.

Verified result

Pockit gives a minimum travel time of \(0.824338670697\ \mathrm{s}\). The analytical cycloid gives \(0.824338669439\ \mathrm{s}\), a difference of about \(1.3\times10^{-9}\ \mathrm{s}\). The plotted numerical path is visually coincident with the cycloid.

Source code

See the complete runnable example: examples/brachistochrone.py.