Solving the Brachistochrone Problem using Pontryagin’s Method and the GNU Octave Tool
.
1. Introduction
.
1.1. The Brachistochrone Problem
Let’s consider 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 possible time from the point:
to the point
The friction and resistance of the medium are neglected. The initial velocity of the point O(0,0) is zero. See the Figure 1-1:
Figure 1-1: Brachistochrone Problem
.
According to the law of conservation of energy:
Then from (1) the velocity is calculated:
Time for the ds section:
Summing up the time over all ds:
with boundary value problem:
Formula (2) represents the functional that needs to be minimized by selecting the function y(x) for this purpose.
Notes:
-
A functional is an operator that assigns (calculates) a number to each function y(x)
-
Multidimensional functions are also used in calculus of variations problems:
and of course more complex functionals. In this article, we will consider the case of a one-dimensional function y(x)
.
1.2. Classical Calculus of Variations and the Solution of the Brachistochrone Problem
The functional:
Boundary value problem:
Euler–Lagrange equations:
Notes:
-
In this article I will use a functional of the form (5). Calculating the integral (5) for some function y(x) we get a number. The task of the calculus of variations is to find a function y(x) that gives a minimum for the functional (5)
-
Obviously the (2) functional of the Brachistochrone Problem corresponds to (5)
-
In the theory of calculus of variations computation, a necessary condition is proved for the function y(x), which gives a minimum for (5). This is the Euler–Lagrange formula (7). If the function Φ(y, y’) does not depend on x, then the Euler–Lagrange formula is simplified and takes the form (8)
A detailed solution to the Brachistochrone Problem using classical calculus of variations is given in the article “Solving the Brachistochrone Problem using Dynamic Programming and the GNU Octave Tool”. Here I will give the finished result:
Thus (12) and (13) provide a general solution to the Brachistochrone problem. The resulting curve is a cycloid.
From the boundary conditions (3):
By solving equation (15) with respect to r_end, the constant C1 can be calculated:
Example:
Let us consider the solution of the brachistochrone problem for boundary conditions:
Note: the example (18) is taken from [2]
Then using (12) – (17):
Figure 1-2 shows the (19) solution:
Figure 1-2: Brachistochrone Problem with Boundary Condition: y(0)=0 and y(2)=1.2
.
The optimal time is determined by the formula (2):
.
1.3. Pontryagin’s Method
Let us now reformulate the problem of the calculus of variations into the language of the optimal control.
The functional in this case will take the form:
where u(x) = y’(x) – a function called control
The system is described by the differential equation:
Boundary value problem:
The task is to find a control u(x) that provides a minimum of functional (21) and hits y(x) at the endpoint (24). Thus we select the slope of the tangent to the function at each point of the graph y(x), which provides the optimal solution.
The principle of the Pontryagin method.
Let there be an optimal control u(x)=u_opt(x), which corresponds to the optimal trajectory y(x)=y_opt(x). In this case it is necessary to have such a nonzero function ψ(x), which is the solution of the differential equation:
using condition:
An auxiliary function H(ψ,y,u) is introduced, called the Hamiltonian of the Task:
This function reaches a minimum for each value of x (x_0 ≤ x ≤ x_end) at u(x)=u_opt(x) and the corresponding y(x)=y_opt(x):
It can also be proved (see the [3]) that for an optimal solution:
Notes:
-
To solve the problem the (21) – (29) are used
-
The above theorem is a necessary condition. Therefore before solving the task, it is to prove the existence of an optimal solution
-
In the classical theory of the Pontryagin’s method (see the [3], [4]) speaks of the Pontryagin maximum. Here I use the function (27), written in other form (see the [1]), which leads to a Hamiltonian minimum
-
Let’s write down equations (21) – (29) for the brachistochrone problem:
.
1.4. Gradient Method for Solving the Pontryagin’s Minimum Problem
Let’s solve problem (21) – (24) using the gradient method to find optimal control. See the [2] for a detailed method description. To do this, the Frechet derivative of the functional shall be calculated:
it can be proved for the (21) functional (see the [2], pages 29-34):
Thus for (21):
It is also convenient to represent the boundary condition (24) in the form of the (24a) functional:
From (25) we find ψ_1(x):
Then:
From (26) and (44):
Then:
Using (39) the Frechet derivative:
Using w0(x), w1(x) and knowing the perturbation of the control element δu, the perturbation of the functional δF can be calculated:
The (49) can be rewritten using the (40) and (47):
Now let’s choose the perturbation δu:
Substituting (51) into (50):
Then λ can be calculated:
Task solution plan (21) – (24):
-
Choose an arbitrary solution u(x), y(x), which translates the system from (x_0, y_0) to (x_end, y_end)
-
Solving equation (25) with (26), we find ψ(x). Next the functions u(x), y(x) are calculated
-
Next the w0(x) shall be calculated using the formula (39a)
-
The λ – Lagrange coefficient shall be calculated by the formula (53)
-
Calculate the R(s) function by formula (54) for various s, and then find the minimum of R(s) by formula (55)
In this case the δu(x) is calculated by the formula (51)
-
Knowing the optimal u_opt(x) = u(x) – s_opt δu(x), the optimal function y_opt(x) should be calculated by solving equation (22) using (23)
Note
The descent procedure is sometimes repeated several times. This article uses only one gradient descent for simplicity
.
2. Pontryagin’s Minimum for the Brachistochrone Problem using Gradient Method
Let us now solve the problem of Brachistochrone (30) – (37) with boundary conditions (18). Let’s write out once again the main equations of the problem:
Now let’s choose an arbitrary solution u(x), y(x), which translates the system from (0, 0) to (2, 1.2):
For this solution (70) we can write:
Now we will look for a solution for various s = 1; 1.4; 1.52; 1.6; 2
Table 2-1: Brachistochrone Problem with y(0)=0, y(2)=1.2 using Gradient Method
.
The exact solution of the problem is t = 0.8006. See (20). Our approximate solution gives t = 0.8051s.
Approximation error:
Figure 2-1: Brachistochrone Problem with Boundary Condition: y(0)=0 and y(2)=1.2 and Functional R(s)
.
Figure 2-2: Brachistochrone Problem with Boundary Condition: y(0)=0 and y(2)=1.2. Pattern and Approximation
.
3. Improving the Accuracy of the Calculation by Eliminating the Singular Point from the Solution
The obtained accuracy of the solution in the previous section is low, since at the start point O(0,0) the value of Φ(x, y, y’) in formula (4) is not defined. If the point y => 0 the value of the Φ(x, y, y’) function goes to infinity, which leads to an error in the approximate solution. Let’s exclude this singular point from consideration. Using the exact solution of the problem (19) and Figure 1-2, we select the other starting point:
The exact solution of the problem is calculated using the formulas:
The optimal time is determined by the formula (2):
Now let’s repeat the calculations given in section 2.
.
3.1 Pontryagin’s Minimum for the Brachistochrone Problem using Gradient Method and Start Point O(0.05, 0.183088)
Thus we have the task with boundary conditions:
The functional will look like:
Now let’s choose an arbitrary solution u(x), y(x), which translates the system from (0.05, 0.183088) to (2, 1.2):
Now we will look for a solution for various: s = 1; 2; 3; 3.4; 4; 5
Table 3-1: Brachistochrone Problem with y(0.05)=0.183088, y(2)=1.2 using Gradient Method
.
The exact solution of the problem is t = 0.6021s. See the (72). Our approximate solution gives t = 0.6023s. Approximation error:
Note
The accuracy of calculations has increased by more than an order of magnitude after the exclusion of the singularity point!
Figure 3-2: Brachistochrone Problem with Boundary Conditions: y(0.05)=0.183088, y(2)=1.2. Pattern and Approximation of the y(x)
.
3.2. Numerical Pontryagin’s Minimum for the Brachistochrone Problem using Gradient Method and Start Point O(0.05, 0.183088)
To solve the problem numerically, we must move from continuous functions to discrete ones. The functions u(x), y(x) and ψ(x):
Using a linear approximation between discrete values of functions (see the [2]):
Notes:
-
The Octave script PontryaginsMaximumForBrachistochroneArticle.m is used to numerically solve the problem
-
The solution (75) is chosen as the initial one
-
The integral functional (74) is calculated by the trapezoid method. See the Octave script function BrachistochroneFunctional0( )
-
The calculation of ψ(x) is performed in the opposite direction along the x axis: from x_end to x_0, using formulas (61) and (62). See the Octave script function: BrachistochronePsiW0Function( )
-
Calculation of w0(x) is performed according to the formula (63). See the Octave script function: BrachistochronePsiW0Function( )
-
The value of λ is calculated by the formula (67). The integral is calculated using the trapezoid method. See the Octave script function: BrachistochronePsiW0Function( )
-
The (68), (69) are calculated for s = -1:0.1:10 and the optimal ones are found: s_opt, u_opt(x) and y_opt(x). See the Octave script function: Brachistochrone_MinimumSearch( ). See also the Figure 3-3: s_opt=3.4 and Figure 3-4: y_opt(x)
Figure 3-3: Brachistochrone Problem with Boundary Condition: y(0.05)=0.183088 and y(2)=1.2: Functional R(s)
.
-
The numerical solution approximately coincides with the results of the previous section. See the Figure 3-2 and Figure 3-4
Figure 3-4: Brachistochrone Problem with Boundary Conditions: y(0.05)=0.183088, y(2)=1.2. Pattern and Approximation of the y(x)
.
4. Octave GNU file PontryaginsMethodForBrachistochroneArticle.m
The script generates graphs for the article, as well as solving the Brachistochrone problem using the numerical Pontryagin’s Method.
% Safe Division with the space: -realmax … realmax
% Input parameters:
% dividend
% divisor
% Output parameter:
% quotient – dividend/divisor
function quotient = safe_division(dividend, divisor)
% Calculate the functional: F0(u(.))=integral(Φ(y(x),u(x))dx)=integral(sqrt((1+(u(x))^2)/2gy(x))dx)
% ψ – psi, Φ – phi, δ – delta, λ – lambda
% Φ(y(x),u(x))=sqrt((1+(u(x))^2)/2gy(x))
% y'(x)=dy/dx=f(y,x)=u(x)
% For the one interval (trapezoidal rule):
% F0_i=0.5*sqrt(1/2g)*(sqrt((1+(u_i)^2)/y_i)+sqrt((1+(u_i-1)^2)/y_i-1))
% Input parameters:
% y_func – y function
% u_func – u function
% x_val – x values
% Output parameter:
% F0_value – Functional0 output
function F0_val = BrachistochroneFunctional0(y_func, u_func, x_val)
% Calculate: the ψ(x), w0(x), lambda
% ψ – psi, Φ – phi, δ – delta, λ – lambda
% psi(end) = 0
% dpsi_dx(n) = sqrt(1 + (u(n))^2)/(2*sqrt(2*g)*(y(n))^(3/2))
% psi_res(n) = psi(n-1) + delta_x_in_meter*(dpsi_dx(n) + dpsi_dx(n-1))/2
% w0_res = u(n)/sqrt(2*g*y(n)*(1+(u(n))^2))+psi_res(n)
% lambda = (1/x_end)*integral_from_0__to_x_end (w0*dx) with trapezoidal rule
% Input parameters:
% y_func – y function
% u_func – u function
% x_val – x values
% delta_x_in_meter – x interval in meter
% Output parameters:
% psi_res – ψ array result
% w0_res – w0 array
function [psi_res, w0_res, lambda] = BrachistochronePsiW0Function(y_func, u_func, x_val, delta_x_in_meter)
% Calculate: the y_next_func, u_next_func
% ψ – psi, Φ – phi, δ – delta, λ – lambda
% u_next(x)=u(x)-s*(w0_res(x)-λ)
% y_next(x)=integral(u_next(x))
% Input parameters:
% y_func – y function
% u_func – u function
% w0_res – gradient function
% lambda – λ Lagrange coefficient
% s_val – s parameter for s*δu=s*(w0_res-λ)
% delta_x_in_meter
% Output parameters:
% y_next_func – next next iteration of the y function
% u_next_func – iteration of the u function
function [y_next_func, u_next_func] = Brachistochrone_Y_U_Function(y_func, u_func, w0_res, lambda, s_val, delta_x_in_meter)
% Calculate: y_func_array, u_func_array, F0_array, optimal index in array
% ψ – psi, Φ – phi, δ – delta, λ – lambda
% Input parameters:
% y_func – y function
% u_func – u function
% w0_res – gradient function
% s_val – s parameter for the s*δu=s*(w0_res-λ)
% delta_x_in_meter
% Output parameters:
% y_func_array – y function array for each s value
% u_func_array – u function array for each s value
% x_val – x values
% F0_array – F0 functional array for each s value
% s_array – s values
% optim_index – optimal index in array (minimal functional)
function [y_func_array, u_func_array, F0_array, optim_index] = Brachistochrone_MinimumSearch(y_func, u_func, x_val, s_array, delta_x_in_meter)
% Solving the equation by the iterative method
% The Equation: x=wait_value*2/1.2025 + sin(x)
% Iteration: x_next = (wait_value*2/1.2025 + sin(x_previous))
% Input parameters:
% x_start – start value
% iteration_nmb – iteration number
% wait_value – wait value
% Output parameter:
% x_end – result value
function x_end = IterationFunction(x_start,iteration_nmb,wait_value)
.
5. Download the PontryaginsMethodForBrachistochroneArticle.m
You can download the file:
PontryaginsMethodForBrachistochroneArticle.m
with the button:
.