We Don't Need no Stinkin' Pitot Tubes

Work in Progress…

So I think we can all agree that (in Australia at least), road bicycles are a disposable middle class trinket, to be cast away like detritus as soon as a new shiny thing comes along. This makes me sad, they are so much fun, good for your health, have a lower impact on the environment than other modes of transport, can be a ludicrously cheap way to get around and they have interesting dynamics. But you know, shiny things and status symbols and all that jazz. You cannot tame the consumerist urge of your average MAMIL. Don’t get me wrong, I like useless technology, I just wouldn’t be caught dead paying a grand for some strain gauges that get called a Power Meter, However, I am also a lazy creature and wouldn’t hack together my own power meter either. So no power meter for me, in fact no road bike for me. I have some skinny tyres on a 10 year old hardtail with mix and match wheels and brakes to get me around the city. But if you like that sort of thing, knock yourself out. Toys are fun. But then there is this. Which seems like the next level. Maybe there are people out there spending money on wind tunnel time so this looks like a cheap option? Who knows, not I. But it got me thinking, surely if you already have a power meter, GNSS measurements and accelerometers like your average gear freak cyclist, I don’t think a pitot tube is going to help you that much. I am going to go through the forces involved and we can see what we can do with what we have already.

Fluid Dynamics is Really, Really Hard

I’m not going to say too much here for fear on embarrassing myself, but at high Reynolds Number the resistive force on a body in a fluid, $F_{drag}$, is reasonably approximated by the following function

$$F_{drag} = 1/2 \rho A c_{d} v_{w}^2$$

where $\rho$ is the density of the fluid, $A$ is the frontal area, $c_{d}$ is the “drag coefficient” and $v$ is the free stream velocity of the fluid (the speed far enough away to be undisturbed by having the thing in the flow). Cool beans. You can imagine that there is this magic coefficient that determines your drag based on your shape, except it sort of is more like a function of speed than a coefficient (ewww). But anyway, let’s pretend. You can see straight away that if you want less force of drag for a fixed speed you can alter either the air density (not practical), your frontal area (very, very practical, that’s what tri-bars are all about) or you can lower your drag coefficient. It is often practical to lump the drag coefficient and frontal area into a “drag area”. Bicycles (and Motorbikes) are an aerodynamic nightmare, they have a terrible drag coefficient. Mostly due to the giant monkey hanging off the back of them. There has been a lot of experimentation with various fairings on motorcycles historically, dustbin fairings had a lot of success in reducing the drag coefficient but were banned due to concerns about stability caused by those interesting dynamics I alluded to earlier, stability in cross-winds and lack of range of motion for the steering. So your best bet for going fast is getting a small frontal area.

Contact Mechanics is Really, Really Hard Also

When two things touch they squish into one another to form a contact patch. If the two materials are elastic then Hertz has you covered. He solved the 2D half space elastic contact problem (along with many other things). He was interested the optical problems caused when lenses were touching due to the elastic deformations in the glass. He casually opened up a whole area of applied mathematics then waltzed away to do other things for which he is more famous. Ok, that is touching solved, what happens with rolling? Say we put a force on one of the bodies in contact, what happens next? Now things get murky. If you are still in 2D (really big radius of curvature in one direction, like a steam roller), still elastic and quasi-static and low tangential force compared to the normal force pushing the bodies together, then you can jump to 1926 and the work of Carter who was interested in stability problems for railway dynamics (a subject close to my heart). Break any of those assumptions and things get really messy. To go to 3D Elastic Contacts, then the most important stuff is probably the work by Kalker. For including non-stationary effects, plasticity, heating, visco-elastic behaviour first consult Johnson and then be surprised that this is still an active area of research! So let’s be properly dumb applied mathematicians (yet aware of our ignorance) and say that for our bicycle, we are going to take a Taylor series for this complicated dynamical phenomenon, just to two terms too!

$$F_{rolling} = a_{stiction} + b_{rolling}v$$

We may as well just chuck all that drive train stuff in there as well, because that is just a bunch of rolling contacts too. So we have a static force to get going and then a slowly rising force that is linear with the rolling speed. If the slip between the wheel and the ground isn’t too big let’s just pretend that this is the bicycle speed.

Good Old High School Physics

What more can I say, this is the old high school physics favourite of a mass on an incline.

$$F_{fatarse} = mgsin(\theta)$$

The Rider

The rider is putting force into the system with their legs. If we have access to a power meter we know the sort of trivial relationship that

$$P = Fv$$

So that we can say

$$F_{rider} = P_{rider}/v$$

Putting it all Together

So if we want to balance all these forces we can pit the driving force against the resistive forces to yield

$$m\dot{v} = P_{rider}/v - \left(1/2 \rho A c_{d}v_{w}^2 + a_{stiction} + b_{rolling}v + mgsin(\theta)\right)$$

Cool, now we have a (crude) 1D model for the longitudinal dynamics of a bicycle. Also, we have all I think we need to estimate the drag area from time series measurements of the longitudinal acceleration $a$, the longitudinal velocity $v$, the rider power $P_{rider}$ and the inclination of the slope $\theta$. You could get the incline angle $\theta$ in many dubious ways, you could get it from the average accelerations deviation from straight down over a time window, you could use the 3D velocity vector from the GNSS and assume the bike is always travelling in the local tangent plane to the ground. Ideally you would do a proper sensor fusion of the inertial measurements with the GNSS (with a vehicle stuck to the ground model) to give you the incline angle, I assume you can get this from your smartphone, if not then get off your arses phone developers! The next section will show how we can put this estimation problem together.

