Brachistochrone on a Velodrome

TL;DR

This paper extends the classical Brachistochrone problem to find optimal cycling trajectories on velodrome surfaces, combining analytical and numerical methods, and considers effects like fatigue to improve lap time predictions.

Contribution

It introduces a novel approach to determine the fastest cycling path on velodromes by modeling the cyclist as an active particle and solving the problem on real track geometries.

Findings

Analytical solutions for simple geometries like slopes and cones.

Numerical solutions on actual velodrome surface.

Discussion of fatigue effects on optimal trajectories.

Abstract

The Brachistochrone problem, which describes the curve that carries a particle under gravity in a vertical plane from one height to another in the shortest time, is one of the most famous studies in classical physics. There is a similar problem in track cycling, where a cyclist aims to find the trajectory on the curved sloping surface of a velodrome that results in the minimum lap time. In this paper we extend the classical Brachistochrone problem to find the optimum cycling trajectory in a velodrome, treating the cyclist as an active particle. Starting with two canonical cases of cycling on a sloping plane and a cone, where analytical solutions are found, we then solve the problem numerically on the reconstructed surface of the velodrome in Montigny le Bretonneux, France. Finally, we discuss the parameters of the problem and the effects of fatigue.

Equations

$$T={\int }_0^{s_0}\frac{1}{v(s)}\mathrm{d}s,$$

$$\frac{1}{2}m{v}^2+mgz=0,$$

$$T=\frac{1}{\sqrt{2g}}{\int }_0^L\sqrt{\frac{1+(\mathrm{d}z/\mathrm{d}x{)}^2}{-z}}\mathrm{d}x.$$

$$2z\frac{{\mathrm{d}}^{2}z}{\mathrm{d}{x}^2}+{\left(\frac{\mathrm{d}z}{\mathrm{d}x}\right)}^2+1=0,$$

$$z(0)$$

$$z(L)$$

$$\frac{\mathrm{d}z}{\mathrm{d}x}(L)=0.$$

$$x(\lambda )=\frac{L}{\pi }(\lambda -\sin \lambda ),z(\lambda )=-\frac{L}{\pi }(1-\cos \lambda ),$$

$$\frac{1}{2}m\left({\dot{x}}^2+{\dot{Y}}^2\right)+mgY\sin \alpha =0.$$

$$T=\frac{1}{\sqrt{2g\sin \alpha }}{\int }_0^L\sqrt{\frac{1+(\mathrm{d}Y/\mathrm{d}x{)}^2}{-Y}}\mathrm{d}x,$$

$$\frac{1}{2}m{v}^2+mg(R-R_0)\sin \alpha =0,$$

$$T=\frac{1}{\sqrt{2g\sin \alpha }}{\int }_{{\theta }_0}^{{\theta }_1}\sqrt{\frac{R^2+(\mathrm{d}R/\mathrm{d}\theta {)}^2}{R_0-R}}\mathrm{d}\theta ,$$

$$\frac{{\mathrm{d}}^{2}R}{\mathrm{d}{\theta }^2}+\frac{3R-4R_0}{2(R_0-R)R}{\left(\frac{\mathrm{d}R}{\mathrm{d}\theta }\right)}^2+\frac{(R-2R_0)R}{2(R_0-R)}=0,$$

$$R({\theta }_0)$$

$$\frac{\mathrm{d}R}{\mathrm{d}\theta }({\theta }_1)$$

$$m\ddot{x}$$

$$m\ddot{y}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{2}m{v}^2+mgY\sin \alpha \right)=F_{\parallel }v.$$

$$F_p=\frac{2\pi T_{max}}{D}\left(1-\frac{2\pi v}{D{\omega }_{max}}\right).$$

$$v_{eq}=\frac{-4{\pi }^2{T}_{max}+2\sqrt{\pi T_{max}(4{\pi }^3{T}_{max}+\rho C_{d}A{D}^3{\omega }_{max}^2)}}{\rho C_{d}A{D}^2{\omega }_{max}},$$

$$P=P_0+(P_{\max }-P_0)e^{-t/\tau },$$

$$t_{\mathrm{sprint}}=\tau \log \left(\frac{P_{\max }-P_0}{1/2\rho C_{d}A{v}_{\mathrm{peak}}^3-P_0}\right).$$

$$\mathcal{L}=\frac{1}{2}m\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)-mgz.$$

