BrachistochroneProblem

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

.

1. Abbreviation

The functional:

BR_formula1

Euler–Lagrange equations:

BR_formula2_5

.

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:

BR_formula5_1

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

BR_formula5_2

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:

BrachistochroneProblem

Figure 2-1: Brachistochrone Problem

.

According to the law of conservation of energy:

BR_formula6

Then from (6) the velocity is calculated:

BR_formula6a

Time for the ds segment:

BR_formula6_1

Summing up the time for all ds:

BR_formula7

with boundary value problem:

BR_formula8

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

BR_formula9

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

BR_formula10

Then substituting (10) into (9):

BR_formula11

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

BR_formula11_1

Next the integral is calculated:

BR_formula12

Introducing a new notation:

BR_formula13

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

BR_formula14_15

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

BR_formula16_17

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

BR_formula18_19

Example
Consider the solution of the brachistochrone problem for boundary conditions:

BR_formula20

Then using (14) – (20):

BR_formula21

Figure 2-2 shows the solution (21):

BrachistochronePlot_1_1

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

.

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

BR_formula22

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:

BR_formula23_24

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

BR_formula25_26

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

BR_formula27

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

BR_formula28

And at the last N-th step:

BR_formula29

.

3. Dynamic Programming for the Brachistochrone Problem

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

BR_formula30

From (30), (30 a) and (30b) it is possible to write:

BR_formula31

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

BR_formula32

It’s being recorded:

BR_formula33

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

BR_formula34

Formula (34) can be written in other versions:

BR_formula34ab

BR_formula35

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

BR_formula35a

Comparing (34) and (35), it can be seen that the values calculated by the formula (34) are always greater than the values from (35). Thus using (35) will lead to an incorrect solution, which will be reduced to an option:

BR_formula35_1

The error is due to the fact that an approximate linear approximation is used, so for numerical calculations we will use only formula (34) for the case u != 0.
Let’s consider a Brachystochronous Problem using dynamic programming and formulas (7), (23) – (29), (34), ( 34a,b).
See Figure 3-1:
DynamicProgrammingBrachistochroneProblem

Figure 3-1: Dynamic Programming for the Brachistochrone Problem

.

A grid is set for the x,y coordinates:

BR_formula36

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

BR_formula37

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

BR_formula38

Step 1:

BR_formula39

Note:
  • Next (in the following steps) optimal strategies are sought for the functions y(x) obtained in the first step.
.
Step 2:

BR_formula40

Step N-1:
BR_formula41
Last Step N:

BR_formula42

The Octave script Brachistochrone.m is used to implement these calculations with one of two similar functions: BrachistochroneCurveWithShortDynamicProgramming() and BrachistochroneCurveWithDynamicProgramming().
Figure 3-2 shows the solution of the problem for a 1 meter x 1 meter area with a coarse grid step:

BR_formula43

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

DynamicProgrammingForExampleOfBrachistochroneProblem

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

.

Unfortunately using a smaller grid step leads to an unstable solution to the problem, since the deviation from the exact solution (cycloid) increases.

.

4. Conclusions

  • The solution of the Brachistochrone problem by dynamic programming was obtained acceptable only for the coarse step of the x-y grid. As the grid step decreases, the error increases. It means, that the solution is unstable. Further analysis is needed to select a grid that will provide a sustainable solution.
  • 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.
  • Often the formula (25) is calculated approximately as follows:

    BR_formula44

.

5. Octave GNU file DynamicProgrammingForBrachistochrone.m

The script generates 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
% Using for the next step only with one delta: y_i+1 = y_i +/- delta_y
% Minimum of the functional: integral(sqrt((1+(y'(x))^2)/2gy(x))dx)
% One step (u_i=y_i+1-y_i)!=0: sqrt(2/g)sqrt(delta_x^2+(u_i)^2)(sqrt(y_i+1)-sqrt(y_i))/u_i)=sqrt(2/g)sqrt(delta_x^2+(y_i+1 – y_i)^2)(sqrt(y_i+1)-sqrt(y_i))/(y_i+1 – y_i))
% One step (u_i=y_i+1-y_i)==0: delta_x/sqrt(2 g y_i)
% Start Point A(0,0)
% End Point B(delta_x_in_meter*x_step_nmb,delta_y_in_meter*Y_end_index)
% Input parameters:
% delta_x_in_meter – x step in meter
% delta_y_in_meter – y step in meter
% x_step_nmb – x step number
% y_step_nmb – y step number
% Y_end_index – index of the end point B
% g = 9.81 m/s^2
% Output parameter:
% f_output – functional output
function [y_output f_output x_grid y_grid fy_nmb_opt] = BrachistochroneCurveWithShortDynamicProgramming(delta_x_in_meter, delta_y_in_meter, x_step_nmb, y_step_nmb, Y_end_index)

% Brachistochrone curve with the dynamic programming
% Using for the next step all possible deltas: y_i+1 = y_i +/- k*delta_y
% Minimum of the functional: integral(sqrt((1+(y'(x))^2)/2gy(x))dx)
% One step (u_i=y_i+1-y_i)!=0: sqrt(2/g)sqrt(delta_x^2+(u_i)^2)(sqrt(y_i+1)-sqrt(y_i))/u_i)=sqrt(2/g)sqrt(delta_x^2+(y_i+1 – y_i)^2)(sqrt(y_i+1)-sqrt(y_i))/(y_i+1 – y_i))
% One step (u_i=y_i+1-y_i)==0: delta_x/sqrt(2 g y_i)
% Start Point A(0,0)
% End Point B(delta_x_in_meter*x_step_nmb,delta_y_in_meter*Y_end_index)
% Input parameters:
% delta_x_in_meter – x step in meter
% delta_y_in_meter – y step in meter
% x_step_nmb – x step number
% y_step_nmb – y step number
% Y_end_index – index of the end point B
% g = 9.81 m/s^2
% Output parameter:
% f_output – functional output
function [y_output f_output x_grid y_grid fy_nmb_opt] = BrachistochroneCurveWithDynamicProgramming(delta_x_in_meter, delta_y_in_meter, x_step_nmb, y_step_nmb, Y_end_index)

% The Equation: x=1-cos(x)+sin(x)
% Help variable y = x/factor
% Iteration: y_next = (1-cos(factor*y_previous)+sin(factor*y_previous))/factor
% Input parameters:
% x_start – start value
% factor – new variable factor
% iteration_nmb – iteration number
% Output parameter:
% x_end – result value
function x_end = IterationFunction(x_start,factor,iteration_nmb)

.

6. Download the DynamicProgrammingForBrachistochrone.m

You can download the files:
DynamicProgrammingForBrachistochrone.m
with the button:

.

7. 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