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:
Then using (14) – (20):
Figure 2-2 shows the solution (21):
Figure 2-2: Brachistochrone Problem with Boundary Condition: y(0)=0 and y(1)=1
.
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:
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):
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 N-th step:
.
3. Dynamic Programming for the Brachistochrone Problem
We introduce the following notation using linear approximation, see the [3]:
From (30), (30 a) and (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:
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:
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:
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:
Step 1:
Note:
-
Next (in the following steps) optimal strategies are sought for the functions y(x) obtained in the first step.
.
Step 2:
…
Step N-1:
Last Step N:
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:
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).
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:
.
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:
.