$$\mathcal{L}=\frac{1}{2}m\left({\dot{x}}^2\left(1+{\frac{\partial f}{\partial x}}^2\right)+{\dot{y}}^2\left(1+{\frac{\partial f}{\partial y}}^2\right)+2\frac{\partial f}{\partial x}\frac{\partial f}{\partial y}\dot{x}\dot{y}\right)-mgf.$$

$$m\ddot{x}$$

$$m\ddot{y}$$

$$\frac{1}{\sqrt{1+{\frac{\partial f}{\partial x}}^2+{\frac{\partial f}{\partial y}}^2}}\left(\begin{array}{c}-\frac{\partial f}{\partial x}\\ -\frac{\partial f}{\partial y}\\ 1\end{array}\right)\cdot \left(\begin{array}{c}\ddot{x}\\ \ddot{y}\\ \ddot{z}\end{array}\right)=0,$$

$$\ddot{z}=\ddot{x}\frac{\partial f}{\partial x}+\ddot{y}\frac{\partial f}{\partial y}+\frac{{\partial }^{2}f}{\partial x^2}{\dot{x}}^2+2\frac{{\partial }^{2}f}{\partial x\partial y}\dot{x}\dot{y}+\frac{{\partial }^{2}f}{\partial y^2}{\dot{y}}^2.$$

$$m\ddot{x}$$

$$m\ddot{y}$$

$$\widehat{\dot{\mathbf{x}}}$$

$$\hat{\mathbf{p}}$$

$$m\ddot{x}$$

Full text

{fmtext}

1 Introduction.

In 1696 Johann Bernoulli posed a problem to the scientific community which, after a year and half, had only been solved by a handful of individuals, including Newton and Leibniz [1, 2]. Along with Newton’s minimal resistance problem, it was one of the first mathematical studies that pioneered the field of the variational calculus, and so had an immense influence thereafter [3]. The Brachistochrone problem, whose etymology comes from the ancient Greek for shortest time [1], describes a curve that carries a particle under gravity in a vertical plane from one height to another in minimal time. The solution is equivalent to the path traced out by a rolling circle, also known as a cycloid. Other variations of this problem have included the effects of friction [4, 5, 6], the motion of a disc on a hemisphere [7], the elastic Brachistochrone (elastochrone) [8], and even the quantum Brachistochrone problem [9].

In sports, a natural extension of the Brachistochrone problem is the motion of a cyclist in a velodrome. Track cyclists compete to move around the sloped velodrome surface as fast as possible, much like a higher dimensional Brachistochrone. Whilst there are many different types of velodrome races, the type of race that lends itself most obviously is the qualification time trial [11]. In this case cyclists complete three and a half laps of the velodrome, where only the time for the final $200\,\mathrm{m}$ of the last lap is measured. The cyclists build up speed on a single high gear roughly over the first two laps [12], staying as high up on the velodrome slope as possible to maintain large potential energy. Then, on the final lap they descend the slope and sprint around the track as fast as they possibly can. Similarly to the Brachistochrone, the choice of the descent trajectory is critical. Time is lost if the descent trajectory is too steep or too shallow. The optimal trajectory must balance the exchange between potential and kinetic energy perfectly.

In Fig. 1(a-d) we show four early time snapshots taken from a velodrome qualification race at the velodrome of Montigny le Bretonneux, France. As the cyclist enters the final lap (at an initial speed of around 58 km/hr), they descend into the sharp corner of the velodrome over a period of around 1 s. After the descent, cyclists typically remain at the black line, which is at the bottom of the slope, within around $1\,\mathrm{m}$ of the blue côte d’azur lane (see Fig. 1(d)), for the duration of the final lap. Hence, the initial descent is of critical importance because it is the only period of the race where variation is observed between different cyclists. In Fig. 1(e) we also display the final times of the top five cyclists from the Rio 2016 Olympic games. The difference in time between cyclists is usually of the order of one tenth, and sometimes one hundredth of a second (corresponding to $0.1-1\%$ of the total time). Hence, if the track time can be reduced even slightly by choosing a better descent trajectory, this can have a significant impact on the final ordering of the athletes.

