4
$\begingroup$

I am currently building a simplistic video game of bicycle racing, a similar to the idea of Zwift. In this game, the user has a physical bicycle that is connected with various hardware that will let me know the power output, in my case measured at the cranks (using a power2max power meter).

My application receives an instantaneous power reading 4 times a second from this power meter. From this power and a few known constants, I am able to model the velocity of the rider given this power output:

$$ P = \frac{1}{1 - \text{drivetrain loss}} \cdot ( F_{\text{gravity}} + F_{\text{rolling}} + F_{\text{drag}} ) \cdot V $$

An expansion of the classic $P = FV$ formula. In my application, I use Newton's method to solve this for $V$, as I know the other values.

My problem is that just solving for $V$ whenever I receive data from the power meter is incorrect, as that would correspond to instantaneous velocity changes. For example, the rider starts from rest, and if I apply 600W at the pedals, in reality I won't suddenly be travelling at 50km/h - I have to accelerate from rest.

I'm at a loss as to how to correctly calculate this acceleration. My current attempt is a rearrangement of $F=ma$ - $a = \frac{F}{m}$. In this scenario I have $F = \frac{P}{\text{modelled velocity for power}} - \frac{1}{1- \text{drivetrain loss}} ( F_{\text{gravity}} + F_{\text{rolling}} + F_{\text{drag}})$. What I'm doing is calculating the force from the model that is above the current forces the rider is experiencing when they don't apply any force to the pedals. This results in acceleration if the rider is travelling under the modelled velocity for a given power, and decelleration if they are travelling at a velocity greater than the modelled velocity (as they are unable to overcome drag).

While this does seem to converge to the correct values, the calculated acceleration seems extremely low. Returning to the example of 600W, I know that if I go out and do this on the roads I'll be travelling at around 20km/h within a few seconds. Yet under the above if I step the model 4 times a second with 600W (simulating a constant power reading from the power meter), after 5 seconds my velocity would be 3km/h, which is most certainly wrong!

Here is the Haskell code I am using for my work:

import Numeric.AD (auto)
import Control.Monad.Fix (MonadFix)
import Numeric.AD.Rank1.Newton.Double (findZero)

velocityForPower :: Double -> Double -> Double -> Double
velocityForPower grade weight power =
        last (take 10
                   (findZero (\v ->
                                calculatePowerOutput (auto grade)
                                                     (auto weight)
                                                     v -
                                auto power)
                             30))
grade, weight :: Double
grade = 0 
weight = 62

accelerationForPower :: Double -> Double -> Double -> Double -> Double
accelerationForPower grade weight velocity knownPower =
  let forceProvided = knownPower / velocityForPower grade weight knownPower
  in (forceProvided - resistingForces grade weight velocity) / weight

calculatePowerOutput :: (Floating a, Ord a) => a -> a -> a -> a
calculatePowerOutput grade weight velocity =
  resistingForces grade weight velocity * velocity

resistingForces :: (Floating a, Ord a) => a -> a -> a -> a
resistingForces grade weight velocity =
  recip (1 - drivetrainLoss) * (fGravity + fRolling + fDrag)
  where fGravity = g * sin (atan grade) * weight
        fRolling
          | velocity > 1.0e-2 = g * cos (atan grade) * weight * rollingResistanceCoeff
          | otherwise = 0 -- There is no rolling resistance if we are stationary
          where rollingResistanceCoeff = 5.0e-3
        fDrag = 0.5 * dragCoeff * frontalArea * rho * velocity * velocity
          where dragCoeff = 0.63
                frontalArea = 0.509
                rho = 1.226
        g = 9.8067
        drivetrainLoss = 3.0e-2

And here is an example run of the model:

mapM_ (print . (/ 1000) . (* 60) . (* 60))
      (take (round $ 5 * recip 0.25)
            (unfoldr (\v -> Just (v, v + accelerationForPower 0 68.2 v 120 * 0.25)) 0))

0.0

