matlab research example
Written by Daniel Herber on June 28, 2017.
---
### Problem Definition
This classic calculus of variations problem is stated as [[link]](http://www.math.uiuc.edu/documenta/vol-ismp/11_knobloch-eberhard-3.pdf):
> "Given two points A and B in a vertical plane, find the path AMB down which a movable point M must by virtue of its weight fall from A to B in the shortest possible time."
The problem can be posed as a dynamic optimization problem [[link]](http://www.gpops2.com/Examples/Brachistochrone.html):
<div>
\begin{align}
\min_{u(t)} \quad & t_f \\
\text{subject to:} \quad & \dot{x} = v \sin(u) \\
& \dot{y} = v \cos(u) \\
& \dot{v} = g \cos(u) \\
& x(0) = y(0) = v(0) = 0 \\
& x(t_f) = x_f > 0 \\
& y(t_f) = y_f > 0 \\
\end{align}
</div>
where $u(t)$ parameterizes the descent curve. Note that the true $y$ needs to be reflected about the $x$-axis.
---
### Parametric Solution
It can be shown that the solution to this problem can be expressed as parametric equations that define a cycloid (see [[link]](http://mathworld.wolfram.com/BrachistochroneProblem.html)):
$$ x(\theta) = \frac{1}{2} k^2 (\theta - \sin(\theta)), \qquad y(\theta) = \frac{1}{2} k^2 (1 - \cos(\theta)) $$
where $k>0$ and $\theta_f \in (0, 2\pi]$ are constants that can be numerically determined by solving the two final boundary conditions on $x$ and $y$. Unfortunately, this solution is lacking since it is not described in $t$ nor does it provide the trajectories for $v$ and $u$.
---
### Scaling the Solution
To find time-based solution to this problem, we will utilize scaling (see [[link]]() for a recent paper on scaling in dynamic optimization). The cycloid trajectories are correct with respect to linear scaling of the time horizon. We introduce this relationship with:
$$ t = \alpha_t \theta $$
where $\alpha_t$ is a positive constant that needs to be determined. The chain rule and linearity of differentiation helps provide the relationship between the derivatives:
<div>
\begin{align}
\frac{dx(t)}{dt} &= \frac{1}{\alpha_t} \frac{dx(\theta)}{d\theta} \\
\frac{d^2x(t)}{dt^2} &= \frac{1}{\alpha_t^2} \frac{d^2x(\theta)}{d\theta^2}
\end{align}
</div>
A similar procedure can be performed with $y$. We also are concerned the effect of scaling on the magnitude of the velocity and acceleration. These are defined as (noting that $| \vec{v}(t) | := v$ above):
<div>
\begin{gather}
\vec{v}(t) = \left( \frac{dx}{dt}, \frac{dy}{dt} \right), \quad | \vec{v}(t) | = \sqrt{ \left( \frac{dx}{dt} \right)^2 + \left( \frac{dy}{dt} \right)^2} \\
\vec{a}(t) = \left( \frac{d^2x}{dt^2}, \frac{d^2y}{dt^2} \right), \quad | \vec{a}(t) | = \sqrt{ \left( \frac{d^2x}{dt^2} \right)^2 + \left( \frac{d^2y}{dt^2} \right)^2}
\end{gather}
</div>
Using the scaling laws we have:
<div>
\begin{gather}
| \vec{v}(t) | = \frac{1}{\alpha_t} | \vec{v}(\theta) |, \qquad
| \vec{a}(t) | = \frac{1}{\alpha_t^2} | \vec{a}(\theta) |
\end{gather}
</div>
Now we can determine the value of the constant. We know that the magnitude of the acceleration should be equal to $g$:
$$ \| \vec{a}(t) \| = g $$
but the magnitude of the acceleration using the cycloid parametrization is:
$$ \| \vec{a}(\theta) \| = \frac{k^2}{2} $$
Combining these values with the scaling law gives the following value for $\alpha_t$:
$$ \alpha_t = \sqrt{\frac{k^2}{2g}} $$
Now the solutions for $x$ and $y$ can be scaled to the appropriate time horizon with $t_f = \sqrt{\frac{k^2}{2g}} \theta_f$.
---
### Completing the Solution
Now we can determine $v$ using the cycloid equations and the scaling law:
<div>
\begin{align}
v(\theta) &= | \vec{v}(\theta) | = \frac{k^2}{\sqrt{2}} \sqrt{1 - \cos(\theta)} = k^2 \sin\left( \frac{\theta}{2} \right), \quad \theta \in [0,2\pi]
\end{align}
</div>
$$ v(t) = \frac{k^2}{\alpha_t} \sin\left( \frac{t}{2 \alpha_t} \right) = k \sqrt{2 g} \sin\left( \sqrt{\frac{g}{2 k^2}} t \right)$$
We can compute the *control* using the differential equation for $v$:
<div>
\begin{gather}
\dot{v} = g \cos( u(t) ) = g \cos\left( \sqrt{\frac{g}{2k^2}} t \right) \\
u(t) = \sqrt{\frac{g}{2k^2}} t, \qquad u(\theta) = \frac{\theta}{2}
\end{gather}
</div>
For completeness, we have:
$$ x(t) = \frac{1}{2} k^2 \left(t \sqrt{\frac{2g}{k^2}} - \sin\left( t \sqrt{\frac{2g}{k^2}} \right)\right), \qquad y(t) = \frac{1}{2} k^2 \left(1 - \cos\left( t \sqrt{\frac{2g}{k^2}} \right)\right) $$
where we can now solve for $k$ and $t_f$ numerically by using the two final boundary conditions.
---
### Matlab Implementation
The following Matlab code implements this time-based Brachistochrone solution using the [symbolic toolbox](https://www.mathworks.com/products/symbolic.html) for a high-precision solution. We use [symfun](https://www.mathworks.com/help/symbolic/symfun.html) and [vpasolve](https://www.mathworks.com/help/symbolic/vpasolve.html) to solve for $k$ and $\theta_f$. The dynamic optimization problem is also solved using pseudospectral methods at [[link]](https://github.com/danielrherber/basic-multiple-interval-pseudospectral).
```matlab
function [T,X,Y,V,U] = BrachistochroneSolution(G,XF,YF)
% number of time points
Nt = 1000;
% initialize symbolic variables
syms theta k xf yf
% define symbolic functions
eqnX = symfun( k^2/2*(theta - sin(theta)) == xf , xf);
eqnY = symfun( k^2/2*(1 - cos(theta)) == yf , yf );
% solve for thetaf and k
S = vpasolve([eqnX(XF),eqnY(YF)],[theta,k],[0 2*pi; 0 Inf]);
% extract
thetaf = S.theta; k = S.k;
% scaling constant
at = sqrt(k^2/(2*G));
% final time
tf = at*thetaf;
% theta and time horizon
theta = linspace(0,thetaf,Nt)';
T = linspace(0,tf,Nt)';
% states
X = k^2/2*(theta - sin(theta));
Y = k^2/2*(1 - cos(theta));
V = k*sqrt(2*G)*sin(theta/2);
% control
U = theta/2;
% convert to double
X = double(X); Y = double(Y); V = double(V); U = double(U);
end
```
Solution trajectories with $x_f=2$ and $y_f=2$:
[](blogs/matlab/post_8/solution2.png){data-lightbox="blog_imgs" data-title="State and Control Trajectories"}
[](blogs/matlab/post_8/animate2.gif){data-lightbox="blog_imgs" data-title=""}
Solution trajectories with $x_f=10$ and $y_f=2$:
[](blogs/matlab/post_8/solution1.png){data-lightbox="blog_imgs" data-title="State and Control Trajectories" ;}
[](blogs/matlab/post_8/animate1.gif){data-lightbox="blog_imgs" data-title=""}