In this study we show how to find the optimum trajectory of a cyclist in a velodrome by modifying the classical formulation of the Brachistochrone problem, treating the cyclist as an active particle on a surface. Since the velodrome track is naturally decomposed into straight and curving sections, we start by studying two canonical cases of motion of a cyclist on a plane and a cone, for which analytical solutions can be found using the Euler-Lagrange equations. Then, using geometrical data taken from the velodrome of Montigny le Bretonneux, we reconstruct the velodrome surface and solve the corresponding optimisation problem using a numerical method that extends from the previous examples. After validating the model by comparison against cyclist velocity and power data, we discuss the various parameters of the problem, as well as the effects of fatigue.

2 Brachistochrone on a plane.

Before discussing motion on two-dimensional surfaces, let us first summarise the classic Brachistochrone problem, formulated in the Euler-Lagrange setting [1, 13].

Consider a particle of mass $m$ that moves in the vertical plane $(x,z)$ under gravity $g$. This is equivalent to motion on a two-dimensional plane in the case where the plane makes an angle $\alpha=\pi/2$ with the horizontal, where $\alpha$ is illustrated in Fig. 2(a). We seek to minimise the total time for the particle to move along a trajectory from position $(0,0)$ to $(L,-H)$, which is given by

[TABLE]

where $s$ is the arclength of the trajectory, varying from [math] to $s_{0}$, and $v=\sqrt{\dot{x}^{2}+\dot{z}^{2}}$ is the speed of the particle. Neglecting friction, the total energy of the particle is conserved, such that

[TABLE]

where we have assumed that the particle is initially at rest. Hence, by using (2) and by rewriting (1) in terms of $x$ and $z$, where $\mathrm{d}s=\mathrm{d}x\sqrt{1+(\mathrm{d}z/\mathrm{d}x)^{2}}$, the total time is

[TABLE]

This quantity (3) can be minimised by solving the Euler-Lagrange equation for the function $z(x)$, which is

[TABLE]

together with the boundary conditions

[TABLE]

In the case where the final height $-H$ is included as an optimisation variable, the boundary condition (6) is replaced by

[TABLE]

It is well-known that the solution is a cycloid which, in the case of the latter boundary condition, is given parametrically by

[TABLE]

for $\lambda\in[0,\pi]$.

To extend the above formulation to motion on a sloping plane is relatively straightforward. Let us now consider a coordinate system $(x,y,z)$, in which the particle moves on a plane $z=y\tan\alpha$ that makes a constant angle $\alpha$ with the horizontal. By considering the rotated coordinate $Y=y/\cos\alpha$ that lies in the plane, it follows that the total energy is given by

[TABLE]

Hence, the total time for the descent is

[TABLE]

which is equivalent to the classic Brachistochrone but with a modified gravity $g^{\prime}=g\sin\alpha$, and a trajectory $Y(x)$ that lies in the plane. Hence, the parametric solution is the same as before (8), except with $z$ replaced by $Y$, and a total descent time multiplied by a factor of $1/\sqrt{\sin\alpha}$.

3 Brachisto-‘cone’

The next most canonical case of a two-dimensional surface is a cone. This is of particular interest to track cycling because the sharp corner at each end of the velodrome is approximately conical, as we discuss later. By converting the above formulation to cylindrical polar coordinates, it is possible to write down the energy equation for motion on the surface of a cone, which is

[TABLE]

where $R=r/\cos\alpha$ is the rotated radial coordinate, and $R_{0}$ is the initial position of the cyclist. Unlike the planar case, where the initial position is irrelevant, in the conical case $R_{0}$ is a necessary parameter, and is related to the curvature at the initial position. The resulting time-minimisation problem is written in terms of the integral

[TABLE]

where $\theta_{1}-\theta_{0}$ is the angle traced out by the trajectory. The resulting Euler-Lagrange equation for the trajectory $R(\theta)$ is

[TABLE]

with boundary conditions

[TABLE]

The Brachisto-‘cone’ problem (13) is different from the Brachistochrone problem (4) since it also takes into account the effects of rotation, such as the centrifugal force. Such forces play an important role in track cycling, because they allow the cyclist to perform tight corners at high velocity by tilting their bike dramatically in the direction of the bend.

Owing to the more complicated form of (13) neither an explicit nor a parameterised closed form solution is available (though one constant of integration can be found by considering the Beltrami identity). However, we can solve the boundary value problem (13)-(15) numerically. In the subsequent sections, we refer to the solutions in the two cases of the planar Brachistochrone and the Brachisto-‘cone’ as analytical, but indeed only the planar case is in closed form.

4 Cycling Brachistochrone

