Solving the Brachistochrone Problem using Dynamic Programming and the GNU Octave Tool

Solving the Brachistochrone Problem using Dynamic Programming and the GNU Octave Tool

.

1. Abbreviation

The functional:

Euler–Lagrange equations:

.

2. Introduction

A functional is an operator that displays (calculates) a number for each function y(x). In this article I will use a functional of the form (1). Calculating the integral (1) for some function y(x), we get a number. The task of the calculus of variations is to find the function y(x), which gives the minimum (maximum) for the functional (1).

Note
Multidimensional functions are also used in problems of variations calculus:

and of course more complex functionals. In this article only one-dimensional functions y(x) and a functional of the form (1) will be considered.

In the theory of variations calculus a necessary condition has been proved for the function y(x), which gives a minimum (maximum) for (1). This is the Euler–Lagrange formula (2). If the function F() does not depend on x, then the Euler–Lagrange formula is simplified and takes the form (3).

Now let’s look at an example of a variational problem: the Brachistochrone Problem. This is one of the first variational problems proposed and solved by Johann Bernoulli in 1696. It is required to determine which line the material point will roll down in the shortest time from the point O(0,0) to the point

at the same time the friction and resistance of the medium are neglected. The initial velocity of the point at point O(0,0) is zero. See the Figure 2-1:

Figure 2-1: Brachistochrone Problem

.

According to the law of conservation of energy:

Then from (6) the velocity is calculated:

Time for the ds segment:

Summing up the time for all ds:

with boundary value problem:

Formula (7) represents the functional that needs to be minimized. To find y(x), use the Euler–Lagrange formula (3):

Let’s introduce a new parameter r using cotangent function:

Then substituting (10) into (9):

Using (10) and (11) can be written:

Next the integral is calculated:

Introducing a new notation:

Then (11) and (12) will take the form:

Thus formulas (14) and (15) give a general solution to the variational brachystochrone problem. The resulting curve is a cycloid. From boundary conditions (8):

Solving equation (17) with respect to r_end, it is possible to determine the constant C1:

Example

Consider the solution of the brachistochrone problem for boundary conditions:

Note: This example (20) is taken from [8]

Using the (14) – (20):

Figure 2-2 shows the solution (21):

Figure 2-2: Brachistochrone Problem with Boundary Condition: y(0)=0 and y(2)=1.2

.

The optimal time is determined by the formula (7):

Most problems of the calculus of variations are not solved analytically and require a numerical solution. In this article I want to solve the Brachistochrone problem numerically using dynamic programming methods and compare the result with the exact solution (21). I have already written about dynamic programming in my article: “Dynamic Programming and Traveling Salesman Problem Using Octave GNU Tool”. See also [1], [2].

Dynamic programming is the search for the optimal solution to a given problem, which can be divided into many steps and search for the optimal solution only for the current step, which greatly simplifies calculations. Sometimes we can talk about a purposeful search of possible variants with an assessment of their optimality, while the complexity of the task decreases compared to a complete search of all variants.

We will look for the function y(x) given at discrete points for this we will define a grid on the x-y plane.

X-grid (0…Nx-1):

Y-grid (0…Ny-1):

Then we have a process of (Nx – 1) steps and the win function (or cost function) at each step (i+1) is calculated from the functional (1):

The calculation is performed using the formula from Bellman (27):

The function f(x) is calculated in the first step using the formula:

And at the last (Nx – 1) step:

Moving along the curve from the endpoint y_end(x_end) back, taking into account the optimal values of the function f(), we obtain the optimal curve y(x). See the code in the DynamicProgrammingForBrachistochrone.m and BrachistochroneVsCopilotAttempt.c/h

.

3. Dynamic Programming for the Brachistochrone Problem

We introduce the following notation using linear approximation, see the [3]:

From (30b) it is possible to write:

Then the linear approximation of y(x) on the segment (interval):

It’s being recorded:

Let’s write down the functional (7) for the (i+1) step using the linear approximation of (31) and (33):

Formula (34) can be written in other versions:

Formula (35) can be written in a slightly different form:

Let’s consider a Brachystochronous Problem using dynamic programming and formulas (7), (23) – (29), ( 34a).

See the Figure 3-1:

Figure 3-1: Dynamic Programming for the Brachistochrone Problem

.

A grid is set for the x,y coordinates:

The boundary conditions are set starting point O and ending point B:

In the first step, it is calculated using formulas (28) and (34a) for all values of the y grid:

Note:
The initial point is excluded from consideration, since at this point the velocity is zero which leads to division by zero in (34a).

Step 1:

.

Step 2:

Step Nx – 2:

 

Last Step Nx – 1:

The Figure 3-2 shows the solution of the problem (20) for the area x = 0 … 2 meter and y = 0 … 1.2 meter with a grid step:

The error of the dynamic programming solution for (43) is 0.5%. The movement of the cargo along the trajectory obtained by dynamic programming (red curve) will take 0.5% longer than the ideal cycloid (blue curve).

Figure 3-2: Dynamic Programming for the Example of the Brachistochrone Problem

.

See the DynamicProgrammingForBrachistochrone.m and BrachistochroneVsCopilotAttempt.c/h for the more implementation information

.

4. Conclusions