0.20502901826445583
0.3645543623001854
0.5240609569432696
0.683538303824054
0.8429759112371182
1.002363296901935
1.1616899907186995
1.3209455375177925
1.4801194998013418
1.6392014604753573
1.7981810255709199
1.9570478269529188
2.115791525014844
2.274401811358143
2.4328684114546824
2.5911810872908574
2.7493296399919145
2.907303912425068
3.0650937917800167
Improve this question
$\endgroup$
4
  • $\begingroup$ Now that some of the main typos are out of the way we can focus on solving the real problem. If you have 600 W, and a mass of 60 kg and initial velocity of 1 m/s, you should be developing a force of 600 N and an acceleration of 10 m/s$^2$. That would be very hard to do (wheels might be slipping) but agrees with your idea that you should come up to speed very quickly (600 W is really quite a lot of power, but OK for a short burst). I would print out the intermediate values for the different components of friction, power left, and computed force for acceleration at each step. That should help. $\endgroup$
    – Floris
    Jan 1, 2016 at 1:54
  • $\begingroup$ I am not familiar with Haskell but it looks like in calculatePowerOutput you are passing (velocity*velocity) to a function resistingForces that expects velocity (not velocity squared). But it's possible I just don't understand the language. $\endgroup$
    – Floris
    Jan 2, 2016 at 15:39
  • $\begingroup$ @Floris in Haskell function application binds tighter than almost anything else, so should read that line is: calculatePowerOutput = resistingForces (grade, weight, velocity) * velocity $\endgroup$
    – ocharles
    Jan 2, 2016 at 17:54
  • $\begingroup$ As I said - it's not my language. In general I prefer to use parentheses to make sure there is no possibility of getting the order of operations wrong. Unless Haskell doesn't have those... $\endgroup$
    – Floris
    Jan 2, 2016 at 17:55

3 Answers 3

Reset to default
4
$\begingroup$

You initial equation for the power equation seems wrong as you are multiplying things that should not be multiplied. Instead, think of the force needed to keep the bicycle going at the current velocity; multiply that force by the velocity to get the power needed; and finally, see if there is "too much" or "not enough" power to maintain that.

I assume you have expressions for the various forces on the bicycle - rolling resistance, air drag (including wind), and the effect of gravity (slope). Drag will be a strong (quadratic) function of velocity; rolling resistance and the force of gravity much less so.

You have power applied to the pedals, and some efficiency factor $\eta$ of the drive train; this means that the available power $P_{available} = \eta P_{pedal}$. Looking at your code, you seem to have a drive train loss of 0.03 (3%), so $\eta = 0.97$.

At a given velocity, you needed a force $F_{total}=F_{drag}+F_{rolling}+F_{gravity}$, and thus a power $P_{needed}+F_{total}\cdot v$. The "excess" power is the difference between the available power and the power needed, $P_{excess}=P_{available}-P_{needed}$. This means that there is a net force $F_{accel}=P_{excess}/v$, and the instantaneous acceleration of the bike is $a=F_{accel}/m_{total}$

You can then do a simple calculation of the velocity at a later time:

$$v_{t+\Delta t} = v_t + a\Delta t$$

If the acceleration is changing constantly, there are better ways to do the integration; but if you have power values four times per second ($\Delta t = 0.25 s$), you should be OK just using this equation.

I am sorry that I cannot comment on the Haskell code - it's "not my language". So I wrote a little program in Python that does largely the same thing; note that I borrowed the idea first shown by @Shane to compute the new velocity by looking at the energy gained; this has the advantage that it doesn't lead to a singularity when $v=0$ (which would lead to "infinite force"):

# bike power calculation
import math
import numpy as np
import matplotlib.pyplot as plt

# a few plausible factors from http://www.cyclingpowerlab.com/CyclingAerodynamics.aspx
dragCoeff = 0.88
frontalArea = 0.32 