Let us replace the particle in the above examples by a point cyclist (or an active particle). In this case, the cyclist applies a pedal thrust in the direction of motion and experiences aerodynamic drag, such that the total energy is no longer constant. Thus, the dynamics of the cyclist in the planar case, which we derive in Appendix A, are given by

[TABLE]

where $F_{\perp}(t)$ is the leaning force perpendicular to the direction of motion (which does not contribute to the energy), and $F_{\parallel}(t)$ is the force parallel to the direction of motion, which is divided into a pedal thrust and a drag force $F_{\parallel}=F_{p}(t)-F_{d}(t)$. We model the drag using the parameterisation $F_{d}=1/2\rho C_{d}Av^{2}$, where $C_{d}$ is the drag coefficient, $A$ is the combined surface area of the cyclist and the bicycle, and $\rho$ is the density of air [14, 15, 16, 17]. Note, the only contribution from the normal force is the gravity term in (17). The total energy of the system satisfies

[TABLE]

Clearly, if we set $F_{\parallel}=0$ then (18) leads to the former Euler-Lagrange formulation. If $F_{\parallel}\neq 0$ then we cannot solve the problem analytically, but a numerical solution can be found which we discuss shortly. For the case of the Brachisto-‘cone’ with forcing, the dynamics are given by (50)-(51) in Appendix A, and the energy equation is identical to (18), except with $Y$ replaced by $R$. Before discussing the solutions to these cases, we first need a model for the pedalling force, and hence the cyclist physiology.

5 Cyclist physiology

To model the pedalling force $F_{p}(t)$, there are certain mechanical and physical considerations that must be taken into account. In particular, since track cyclists must choose a fixed gear ratio for the duration of the race, the pedalling force depends strongly on the instantaneous pedalling rate, and this relationship depends on the physiology of the individual cyclist.

As shown by Dorel et al. [11], the pedal torque that a cyclist applies in a sprint is a linearly decreasing function of the pedalling frequency. This linear relationship is characterised by two coefficients $T_{max}$ and $\omega_{max}$, which correspond to the maximum possible torque (occurring at zero pedalling frequency) and the maximum possible pedalling frequency (occurring at zero torque). Each cyclist has a sprint performance characterised by these two parameters, and these are easily measured with a pedalling experiment.

The pedalling torque and frequency are related to the pedalling force $F_{p}$ and speed $v$ via the development $D$, which is the distance travelled by one single rotation of the pedals (analogous to the gear). Hence, the pedalling force is given by the linear relationship

[TABLE]

This relationship is similar to the force-velocity equation related to muscle physiology, sometimes called the Hill equation [18]. Note that (19) is only valid for pure anaerobic respiration, and does not include the effects of fatigue. We will discuss the effects of fatigue later in Section 6.

In Fig. 2 we plot solutions to both the planar Brachistochrone and the Brachisto-‘cone’ problem. Light blue dashed curves correspond to the analytical solution in the case of zero pedalling and drag force $F_{\parallel}=0$. Black curves correspond to the numerical solution to the equivalent optimal control problem, which is achieved by formulating an interior point constrained optimisation [19, 20] using the dynamics (16)-(17) to govern the variables $x(t),y(t)$, and using the forcing $F_{\perp}(t)$ as a control function (see Appendix B). The analytical solution is useful for validating the numerical approach, giving us confidence when applying it to the case of non-zero pedalling and drag force, for which an analytical solution is not available. Such solutions (for $F_{\parallel}\neq 0$) are shown on the same plot with dotted blue curves.

To calculate these trajectories we choose values for the model parameters that correspond to realistic cycling scenarios. We take the combined mass of the cyclist and the bike as $m=86\,\mathrm{kg}$, and the physiological characteristics $\omega_{max}=25\,\mathrm{rad/s}$, $T_{max}=260\,\mathrm{Nm}$, correspond to data taken from an elite athlete [11]. The product of the drag coefficient and the surface area is $C_{d}A=0.22\,\mathrm{m^{2}}$, which is equivalent to a streamlined cycling position and modern equipment [21, 22, 15, 16, 23, 17]. Finally, we choose a development of $D=8.5\,\mathrm{m}$, which is typical for velodrome sprint qualification trials.