• The complexity of the dynamic programming calculation is determined by the number of combinations that to be checked/analyzed. At the first and last step, (Ny – 1) options are calculated and for the remaining intermediate steps (Ny – 1)^2. The number of combinations can be calculated with the formula:

  

  Using the (44) for (43):

  

  Note:
  A complete direct search of all the options can be calculated:

  

  Using the (45) for (43):

  

  Thus dynamic programming provides great computational benefits!

• To more accurately solve the problem using dynamic programming with boundary condition (8), you can introduce an additional constraint with a Lagrange multiplier/coefficient or define an additional functional that will ensure that the optimal solution reaches the endpoint corresponding to the boundary condition. A similar technique is used in the article “Solving the Brachistochrone Problem using Pontryagin’s Method and the GNU Octave Tool”.
  

• For any problem of variations calculus, it is necessary to prove the existence of an optimal solution. The existence of a solution has been proven for the Brachistochrone problem.

• In [3], [4] a numerical solution of the Brachistochrone problem is given using additional techniques.

.

5. Octave GNU file DynamicProgrammingForBrachistochrone.m

The script genrates a solution to the Brachistochrone problem using the numerical method of dynamic programming and the analytical method of variations calculus.

% Brachistochrone curve with the dynamic programming
% Minimum of the functional: integral(sqrt((1+(y'(x))^2)/2gy(x))dx)
% One step: v_average = sqrt(g/2)*(sqrt(y_i+1) + sqrt(y_i))
% t_i+1 = sqrt(delta_x_in_meter^2+((y_i+1)-(y_i))^2)/v_average
% Start Point A(0,0)
% End Point B(delta_x_in_meter*(x_point_nmb-1),delta_y_in_meter*(Y_end_index-1))
% Input parameters:
% delta_x_in_meter – x step in meter
% delta_y_in_meter – y step in meter
% x_point_nmb – x point number
% y_point_nmb – y point number
% Y_end_index – index of the end point B
% Output parameter:
% y_output – optimal output
% f_output – functional output
% x_grid – x axis
% y_grid – y axis
function [y_output f_output x_grid y_grid] = BrachistochroneCurveWithDynamicProgramming(delta_x_in_meter, delta_y_in_meter, x_point_nmb, y_point_nmb, Y_end_index)

% The Equation: x=1-cos(x)+sin(x)
% Iteration: x_next = (1-cos(x_previous)+sin(x_previous))
% Input parameters:
% x_start – start value
% iteration_nmb – iteration number
% Output parameter:
% x_end – result value
function x_end = IterationFunction(x_start,iteration_nmb)

.

6. C Language BrachistochroneVsCopilotAttempt.c/h

The code was generated by artificial intelligence (AI): Copilot, Claude 3.5, Sonnet model. Then the bugs were fixed and the code was optimized.

/*
The function for calculating the moving time between two points
Input
float y1 – y ordinate of the point 1
float y2 – y ordinate of the point 2
float dx – delta x
Return
float – moving time from y1 to y2 and dx
*/
float calculateTime(float y1, float y2, float dx);

/*
The function for solving the brachystochrone problem
Input
point_t * start – start point
point_t * end – end point
float * x_grid – x grid
float * y_grid – y grid
float * y_output – optimal curve
Return
float – optimal time from start point to end point
*/
float solveBrachistochrone(point_t * start, point_t * end, float * x_grid, float * y_grid, float * y_output);

.

7. Download the DynamicProgrammingForBrachistochrone.m and BrachistochroneVsCopilotAttempt.c/h

You can download the files:

DynamicProgrammingForBrachistochrone.m
BrachistochroneVsCopilotAttempt.c/h

with the button:

.

8. Literature / References

[1] R. Bellman “Dynamic Programming“, Princeton University Press, New Jercy, 1957

[2] R. Bellman, S. Dreyfus “Applied Dynamic Programming“, Princeton University Press, New Jercy, 1962

[3] J. Vomlel “Solving the Brachistochrone Problem by an Influence Diagram”, Institute of Information Theory and Automation, Czech Academy of Sciences, 2017

[4] Christopher Raphaely, “Coarse-to-Fine Dynamic Programming”, 2001

[5] L. Elsgolz “Differential equations and calculus of variations”, Publishing House: Nauka, Moscow, 1969

[6] L. Zlaf “Calculus of variations and integral equations”, Reference guide, Publishing House: Nauka, Moscow, 1966

[7] L. Rozov “Optimal Control”, Digest, Publishing House: Znanie, Moscow, 1978

[8] R. P. Fedorenko “Approximate solution of optimal control problems”, Publishing House: Nauka, Moscow, 1978

.

Appendix. Bellman’s Partial Differential Equations

Let there be a functional:

Further:

As a result:

The equation (A1-2) can be represented as two equations:

At the same time:

Equations (A1-2) as well as (A1-3)-(A1-5) are called Bellman’s partial differential equations

.

Notes:

• The Euler-Lagrange equation is easily obtained from formulas (A1-3) and (A1-4).

  To do this, take the partial derivative to y from (A1-3):
  
  And differentiate (A1-4) by x:
  
  Next subtracting equation A(1-7) from A(1-6) and the Euler–Lagrange equation will be obtained:
  

• See the [2] for more information about it

.