rho = 1.2
eta = 0.97 # 1 - drive train loss = efficiency
rollingCoeff = 5.0e-3 # from haskell code
mass = 62 #kg
grade = 0.0
g = 9.81

# functions to compute the various forces:
def fDrag(velocity):
    return 0.5*dragCoeff*frontalArea*rho*velocity*velocity

def fRolling(grade, mass, velocity):
    if velocity > 0.01:
        return g * math.cos(math.atan(grade)) * mass * rollingCoeff
    else:
        return 0.0

def fGravity(grade, mass):
    return g*math.sin(math.atan(grade))*mass

# the actual program:
v=0.0  # initial velocity
power = 600 # constant power in W
dt = 0.1    # time step
va=[0]      # store the results in a list which will grow in each iteration
ta=[0]

# loop over time:
for t in np.arange(0,10,dt):
    totalForce = fDrag(v) + fRolling(grade, mass, v) + fGravity(grade, mass)
    powerNeeded = totalForce * v / eta
    netPower = power - powerNeeded
    # kinetic energy increases by net energy available for dt
    print "t = %.2f; v=%.1f; drag = %.2f N; F roll = %.2f N; F gravity = %.2f N"%(t, v, fDrag(v), fRolling(grade, mass, v), fGravity(grade, mass))
    v = math.sqrt(v*v + 2 * netPower * dt * eta / mass)
    va.append(v)
    ta.append(t+dt)

plt.figure()
plt.plot(ta, va)
plt.xlabel('time (s)')
plt.ylabel('velocity (m/s)')
plt.show()

This gives the following plot:

enter image description here

The values look fairly plausible. Note that drag (23.25) plus rolling resistance (3.04 N) add up to a little over 27 N by the end - accounting for about half the power used. That makes sense - at 26 mph, you need on the order of 300 W to keep going.

I recommend that you check intermediate results from your program against some of these values to see where the difference is coming from.

UPDATE

Here are the values produced by my code (I ran time steps of 0.1 second):