In Fig. 3, we illustrate how each of these parameters affects cyclist performance. We display the trajectory and velocity profile for the planar cycling Brachistochrone (from Fig. 2(a)) which we use as a base case (though extended to an $80\,\mathrm{m}$ track), as well as three other curves for which we perturb each of the three physiological parameters $T_{max}$, $\omega_{max}$ and $D$. In each case, the perpendicular force profile $F_{\perp}(t)$ and total time are fixed at their base case values. The perturbed values of each of these parameters ($T_{max}=200\,\mathrm{Nm}$, $\omega_{max}=16\,\mathrm{rad/s}$, $D=10\,\mathrm{m}$) result in a reduced velocity and a shallower descent trajectory. In the case of $T_{max}$ and $\omega_{max}$ this can be explained by a reduction in the forward thrust and maximum velocity, respectively. The effect of increasing the development $D$ is less obvious, since this simultaneously increases the maximum velocity whilst decreasing the maximum thrust. However, this becomes more clear by considering the equilibrium velocity

[TABLE]

which is equivalent to a balance between the pedalling force and the drag force. As indicated by Fig. 3(b), $v_{eq}$ is non-monotonic function of $D$ for typical parameter values, and has a maximum at $D\approx 11.5\,\mathrm{m}$. Therefore, whilst increasing $D$ from $8.5\,\mathrm{m}$ to $10\,\mathrm{m}$ raises the equilibrium velocity, the pedal force is simultaneously reduced, resulting in a longer time to achieve equilibrium, and hence a slower, shallower trajectory.

6 From Brachistochrone to velodrome

6.1 Validation

The next step is to apply the above method to find the optimum trajectory on a real velodrome track and validate the model by comparison with cyclist velocity and power data. As a case study, we choose the velodrome of Montigny le Bretonneux in France. For this velodrome, the inside lane is composed of two straight lines of length $L=38\,\mathrm{m}$, connected by two half-ellipses of semi-major and semi-minor axes $a=29.8\,\mathrm{m}$ and $b=24.2\,\mathrm{m}$ (see Fig. 4(a)). Since the velodrome slope varies as one moves around the track, it is convenient to make use of the tangent and normal coordinates measured on the inside lane $(s,n)$ (see the schematic diagram in Fig. 4(c)). In terms of these coordinates, the velodrome surface is written simply as $z=n\tan\alpha(s)$, where the width of the surface is a constant $W=7.9\,\mathrm{m}$. Therefore, unlike the previous examples, motion on the velodrome is bounded, such that $0\leq n\leq W\cos\alpha$. In Fig. 4(b) we plot the angle $\alpha(s)$, which was measured at the race track using an angle-metre. The slope angle is approximately sinusoidal, varying between $14^{\circ}$ in the middle of the straight sections to $45^{\circ}$ in the middle of the curved sections. The reconstructed velodrome surface is shown in (a).

Based on the previous canonical examples, we formulate and solve a numerical optimisation problem for the trajectory, where the dynamics in the straight regions are similar to (16)-(17), and for the elliptical sections of the velodrome we consider a small perturbation from the cone example, where the ellipticity $\epsilon=a/b-1=0.23$ is treated as a small parameter in an an asymptotic expansion (see Appendix A).

In Fig. 5 we display the optimal trajectory, as well as plots of velocity $v$ and power $P=vF_{p}$ as functions of time [24]. For this example, we use the same parameter values as before, except with initial position $s(0)=0,$ $n(0)=W\cos\alpha(0)$, and an initial velocity of $58\,\mathrm{km/hr}$ purely in the $x$-direction, which corresponds to the typical speed before descent for Olympic athletes (further discussion on optimising the initial conditions of the problem is given in the next section). The optimal trajectory is very similar to the strategy employed by athletes, descending before the first corner and then hugging the inside lane thereafter. The descent follows closely to the position of the cyclist in the snapshots in Fig. 1(a-d), which we indicate approximately with white stars. The velocity and power in Fig. 5(b,c) are compared to data taken from an elite cyclist in the velodrome of Montigny le Bretonneux over three different race attempts. The close agreement between our model and the data indicates that the cyclist’s trajectory is already nearly optimal, and illustrates the validity of our model. Furthermore, the time to complete the final $200\,\mathrm{m}$ of the lap is around $9.7\,\mathrm{s}$, which is close to high-ranking Olympic performance, as seen in the table in Fig. 1(e).

6.2 Discussion of the parameters and the effects of fatigue