Thomas Bayes to the Rescue

So we ask ourselves yet again, “what is the probability of the various coefficients here (including the drag area) being a certain value, given the data”. Our old friend Bayes Theorem gives us the hint

$$P(coefficients | measurements) \propto P( measurements | coefficients )$$

So if only we had access to a function that relates uncertainty in parameters to these measured values we could find the maximum of this probability to estimate the most likely values for the parameters. Oh wait, we just did that. But we are missing a few bits, what of the wind? And what is the uncertainty in our estimates. First the wind. Let’s just stay in 1D and say that

$$v_{w} = v + v_{wind}$$

i.e. the ground speed plus the (head)wind speed. Now it is probably dumb to assume the wind speed is constant, this defies everyday experience. If we left the wind speed free to be estimated at every time interval we are suddenly less well determined and we will get a batshit crazy answer for the wind speed at every time interval. We could go crazy with prior distributions here, but it is probably enough to just put a smoothness prior on the change in wind speed, let’s say it is a Gauss-Markov random walk thingo. This will have a “differential equation” like

$$\dot{v}_{wind} = -v_{wind}/\tau + N(t)$$

i.e. I have made a dubious proposal that the wind has some correlation time (say 20 seconds I dunno) and it is driven by some noise process. It’s probably fine, you can do better if you want to but it probably won’t make a huge difference to the answer you get. Now is probably a good time to discretize everything in sight. We will get our acceleration, velocity, incline and power output at discrete time intervals. Annoyingly you will get them at different time intervals, but that is a practical detail for clever engineers who want to actually solve this problem, not just fuck around. Even more annoying, you won’t actually measure these values directly, they will be voltages or some digital encoding of the signal. And you won’t even be measuring what you think you are measuring, you will be measuring some change in strain over some time interval in some weird reference frame. Let’s ignore all that stuff. Details are for people who are really trying to do this properly to worry about. So the force balance equation is dead simple (we can just do the finite difference approx to the derivative of velocity), so we just need to discretize the random walk in the wind speed. Because this is completely made up anyway and the signal discretization is smaller than the wind speed correlation time (because I said so), we can do the old Euler approximation to the derivative and saying the driving force is normal of a known variance to get

$$v_{wind_{k+1}} = (1 - \Delta t \tau ) v_{wind_{k}} + N(0, \sigma)$$

There is something messy about Stratonovich Integrals that means there is a $\sqrt(\Delta t)$ scaling on the variance (if you knew the “true” variance and wanted to work with different discretizations), let’s ignore that. Too confusing. All the other measurements will have noise on them as well, which we will asssume is normal as well. Lumping it all together and putting into Bayes Theorem, maximizing the log likelihood etc (exercise left to reader) we will get a regularized (for the wind speed) nonlinear least squares problem.

$$\left\lVert f(c_{d}A, m, a_{stiction}, b_{rolling}, \mathbf{v_{wind}}, \mathbf{v}, \mathbf{a}, \mathbf{P}, \mathbf{\rho}) \right\rVert_{2}^{2} + \lambda \sum_{i=1}^{n-1} (v_{wind_{i+1}} - (1 - \Delta t \tau ) v_{wind_{i}})^2$$

Cool beans. The only gremlin that is immediately apparent is the $P/v$ term, it is likely to be very unreliable when the speed is slow, so let’s not estimate those bits. In the next section I will show how you might solve this in code with a deliberately misleading example.

The Tedious Details

So that was sort of fluffy and hand wavy, if you are anything like me you would have been frustrated by my lack of clear exposition. So let’s be so obvious that we can explain it to a computer. As a starting point we possess regularly sampled time series measurements for the acceleration, $\mathbf{a_{m}}$, the velocity, $\mathbf{v_{m}}$, the incline, $\mathbf{\theta_{m}}$, the rider power, $\mathbf{P_{m}}$ and the air density $\mathbf{\rho_{m}}$ (which you probably get from some altitude, pressure and temperature measurement with a simple model for the atmosphere). The bold font indicates that this is a vector of measurements on a fixed time interval. An individual measurement is indicated with the subscript (e.g. $a_{m_{k}}$). Now we have an option here, we can assume that every measurement has a “true” value, whose distribution we are trying to infer. The problem is that without a dynamic equation to regularize the time series of estimates we will just return the measured value. We can always make fake dynamics (like I did for the wind) but this can be pretty unnatural, especially for things like acceleration. So let’s directly use the measured values for the acceleration, density and power in the 1D dynamics equation. For the speed, $v$ we will make the additional term in the objective function that the measured value should be close to the true value, with normally distributed errors, i.e.

$$v_{error} = \left\lVert \mathbf{v_{m}} - \mathbf{v} \right\rVert_{2}^{2}$$

I remember getting confused seeing people do this sort of simplification when they set up these types of estimation problems. The consistent thing to do would be to treat all the measured variables this way. But you only need to do it if there are dynamics messing with the point estimates coming from measurements. I find this a really interesting way to think of “dynamics” or “physics”, that the equations of motion are really just a regularization term to smooth out the estimate you would get from an individual measurement.

This is getting out of hand, I’m going to park this one until inspiration strikes…