t = 0.00; v=0.0; drag = 0.00 N; F roll = 0.00 N; F gravity = 0.00 N
t = 0.10; v=1.4; drag = 0.32 N; F roll = 3.04 N; F gravity = 0.00 N
t = 0.20; v=1.9; drag = 0.63 N; F roll = 3.04 N; F gravity = 0.00 N
t = 0.30; v=2.4; drag = 0.95 N; F roll = 3.04 N; F gravity = 0.00 N
t = 0.40; v=2.7; drag = 1.26 N; F roll = 3.04 N; F gravity = 0.00 N
t = 0.50; v=3.0; drag = 1.57 N; F roll = 3.04 N; F gravity = 0.00 N
t = 0.60; v=3.3; drag = 1.88 N; F roll = 3.04 N; F gravity = 0.00 N
t = 0.70; v=3.6; drag = 2.19 N; F roll = 3.04 N; F gravity = 0.00 N
t = 0.80; v=3.8; drag = 2.49 N; F roll = 3.04 N; F gravity = 0.00 N
t = 0.90; v=4.1; drag = 2.80 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.00; v=4.3; drag = 3.10 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.10; v=4.5; drag = 3.41 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.20; v=4.7; drag = 3.71 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.30; v=4.9; drag = 4.01 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.40; v=5.0; drag = 4.31 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.50; v=5.2; drag = 4.60 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.60; v=5.4; drag = 4.90 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.70; v=5.5; drag = 5.19 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.80; v=5.7; drag = 5.48 N; F roll = 3.04 N; F gravity = 0.00 N
t = 1.90; v=5.8; drag = 5.78 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.00; v=6.0; drag = 6.06 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.10; v=6.1; drag = 6.35 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.20; v=6.3; drag = 6.64 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.30; v=6.4; drag = 6.92 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.40; v=6.5; drag = 7.20 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.50; v=6.7; drag = 7.48 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.60; v=6.8; drag = 7.76 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.70; v=6.9; drag = 8.04 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.80; v=7.0; drag = 8.32 N; F roll = 3.04 N; F gravity = 0.00 N
t = 2.90; v=7.1; drag = 8.59 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.00; v=7.2; drag = 8.86 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.10; v=7.4; drag = 9.13 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.20; v=7.5; drag = 9.40 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.30; v=7.6; drag = 9.67 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.40; v=7.7; drag = 9.93 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.50; v=7.8; drag = 10.20 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.60; v=7.9; drag = 10.46 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.70; v=8.0; drag = 10.72 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.80; v=8.1; drag = 10.97 N; F roll = 3.04 N; F gravity = 0.00 N
t = 3.90; v=8.2; drag = 11.23 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.00; v=8.2; drag = 11.48 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.10; v=8.3; drag = 11.73 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.20; v=8.4; drag = 11.98 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.30; v=8.5; drag = 12.23 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.40; v=8.6; drag = 12.48 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.50; v=8.7; drag = 12.72 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.60; v=8.8; drag = 12.97 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.70; v=8.8; drag = 13.21 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.80; v=8.9; drag = 13.45 N; F roll = 3.04 N; F gravity = 0.00 N
t = 4.90; v=9.0; drag = 13.68 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.00; v=9.1; drag = 13.92 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.10; v=9.2; drag = 14.15 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.20; v=9.2; drag = 14.38 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.30; v=9.3; drag = 14.61 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.40; v=9.4; drag = 14.84 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.50; v=9.4; drag = 15.07 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.60; v=9.5; drag = 15.29 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.70; v=9.6; drag = 15.51 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.80; v=9.6; drag = 15.73 N; F roll = 3.04 N; F gravity = 0.00 N
t = 5.90; v=9.7; drag = 15.95 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.00; v=9.8; drag = 16.17 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.10; v=9.8; drag = 16.38 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.20; v=9.9; drag = 16.60 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.30; v=10.0; drag = 16.81 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.40; v=10.0; drag = 17.02 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.50; v=10.1; drag = 17.22 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.60; v=10.2; drag = 17.43 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.70; v=10.2; drag = 17.63 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.80; v=10.3; drag = 17.84 N; F roll = 3.04 N; F gravity = 0.00 N
t = 6.90; v=10.3; drag = 18.04 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.00; v=10.4; drag = 18.23 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.10; v=10.4; drag = 18.43 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.20; v=10.5; drag = 18.63 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.30; v=10.6; drag = 18.82 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.40; v=10.6; drag = 19.01 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.50; v=10.7; drag = 19.20 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.60; v=10.7; drag = 19.39 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.70; v=10.8; drag = 19.57 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.80; v=10.8; drag = 19.76 N; F roll = 3.04 N; F gravity = 0.00 N
t = 7.90; v=10.9; drag = 19.94 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.00; v=10.9; drag = 20.12 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.10; v=11.0; drag = 20.30 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.20; v=11.0; drag = 20.48 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.30; v=11.1; drag = 20.66 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.40; v=11.1; drag = 20.83 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.50; v=11.1; drag = 21.00 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.60; v=11.2; drag = 21.17 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.70; v=11.2; drag = 21.34 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.80; v=11.3; drag = 21.51 N; F roll = 3.04 N; F gravity = 0.00 N
t = 8.90; v=11.3; drag = 21.68 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.00; v=11.4; drag = 21.84 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.10; v=11.4; drag = 22.01 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.20; v=11.5; drag = 22.17 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.30; v=11.5; drag = 22.33 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.40; v=11.5; drag = 22.49 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.50; v=11.6; drag = 22.64 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.60; v=11.6; drag = 22.80 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.70; v=11.7; drag = 22.95 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.80; v=11.7; drag = 23.10 N; F roll = 3.04 N; F gravity = 0.00 N
t = 9.90; v=11.7; drag = 23.25 N; F roll = 3.04 N; F gravity = 0.00 N
Improve this answer
$\endgroup$
8
  • $\begingroup$ Apologies, those multiplications are a typo in my question and not present in the Haskell code, I'll correct for that. $\endgroup$
    – ocharles
    Jan 1, 2016 at 1:33
  • $\begingroup$ Still not sure where your $\rm{\frac{1}{drivetrain~loss}}$ comes from when you have loss=0.03. I think you need $\rm{\frac{1}{1-(drivetrain~loss)}}$ $\endgroup$
    – Floris
    Jan 1, 2016 at 1:40
  • $\begingroup$ You are correct. Poor LaTeXing from me. $\endgroup$
    – ocharles
    Jan 1, 2016 at 1:45
  • $\begingroup$ This still feels a little to slow in terms of acceleration to me, but it's the most convincing answer so I've accepted it. Thanks for helping! I'll continue to explore modelling this to see if I can work out why in reality 600W accelerates me much faster. (I know that I can hold 600W for about 5 seconds without running out of energy, and after that I'm well up to 30km/h) $\endgroup$
    – ocharles
    Jan 2, 2016 at 11:05
  • $\begingroup$ Note that my plot has units of m/s - so after five seconds you are at roughly 9 m/s which is (9x3.6=)32 km/h; this agrees nicely with your experience. $\endgroup$
    – Floris
    Jan 2, 2016 at 14:10