Although we have so far only discussed optimising the cyclist trajectory, there is also the question of the optimal choice of the initial position and velocity of the cyclist. In this section, we discuss these initial conditions, including how fatigue plays a role in their optimal values.

Written in terms of the tangent and normal coordinates, there are four initial conditions $s(0)$, $n(0)$, $\dot{s}(0)$ and $\dot{n}(0)$. We find that the optimum value of $n(0)$ is always given by $n(0)=W\cos\alpha(s(0))$, which is in accordance with our observation that track cyclists always begin their descent from the top of the velodrome ramp. This choice of $n(0)$ can be explained by the desire to maximise potential energy before descent. Note, that due to the proximity of the velodrome boundary, this also constrains the initial descent angle of the cyclist, such that $\dot{n}(0)=0$. Hence, the initial speed is simply $v(0)=\dot{s}(0)$.

To determine the remaining optimum parameter values, $s(0)$ and $\dot{s}(0)$, one must consider the effects of fatigue. As discussed by [11], the relationship (19) is only valid whilst the cyclist sprints using pure anaerobic respiration and, in reality, this can only last for around $\sim 5\,\mathrm{s}$. Afterwards, the effects of fatigue take action, and eventually respiration becomes aerobic. After this point, the power exerted by the cyclist deviates towards a near-constant plateau of around $P_{0}/m=6.35\,\mathrm{W/kg}$ [25]. The initial velocity is explained by balancing the pedalling force that corresponds to constant aerobic power ($P=P_{0}$) and the drag force, such that $P_{0}/v\sim 1/2\rho C_{d}Av^{2}$, which for a cyclist of mass $86\,\mathrm{kg}$ gives a solution $v=58\,\mathrm{km/hr}$. If the cyclist were to choose a larger initial velocity than this, they would risk eating up some of their energy budget during the final lap. Hence, they choose to start the descent at the maximum possible velocity without triggering the effects of fatigue prematurely.

The final parameter of interest $s(0)$ is one of the main differences between cyclist strategies. Cyclists typically start the descent between $s=3\,\mathrm{m}$ and $s=-11\,\mathrm{m}$ (where $s$ is measured with respect to the finish line). The descent position $s(0)$ is chosen as a tradeoff between speed and fatigue. If the descent starts too late, then it is not possible to build up enough speed for the final lap. On the other hand, if the descent starts too early then, although the final lap will commence at high velocity, the cyclist will quickly become tired.

According to [25], the pedalling power decays exponentially due to fatigue, where the formula

[TABLE]

gives a good fit with cyclist data using parameter values $P_{\mathrm{max}}/m=17.4\,\mathrm{W/kg}$ and $\tau=38\,\mathrm{s}$. The choice of starting position should be such that the cyclist crosses the finish line before fatigue causes any deceleration. Hence, we expect the peak velocity $v_{\mathrm{peak}}$ to be precisely at the finish line. Therefore, at this moment, we expect a perfect equilibrium between the pedalling power (21) and the drag power. From this balance, we derive the sprint timescale

[TABLE]

During the London 2012 Olympics qualification trials, starting positions varied between $s(0)=3\,\mathrm{m}$ and $s(0)=-11\,\mathrm{m}$, with corresponding total sprint times (including the descent) of $t_{\mathrm{sprint}}=13-15\,\mathrm{s}$. According to (22), these sprint times correspond to peak velocity values of around $v_{\mathrm{peak}}=74-75\,\mathrm{km/hr}$, which are close to observations (e.g. $v_{\mathrm{peak}}=77\,\mathrm{km/hr}$ in Figure 5). Hence, this simple scaling argument provides a good explanation for the choice of starting position. To illustrate how the descent position $s(0)$ affects the shape of the descent trajectory, we display different calculated optimum trajectories for various values of $s(0)$ in Appendix C.

7 Concluding remarks

Whilst the classic Brachistochrone problem is limited to one dimension, we have shown how a similar formulation can extend the problem to motion on a plane or a cone, for which the Euler-Lagrange equations yield analytical solutions. By extending the classic Brachistochrone using these two canonical cases, as well as including the effects of cyclist physiology, we have shown how to optimise the trajectory of a cyclist in the real example of the velodrome in Montigny le Bretonneux, finding very close agreement with measured athlete data. In addition to the optimum trajectory, we have also discussed optimising the initial conditions before descent, where the effects of fatigue play a role.

This work not only paves the way for future studies of cycling, but also has implications for research across the field of physics on the topic of active particles on surfaces. For future work, the dynamics of other cyclists in the race could be included [26], where the effects of slipstreaming and racing psychology must be taken into account [27].

\ack

We thank Tullio Traverso and Benjamin Lallemand for useful discussions. We gratefully acknowledge the support from Ecole Polytechnique for the research program Sciences 2024.

Appendix A Derivation of the governing equations

A.1 Cartesian coordinates

Consider the motion of a particle on a two-dimensional surface $z=f(x,y)$ under gravity without any other external forcing. The Lagrangian of the system is given by the kinetic energy minus the potential energy

[TABLE]

Since the motion is confined to the surface, we can eliminate $z$ to give

[TABLE]

The Euler-Lagrange equations for this system can be rearranged to give two second order differential equations for $x(t)$ and $y(t)$, which are

[TABLE]

These equations simplify by noticing firstly that there is no acceleration in the normal direction to the surface

[TABLE]

and secondly that the vertical acceleration is given by

[TABLE]

Hence, (25)-(26) reduce to

[TABLE]

Now consider that the particle is pushed parallel and perpendicular to the direction of motion by forces $\boldsymbol{F}_{\parallel}$ and $\boldsymbol{F}_{\perp}$, respectively. The unit vector in the direction of motion is given by

[TABLE]

and the unit vector perpendicular to the direction of motion, which also lies in the tangent plane of the surface, is given by

[TABLE]

We take the parallel and perpendicular forces as $\boldsymbol{F}_{\parallel}=F_{\parallel}\hat{\dot{\boldsymbol{x}}}$ and $\boldsymbol{F}_{\perp}=F_{\perp}\hat{{\boldsymbol{p}}}$, without loss of generality. Therefore, the full system of dynamical equations becomes

[TABLE]

Note that rate of change of energy is given by

[TABLE]

By inserting (33)-(34) into (35), we find

[TABLE]

Hence, the contributions from the gravitational and perpendicular forces are zero (as expected), whereas the contribution from the parallel force is equivalent to the applied power.

For the planar Brachistochrone we have $\partial f/\partial y=\tan\alpha$ and $\partial f/\partial x=0$. In this case, the governing equations reduce to

[TABLE]

The corresponding energy equation is

[TABLE]

or written in terms of the rotated variable $Y=y/\cos\alpha$, we get

[TABLE]

A.2 Cylindrical polar coordinates

Following the same steps as before, the Lagrangian in cylindrical polar coordinates is

[TABLE]

The resulting Euler-Lagrange equations, after simplification, are

[TABLE]

The unit vector in the direction of motion is given by

[TABLE]

and the unit vector perpendicular to the direction of motion is given by

[TABLE]

Hence, the governing equations are

[TABLE]

The rate of change of energy is given by

[TABLE]

Again, by inserting (46)-(47) into (48) we get

[TABLE]

In the case of a cone, where $\partial f/\partial\theta=0$ and $\partial f/\partial r=\tan\alpha$, (46)-(47) reduce to

[TABLE]

The corresponding energy equation is

[TABLE]

or written in terms of the rotated variable $R=r/\cos\alpha$, we get

[TABLE]

A.3 Calculating surface gradients

As explained in the main text, the velodrome surface is written very simply in terms of the tangent and normal coordinates, such that $z=n\tan\alpha(s)$. Therefore, to solve the dynamical equations above, we need to know how to convert the surface gradients (given in terms of Cartesian or polar coordinates) in terms of $s$ and $n$. In particular, we have the following chain rule identities

[TABLE]

In the straight regions of the velodrome, these calculations are straightforward, since $s$ and $n$ are linearly related to $x$ and $y$, such that

[TABLE]

for some constants $x_{i},y_{i}$, $i=1,3,5$ (using the same region numbering as in the main text).

In the case of the curved sections, a little more thought is needed. First we note that if the curved part of the velodrome were circular, rather than elliptical, $s$ and $n$ would be linearly related to $r$ and $\theta$. Hence, as a simple approach, we consider a weakly elliptical shape (i.e. a perturbation from a circle). For an ellipse with semi-major and semi-minor axes $a(1+\epsilon)$ and $a$, respectively, the curve is given by

[TABLE]

where $\epsilon\ll 1$ is a small parameter. Expanding in powers of $\epsilon$, and keeping only first order terms, we get

[TABLE]