4
$\begingroup$

You should first focus on energy, then find the kinetic energy at different steps in time and then you know the velocity. From there you can calculate the acceleration if you want.

So at one time $t_1$ you have kinetic energy $1/2mv_1^2$. After one time step this energy changes by $P\Delta t\eta$, where $\eta $ is efficacy.

So your new kinetic energy $1/2mv_2^2 = 1/2mv_1^2 + P\Delta t\eta$

Solving for v: $v_2 = \sqrt{v_1^2 + \frac{2P\Delta t\eta}{m}}$

Acceleration can now be computed but if you just wanted velocity is not necessary.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Something seems to be missing here though, because that would imply that if I provide a fixed power my velocity increases indefinitely - what's missing to bound the velocity? $\endgroup$
    – ocharles
    Dec 31, 2015 at 17:12
  • $\begingroup$ Viscous resistance from the air! $\endgroup$
    – bartavelle
    Dec 31, 2015 at 17:40
  • $\begingroup$ I don't think your solving for $v_2$ correctly either. I think $v_2 = \frac{\sqrt{2 \delta t P + m v_1^2}}{\sqrt{m}}$ $\endgroup$
    – ocharles
    Dec 31, 2015 at 18:56
  • $\begingroup$ Sorry ignore my last comment, I read the equation incorrectly. $\endgroup$
    – ocharles
    Jan 1, 2016 at 11:16
  • $\begingroup$ @ocharles I reformatted the equation to show the calculation (which I borrowed for my demo code) is correct. $\endgroup$
    – Floris
    Jan 2, 2016 at 14:16
2
$\begingroup$

Here is my attempt:

import Data.Foldable

main = do
  let vs = iterate step v0 -- list of velocities every 0.1 s
  forM_ (take 1000 $ zip [dt, (2 * dt) ..] vs) $ \(t, v) -> do
    putStrLn $ show t ++ "\t" ++ show v
  where
  v0 = 0
  step v | v < v' = v + f0 / m * dt - c * rho * a * v * v / m * dt
  step v = v + p / m / v * dt - c * rho * a * v * v / m * dt
  c = 0.5
  m = 70 -- kg
  a = 0.33 -- m^2
  dt = 0.1 -- s
  p = 400 -- W
  rho = 1 -- kg/m^3
  f0 = p / v'
  v' = 7 -- m/s

-- to run:
-- runhaskell main.hs | graph -X "Time, sec" -Y "Velocity, m/sec" | plot -TX

It is based on "Computational physics" book by Nicholas Giordano.

The idea is to use apply constant force f0 initially and switch to constant power p at the moment when the force produces the same power, e.g. f0 * v' = p. The critical velocity v' is chosen semi-arbitrary.

Results:

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.