Now, consider the distance from a point with polar coordinates ($r$,$\theta$), to a point on the ellipse with coordinates ($R(\phi)$,$\phi$), for some angle $\phi$, as illustrated in Fig. 6. The distance between these two points is

[TABLE]

or equivalently,

[TABLE]

Considering that for a circle we have $\phi=\theta$, here we write $\phi=\theta+\epsilon\hat{\phi}$. By expanding (62) in powers of $\epsilon$, we find that $d$ is given by

[TABLE]

The normal coordinate $n(r,\theta)$ is given by the minimum possible value of $d$, which corresponds to some value $\hat{\phi}^{*}$. Clearly, it is not necessary to evaluate $\hat{\phi}^{*}$ to calculate $n$ to first order, since we have

[TABLE]

However, we calculate $\hat{\phi}^{*}$ here since it will be necessary when evaluating $s$ later. To find $\hat{\phi}^{*}$ we solve the equation

[TABLE]

which has solution

[TABLE]

At a given angle $\phi$, the arclength $\mathcal{A}$ along the ellipse is given by

[TABLE]

and the tangential distance corresponds to the arclength where $\phi=\theta+\epsilon\hat{\phi}$, such that

[TABLE]

Hence, using (64) and (69), the surface gradients in the elliptical regions (56)-(57) can be calculated directly.

Appendix B Numerical optimisation

In this section we briefly describe our solution method for numerical optimisation. As described in the main text, for the case of non-zero forcing in the direction parallel to the cyclist motion, no analytical solution for the optimal trajectory is available. However, a numerical solution is found by formulating a time-minimisation problem with dynamical system constraints, and following the interior point method [20, 19].

In the full velodrome case, the dynamical equations are given by (33)-(34) in the case of Cartesian coordinates and (46)-(47) in the case of polar coordinates. In the simplified cases of motion on a plane or a cone, the dynamical equations are (37)-(38) and (50)-(51). To formulate the optimisation problem, we first discretise the state variables in time, which are either $x(t)$ and $y(t)$ in the Cartesian regions or $r(t)$ and $\theta(t)$ in the polar regions, as well as the control variable $F_{\perp}(t)$. Note that $F_{p}(t)$ is neither a state variable nor a control variable, since it is given by (19) in the main text. We denote the discretised state variables by $X_{i}$ and the control variable by $F_{i}$. We use a uniform discretisation in time $\delta t$ with $N$ points, such that the total time is given by $T=N\delta t$. The governing differential equations are discretised using a second order forward Euler scheme, producing a system of algebraic equations

[TABLE]

There are $2N$ equations in total since there are always two state variables (either $x$ and $y$ or $r$ and $\theta$). It should be noted that for $j=1,2$, the corresponding equations (70) are precisely the initial conditions

[TABLE]

where $X_{0}$ and $U_{0}$ represent the initial position and velocity of the cyclist.

In the case of the velodrome, as discussed in the main text, the cyclist enters between Cartesian and polar regions. Therefore, it is convenient to split (70) into several partitions corresponding to the dynamical equations in each of these regions. To impose continuity between regions, the initial position and velocity (71)-(72) at the beginning of each region must be imposed to be the same as the position and velocity at the end of the previous one. Hence, instead of solving one optimisation problem, we solve several coupled together via their initial conditions.

Another consideration that must be taken into account is that the cyclist has to remain within the bounds of the velodrome at all times. The latter constraint is written in terms of the normal coordinate $n$, which is given by (58) in the case of Cartesian coordinates, and (64) in the case of polar coordinates, such that

[TABLE]

To enforce the equality constraints (70) we use the quadratic penalty method. This involves placing the equality constraints as a penalty in the objective function to be maximised, which is written as

[TABLE]

The penalty parameter $\mu$ is chosen to be sufficiently large that (70) is imposed accurately, but not too large that the problem becomes ill-conditioned. More details on the choice of $\mu$ are discussed by [19]. To enforce the inequality constraints (73), we use the interior point method. This involves using logarithmic barrier functions, similar to the quadratic penalty in (74). We make use of the IpOpt implementation of the interior point method, as discussed by [20].

The optimisation problem is solved numerically using Newton’s method [19], where gradients are calculated using automatic differentiation in the JuMP package [28] of the Julia programming language [29]. For $N=200$, computation time on a laptop computer is around $10\,\mathrm{s}$, proving a very fast method.

Appendix C Optimum trajectories for different starting positions