dvmarinoff / Flux

Indoor Cycling App for Structured Training
https://flux-web.vercel.app
GNU Affero General Public License v3.0
506 stars 86 forks source link

Ability to create workouts off Distance #107

Closed PhatWheZ closed 8 months ago

PhatWheZ commented 2 years ago

Putting together a workout file based on estimated durations per data point for a upcoming event. Would it be possible to generate a workout_file using Distance rather than Duration?

`<workout_file>
    <author>PhatWheZ</author>
    <name>Karoo2Coast</name>
    <category>Actual Route</category>
    <description> Karoo2Coast route based on Data points provided by organisers</description>
    <sportType>bike</sportType>
    <tags>
    </tags>
    <workout>
        <Warmup Duration="120" PowerLow="0.32" PowerHigh="0.39"/>
        <FreeRide Duration="300" Sim="0.0"/>
        <FreeRide Duration="2880" Sim="6.5"/>
        <FreeRide Duration="1380" Sim="-1.9"/>
        <FreeRide Duration="1260" Sim="2.9"/>
        <FreeRide Duration="2520" Sim="-5.1"/>
        <FreeRide Duration="1560" Sim="4.8"/>
        <FreeRide Duration="1320" Sim="-3.2"/>
        <FreeRide Duration="2340" Sim="9.0"/>
        <FreeRide Duration="600" Sim="0.0"/>
        <FreeRide Duration="1320" Sim="2.6"/>
        <FreeRide Duration="1380" Sim="-3.1"/>
        <FreeRide Duration="1140" Sim="0.6"/>
        <FreeRide Duration="480" Sim="6.3"/>
        <FreeRide Duration="1320" Sim="-4.2"/>
        <FreeRide Duration="1200" Sim="4.6"/>
        <FreeRide Duration="600" Sim="-4.8"/>
    </workout>
</workout_file>

image

PhatWheZ commented 2 years ago

A program I have come across is ridewithgps.com which allows me to export the route as a .gpx file.

I have ready the threads of the "inaccuracy of conversion" so I'm not too worried about the specific experience, but rather the general overview.

Will definitely see what i can find to convert on the opensource / free side of things and give that a bash

dvmarinoff commented 2 years ago

@PhatWheZ using Distance instead of Duration is not supported right now, but it's something I've been wanting to add. I'll come up with a couple of options and post here.

The main obstacles to tackle:

PhatWheZ commented 2 years ago

@dvmarinoff , keen to see what you can come up with. Was able to get the .gpx converted into a .zwo file, however the duration of the workout is almost 1 hour quicker, so highly unlikely ill survive that :)

mjunkmyjunk commented 2 years ago
  • calculating speed from power

You can get some reference here https://github.com/WouterJD/FortiusANT/blob/31c210799c2914f22e187cb389d42a880a34b368/pythoncode/usbTrainer.py#L752

https://github.com/GoldenCheetah/GoldenCheetah/blob/cabe078453cbe8e136b8adf8d4187a5be878671b/src/Train/BicycleSim.h

GCuser99 commented 2 years ago

@dvmarinoff, I like the slope-mode with respect to distance idea as well. I remember training for mountain bike races and turning my activity gpx (with power meter) into an erg representation of the race course but it sucked for training. I realized then that erg mode is important and excellent for interval training but not much else. I think being able to ride a gpx course (or slope-distance simplification thereof) would be a nice feature. Oh - and then we'll ask you to output virtual elevation in the fit file as well so mountain bikers like me can get their Strava elevation stats!

On the physics behind speed calculation - one of the most commonly referenced papers is "Validation of a Mathematical Model for Road Cycling Power, Martin et al". I can provide pseudo code... I'm not 100% sure as I only glanced at the two links by @mjunkmyjunk but it didn't appear that the Gribble one included the KE term - including KE should be given strong consideration...

The link based on Gribble uses a numerical solution to pulling speed out of the power() equation. However, speed can be solved directly using one of many cubic root solvers. Depending on the desired accuracy, one method may be significantly faster than the other. It also uses the small angle approximation for slope - I'm not sure what the cost is for calculating sin(atan(slope)) and cos(atan(slope) in your implementation but that maybe should be evaluated too.

The link to Golden Cheetah is interesting - I think it's based on http://www.kreuzotter.de/english/espeed.htm#forml which I've seen cited on cycling forums. The power relation includes a "dynamic rolling resistance" term which I have yet to find an associated paper reference for (?). It might complicate usage as the total or "effective" rolling resistance is broken into static and dynamic components, of which user or calling program specifies only the static portion. I'm not asserting that there is anything wrong with this, but it might be confusing as most other simulators don't model a dynamic component, so user or calling program specifies the total rolling resistance (Crr), which is static. I've done comparisons between this model and the one from Martin et al - one has to specify a lower Crr for the former (e.g. Crr=.0033) versus a higher Crr for the later (e.g.=.004) to calculate the same (approximate) power. Sufficiently confused? Read very bottom topic here.

I have verified through model fitting of the output speed and power, that my Wahoo Kickr 2016 (Ble) uses a Sim mode physics model that is closely in line with Martin et al and includes the KE term. It does not include the dynamic rolling resistance term.

Thanks again!

mjunkmyjunk commented 2 years ago

https://github.com/GoldenCheetah/GoldenCheetah/blob/ec6d7838d158816dfcece3a1e33b0831cf4d779c/src/Train/BicycleSim.cpp#L190

@GCuser99 I see reference to the KE term

the first link also has a 2nd calc. https://github.com/WouterJD/FortiusANT/blob/31c210799c2914f22e187cb389d42a880a34b368/pythoncode/usbTrainer.py#L777

I don't know how to read the link as it's not in English

GCuser99 commented 2 years ago

@dvmarinoff in case it's of any use, attached is an MS Excel speed-power calculator (zipped) with VBA code based on Martin et al. There are two sheets - one to verify my VBA code using example calculations given in the paper. The other sheet is a general calculator designed to let you test sensitivity of some of the various perturbations of the basic speed-power model (such as how big of a difference does the small angle approximation make in calculating the P_grav and P_roll components).

PS: I apologize in advance for my poor programming practices - I am a scientist by background with minimal knowledge of software design!

power test model.zip

dvmarinoff commented 2 years ago

@GCuser99 Brilliant!

will check these out, hope Open Office has support for VBA.

dvmarinoff commented 2 years ago

It opens with Open Office! Will try to port it to javascript.

GCuser99 commented 2 years ago

After having some time in surgery recovery to further reflect about this, here are my thoughts about some of the issues...

I think it would be nice to have an implementation that will yield reasonable speed-distance (and elevation?) results for both Erg and Sim modes.

There is the decision about what power-speed model to use. For example, do you include:

My personal preference is to keep the model on the simple side and be able to back it up with research reference(s).

I think that the model should account for Pke in Sim mode.

And then there are the implementation details...

I realized that the method for calculating speed that I provided above in this thread using the root solver is not useful in real-time as you don't know acceleration at the current state that you are trying to solve for velocity. In real-time what you have is power measured at the current state, velocity calculated at previous state, and delta-time between states. So needing acceleration as a direct input is an obvious no-go. Sorry that I didn't think of this earlier. I was using my code to fit the output of my trainer power and speed, in order to verify the underlying model and parameters being used, and so it was fine for that purpose.

I think a workable implementation is something along the lines of this here: https://physics.stackexchange.com/questions/226854/how-can-i-model-the-acceleration-velocity-of-a-bicycle-knowing-only-the-power-ou/226856.

To summarize that solution - given the previous calculated velocity, you calculate the power needed to sustain that velocity (Pzero_accel). Then calculate the difference between the current power and Pzero_accel to get the Kinetic Energy component (Pke) at current state. A non-zero Pke implies an acceleration and thus a change in velocity, everything else being held equal*. From the calculated Pke and delta time one can then back out the current state's velocity.

It appears that GC implementation in BicycleSim by Eric Christoffersen mentioned above takes this approach. It also appears to me that it is being used for both Erg and Sim modes (?).

I will do some more digging and testing and report back later....

Cheers...

*What happens if not all is held equal, such as change in slope between states?

dvmarinoff commented 2 years ago

@GCuser99 Best wishes and speedy recovery! Hope everything is well now!

dvmarinoff commented 2 years ago

About the model, ideally I want to have speed represented as an accumulation of accelerations over time. Something in the lines of this article and the physics exchange answer you linked. I am gonna code something today and see how far I get.

dvmarinoff commented 2 years ago

Started work on this feature. Code is going to be on the distance branch and is hosted here: https://flux-test.vercel.app/. I've added a new Virtual Speed field in the place (just for now) of the Distance field next to the regular Speed field.

The model is in here: https://github.com/dvmarinoff/Flux/blob/distance/src/physics.js,

Where I have included a direct translation to javascript of the VBA calcVel and calcPower functions, a Cubic solver function I quickly grabbed from stackoverflow, and the Model object function, which has a refactored version of calcVel called powerToMaxSpeed, and a virtualSpeed function. It's just a first attempt, needs a lot of refining.

UPDATE: added some fixes to slope

GCuser99 commented 2 years ago

Great start on this @dvmarinoff. I think there is an error in the following line in virtualSpeed:

const force = (power * (1 - drivetrainLoss)) - (totalResistance * speed);

The result of what is on the right of the assignment is power (force*speed). Below is one possible refactor to get speed and acceleration (please run through debugger as I have not!). But note that the speed and acceleration are approximate in that they are discrete (not instantaneous) and lag the current time. To reduce error, dt's will need to be small (hopefully a lot less than 1s?) for distance integration and sampling.

 function virtualSpeed(args = {}) {
        const powerMeasured  = args.power; // W
        const slope          = args.slope ?? defaults.slope; // 0.01 is 1%
        const mass           = args.mass ?? defaults.mass; // kg
        const windSpeed      = args.windSpeed ?? 0; // m/s
        const draftingFactor = args.draftingFactor ?? defaults.draftingFactor; // 0..1
        const dt             = args.dt ?? 1; // s

        let speedPrev        = args.speed ?? 0; // m/s

    // let acceleration     = args.acceleration ?? 0;

    // calculate power needed to sustain previous speed
        // so we set acceleration = 0 (no change in speed)
        // since power measured at fly wheel, dtloss is zero

        let gravitationalResistance = g * mass * sinBeta;
        let rollingResistance = g * mass * cosBeta * crr;
        let windResistance = 0.5 * rho * (CdA + spokeDrag) * Math.pow((speedPrev + windSpeed), 2) * draftingFactor;
    let keResistance = 0; // zero acceleration

        let totalResistance = gravitationalResistance + rollingResistance + windResistance + keResistance;
        let powerSteadyState = totalResistance * speedPrev;  // zero-acceleration power

        // calculate pke in this time step 
        let powerKE = powerMeasured - powerSteadyState;

    // if pke is non-zero, then this implies a change in speed - solve for it
    // note that speed and acceleration below are discrete (NOT instantaneous) and lag in time
    // so small dt's will be needed to reduce the error in integration and sampling for distance
        // don't forget to include wheelInertia mass since we are using pke to back out speed

        let speedNow = Math.sqrt(Math.pow(speedPrev, 2) + 2 * powerKE * dt / (mass + wheelInertia));

    if(speedNow < 0) speedNow  = 0;

        let acceleration = (speedNow - speedPrev)/dt;

        return {acceleration, speedNow };
    }

Below is the result of a quick experiment to check if the behavior of above method looks reasonable (performed in Excel). You can see the acceleration of speed at the beginning, and deacceleration when the power changes from 150 to 50. The magnitude of the numbers appears to be reasonable as well. But this was a very idealized experiment where the time sampling is very fine and at a constant sample rate. I will do some more realistic simulations next - will compare notes as we go....

Also, just as an interesting observation, of potentially little or no use - if one was only interested in Erg mode distance, then speed and distance, given rider weight, could be pre-computed before the workout ever starts because you already know the power output profile!

dvmarinoff commented 2 years ago

@GCuser99 I can't get this version of the model to work will need help. Created a notebook-like executable version of the physics.js file with observablehq here: https://observablehq.com/d/17cf1ccfcc8cfa56, maybe we can collaborate over there to get this working.

Main question I think is how to update keResistance when acceleration in not 0?

GCuser99 commented 2 years ago

Ok no problem. First a couple of questions:

1) should I use the defaults in function Model if argument not specified in Function test? 2) did you mean to set dynamicCrr=true and spokeDrag=true in Function test?

In virtualSpeed we never update keResistance - it's always zero because we are calculating steadystate power, ie, power needed to sustain last state's speed. Acceleration is an output from virtualSpeed, not an input.

One potential problem is that when I refactored virtualSpeed I assumed from your original that you had decided not to include the modeling options other than wheelInertia. So I left those out of my refactor... I see now that because you are comparing results to powerToMaxSpeed that maybe my assumption was wrong. I will fix that.

Your velocity calculation mismatch is likely due in part to the issue that velocity and acceleration are discrete in your implementation (as they must be), but in powerToMaxSpeed the inputs need to be instantaneous (which you don't have). The results should converge as the sample rate gets very small, but I doubt 1s is small enough.

Anyway, I don't think the discrete vs instantaneous issue is a problem with virtualSpeed - it's a problem with powerToMaxSpeed.

So another question - how many times a second can you get power measurements from the trainer?

Anyway, I think you are closer than you think!

dvmarinoff commented 2 years ago

I just factored out of the main function some of the parameters, so they can be reused between functions as some of them doesn't exactly need to update in real time (weight, etc.). The defaults are just fallback values as specified in the ANT+ spec, you can use any you wish and the use object in defaults is used to switch on and off the various boolean parameters, so dynamicCrr and spokeDrag are set to true for this test, but they can be set anyway.

The previous naive version of the code in the worst case is about 0.44 km/h out of the gribble and kreuzotter calculators.

Most trainers will transmit power at 4 times per second. But it is not exact it will average 4 times per second over the course of a minute or so. I can sample speed at any rate though.

dvmarinoff commented 2 years ago

I added grade and altitude to the .fit file. Grade is recorded successfully, but altitude is set to 0. I will need to think on what would be the most efficient way to calculate the updates to altitude and where to fit them in the code.

GCuser99 commented 2 years ago

Ok the first order problem is my very bad vba to js translation! I apologize. In virtualSpeed change to the following:

speed = Math.sqrt(Math.pow(speedPrev, 2) + 2 * powerKE * dt / (mass + wheelInertia));

Now we are a lot closer. I will work on the other problems while you are sleeping :-)

GCuser99 commented 2 years ago

I'm finished with this pass.

Aside from the change I mentioned above (which BTW was correct in my original version of virtualSpeed above so I don't know what happened!), here are some other issues I fixed:

With the following settings:

dtloss=0,             
spokeDrag: true,
bearingLoss: true,
wheelInertia: true,
dynamicCrr: false,
smallAngleApprox: false,

I get results

{speedReached: 28.42939689674408, speedPredicted: 28.288156916893566, error: 0.14123997985051417}

If I set dt=.1s and run through 300 time steps I get:

{speedReached: 28.30434147581183, speedPredicted: 28.290233113804895, error: 0.014108362006936659}

Some general comments:

Here is the modified code: physics.txt

Thanks for letting me help - cheers!

dvmarinoff commented 2 years ago

@GCuser99 Thanks for the help! It was really a set of rogue parens that caused the issue. The new models is giving very very small error, however it takes some time to decelerate to 0 km/h which is about 120 seconds independent of the initial speed or the time resolution (dt). Will do more experiments and report later. Is it worth the time to dive in the other integrators, it appears that Golden Cheetah is using something called KahanLi8 from Julia (probably the language)?

GCuser99 commented 2 years ago

Using Martin et al modeling options (true, true ,true ,false ,false), crr=.004, CdA=.4, Rho=1.275, mass=85 kg, v0=7 m/s, slope=0, I get 99.2 secs and 255 m stopping time and distance (which does seem long, intuitively). But in the last 20 secs of that, I only traveled 7 m (practically track standing)!

I did see a dependence on v0. But as v0 gets higher and higher, the increment in stopping time/distance gets smaller. I'm thinking that is due to wind resistance being proportional to v^3, and thus the higher v0 the more energy gets zapped out of the system causing a dampening effect.

Unfortunately, I can't get on a bike right now, otherwise I would run a few tests outdoors and in Zwift. I will look back at some virtual rides I have done in the past and check the model against those and report back...

I haven't yet looked at integration and sampling, but will do that as well when I get some time.

Cheers!

dvmarinoff commented 2 years ago

No worries, the model is pretty good right now. Worst case we implement 'virtual brakes' or something.

Cheers!

GCuser99 commented 2 years ago

@dvmarinoff - on the integrators issue - can you explain how you are integrating/sampling now? From what you said above I understand that you get power-time messages from the trainer at a non-constant sample rate of about .25 seconds plus or minus. Does the trainer also output velocity and distance too (even if meaningless)? How do you sample the data stream onto a 1s grid for the fit file?

dvmarinoff commented 2 years ago

Over Bluetooth I get a series of discrete messages from the device. The indoorBikeData message can contain the latest measurement for power, speed, cadence, distance and heart rate if those are supported by the trainer. Most trainers output power and speed, very few do heart rate (by proxying a heart rate monitor over your trainer), and I haven’t seen yet a trainer that transmits distance. Speed is almost always a function of flywheel speed and is suppose to map closely to speed if riding on flat and doesn’t take into account simulation parameters like slope and wind speed (Wahoo trainers are an exception here).

I accumulate the incoming data, and run a loop every second that takes an avg of the accumulation and reports it as power1s or 1 second power, this is what goes into the .FIT file. There is another loop that collects 30 second power, that I planning to use for a live TSS estimation. Also there are other ‘meta-properties’ like powerLap and powerAvg, which accumulate until the corresponding lap or workout event occurs.

At the moment I am using the real time power stream as source for virtualSpeed, but I can easily construct a loop that samples whatever is accumulated if any for a specific dt and not rely on the device transmission rate.

All of this is implemented through a proxy meta object called db that dispatches notification events when there is a write on any of its properties. This object serves as a sink for data incoming from various connected devices. MetaProperties use this data to construct complex data like one second power and lap average power.

GCuser99 commented 2 years ago

@dvmarinoff , thanks for that explanation - very useful. I think using the real time power stream as a source for virtualSpeed is the correct approach, and then resample for .FIT on a separate loop. That way you are getting the finest sample rate possible to reduce error.

Quite frankly I am having difficulty with why the choice of integrator is so important in this application. Firstly, I don't understand the trainer mechanics of calculating power. For example does the trainer accumulate work over a given time interval and output P=dW/dt? In that case power would represent a lagging average interval power (block average) and not instantaneous power. I would think that knowing the answer to that would drive the integration decision. Also, if you can use the real time power stream for calculations, the errors in integration will be minimized due to the fine sample rate. I have a feeling that this is not going to be a big issue, and that for now dd=speednow*dt is fine. I will do further testing to verify this though.

On the virtualSpeed calcs - I've done some further checking and have made the following observations:

Although GC's model is different (eg. it includes dynamic rolling resistance), the technique for calculating speed looks equivalent to the Flux method. For example, the important Flux code line:

speed = Math.sqrt(Math.pow(speedPrev, 2) + 2 * powerKE * dt / (mass + wheelInertia));

Is precisely equivalent to this in GC:

speed = VFromKE(KEFromV(speedPrev, mass + wheelInertia) + powerKE * dt, mass + wheelInertia);

Also, if you want to add a paper reference for that code line, it could be Martin et al:

Pke=.5*(m+I/r)*(vf^2-vi^2)/dt

If one solves the above for vf (speed final), one gets the same that you are using.

I can't yet get on a bike or trainer, so instead I took a previous Flux Sim mode test and calculated virtualSpeed for comparison to trainer speed. I used all the same inputs, with Martin modeling options. I had to shift forward the slope by 2 secs to line up changes in velocity to changes in slope. Note that the overall match is pretty good, except on the first acceleration in the beginning. The overall distance calculated matches pretty good as well, within a half of a percent, using the simple integrator dd=speednow*dt.

image

When I am able to ride again, I want to do more empirical testing.

To that end, I'm wondering if it would be easy and clean for you to add an undocumented option to use trainer speed instead of virtual speed so that I can continue to compare the two? If not, then no problem - I will find other ways to verify the physics...

Cheers!

GCuser99 commented 2 years ago

@dvmarinoff, at the risk of you tiring of me... I realized yesterday that the virtualSpeed calculation, which is based on the physics.stackexchange link above, and appears to be same as GC's version of simulated speed, uses an approximation that actually has an exact closed-form solution. As such, I refactored virtualSpeed into virtualSpeedAlt! I think this is a unique approach (ie. I haven't found a similar solution elsewhere) and should be considered as long as there is not a huge computational cost increment.

In case you are interested, the approximation being made in virtualSpeed logic is to calculate Pke based on the difference between the measured power and power to sustain the previous state's speed. If the difference is non-zero then that implies an (de)acceleration. However, if there is an (de)acceleration, then all of the other power components need to be updated with the new speed. In fact, that does not happen in virtualSpeed. This results in an overestimate the magnitude of acceleration and underestimate of magnitude of speed. This could be corrected using an iterative approach but using your qubic root solver provides a closed form solution instead.

All this said, the error introduced by the virtualSpeed approximation does not look very large to me...

To demonstrate the correctness of virtualSpeedAlt, I am including the sister function virtualPowerAlt. With these two, you will be able to test that virtualPowerAlt will return the exact power input to virtualSpeedAlt, as in the following pseudo code:

powerRoundTrip=virtualPowerAlt(virtualSpeedAlt(power)) //returns input power almost exactly

Below are the new functions. I do not think you need powerToMaxSpeed any longer (but if you keep it, be sure to add drafting factor to the code).

    function virtualSpeedAlt(args = {}) {
        // This function is the "exact" version of virtualSpeed (which uses an approximation),
        // However, it yields results very close to virtualSpeed so the error in approximation is small
        // The original VirtualSpeed will over-estimate the magnitude of acceleration
        // and under-estimate the magnitude of speed.
        // This function is compatible with virtualPowerAlt...
        // Because it uses the Qubic root solver, it may (or may not!) be significantly more computationally
        // expensive than virtualSpeed. 
        // If the cost is not significant, then this should be used instead of virtualSpeed
        // To minimize error, dt should be made as small as possible

        const power          = args.power; // W
        const slope          = args.slope ?? defaults.slope; // 0.01 is 1%
        const mass           = args.mass ?? defaults.mass; // kg
        const windSpeed      = args.windSpeed ?? 0; // m/s
        const draftingFactor = args.draftingFactor ?? defaults.draftingFactor; // 0..1
        const dt             = args.dt ?? 1; // s
        const speedPrev      = args.speed ?? 0; // m/s
        let speed            = args.speed ?? 0; // m/s
        let distance         = args.distance ?? 0; // m
        let altitude         = args.altitude ?? 0; // m

        const cosBeta = CosBeta(slope);
        const sinBeta = SinBeta(slope, cosBeta);

        const c1bl = c1bearingLoss;
        const c2bl = c2bearingLoss;
        // const crv  = use.dynamicCrr  ? 0.1    : 0;

        // Here we expand Pke=.5*m*(vf^2-vi^2)/dt from Martin et al to get zero'th and 2'nd order ke terms
        const c0ke      = -.5 * (mass + wheelInertia) * Math.pow(speedPrev, 2) / dt;
        const c2ke      = .5 * (mass + wheelInertia) / dt;  

        const c1grav    = g * mass * sinBeta;
        const c1roll    = g * mass * crr * cosBeta;
        const c1air     = 0.5 * (CdA + spokeDrag) * rho * (Math.pow(windSpeed, 2)) * draftingFactor; // vWind ^ 2
        const c2air     = (CdA + spokeDrag) * rho * windSpeed * draftingFactor;
        const c2dynroll = crv * cosBeta;
        const c3air     = 0.5 * (CdA + spokeDrag) * rho * draftingFactor;

        const c0 = -power * (1 - drivetrainLoss) + c0ke;
        const c1 = c1grav + c1roll + c1air + c1bl;
        const c2 = c2air + c2bl + c2dynroll + c2ke;
        const c3 = c3air;

        const roots = Qubic(c3, c2, c1, c0); // roots = Qubic(c2 / c3, c1 / c3, c0 / c3);

        speed = roots[0]; // original index is 1, maybe it's 0
        let thisSpeed;

        for(var root of roots) {
            thisSpeed = root;
            if(speed > 0) {
                if((thisSpeed > 0) && (thisSpeed < speed)) {
                    speed = thisSpeed;
                }
            } else {
                if(thisSpeed > speed) {
                    speed = thisSpeed;
                }
            }
        }

        if(speed < 0 || isNaN(speed)) speed = 0;
        const acceleration = (speed - speedPrev) / dt;
        const dx = speed * dt;
        const da = dx * sinBeta; // Math.sin(Math.atan(slope));
        distance += dx;
        altitude += da;

        if(altitude < 0) altitude = 0;

        // const acceleration = powerKE / (mass + wheelInertia);
        // speed        = speed + acceleration * dt;
        // if(speed < 0) speed = 0;

        return { acceleration, speed, distance, altitude };
    }

    function virtualPowerAlt(args = {}) {
        // This function calculates power compatible with virtualSpeedAlt.
        // It calculates the kinetic energy power component based on Martin et al 
        // using the parameterization Pke=.5*m*(vf^2-vi^2)/dt.
        // Note that user must pass speedPrev and dt.
        // To minimize error, dt should be made as small as possible.

        const slope          = args.slope ?? defaults.slope; // 0.01 is 1%
        const mass           = args.mass ?? defaults.mass; // kg
        const windSpeed      = args.windSpeed ?? 0; // m/s
        const draftingFactor = args.draftingFactor ?? defaults.draftingFactor; // 0..1
        const dt             = args.dt ?? 1; // s
        const speedPrev      = args.speedPrev ?? 0; // m/s
        const speed          = args.speed ?? 0; // m/s

        const cosBeta = CosBeta(slope);
        const sinBeta = SinBeta(slope, cosBeta);

        const gravitationalResistance = g * mass * sinBeta;
        const rollingResistance       = g * mass * cosBeta * crr + speed * crv * cosBeta;
        const windResistance          = 0.5 * rho * (CdA + spokeDrag) * Math.pow((speed + windSpeed), 2) * draftingFactor;
        const bearingLossResistance   = use.bearingLoss ? (0.091  + 0.0087 * speed) : 0;
    const keResistance            = .5 * (mass + wheelInertia) * (Math.pow(speed, 2) - Math.pow(speedPrev, 2)) / (speed * dt);
        const totalResistance =
              gravitationalResistance +
              rollingResistance +
              windResistance +
              bearingLossResistance +
              keResistance;
        let power = totalResistance * speed/(1 - drivetrainLoss);

        return power;
    }

Below is my testing code from Observable. Sorry this has evolved so slowly... physics.txt

dvmarinoff commented 2 years ago

@GCuser99 Thanks for the help? I think the model is going very well, and I hope you can excuse my below average understanding of physics.

I am working on the suggested improvements and just pushed an update to the distance branch it's not very stable right now, but has the following:

I am going to do some tests with the new function tomorrow, but meanwhile a was thinking of asking what default values do you think would be suitable for the simulation parameters?

I am using the default values that the ANT+ FE-C spec suggests, but they seem to be a bit conservative. They are: dragCoefficient: 1, frontalArea: 0.4, resulting in CdA: 0.4, crr: 0.004, rho: 1.275. Zwift and RGT seem to be using way more 'faster' ones.

Cheers!

GCuser99 commented 2 years ago

You are making awesome progress - thanks!

I'm curious - how do you know Zwift and RGT are using faster parameters?

In the past, I did some work on reverse-engineering Zwift's model using the ZwiftInsider data. They ran a bot through the game at different power, body weights, body heights, bikes, and wheels. I downloaded the actual Strava streams for those tests and did many model fits on the time series data. I found pretty convincing evidence that at that time (2016-2017) they were using either Crr=.004 without dynamic rolling resistance, or Crr=.0033 with dynamic rolling resistance (the later fit slightly better but not enough to be highly confident that the Zwift model includes dynamic rolling resistance). So I think Crr=.004 without dynamic rolling resistance is reasonable based on Zwift dynamics. As I have mentioned before, I have yet to find a paper reference for dynamic rolling resistance and thus would be hesitant to build that in.

Here is what I extracted for CdARho from that study for a body height=183 cm:

<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:x="urn:schemas-microsoft-com:office:excel" xmlns="http://www.w3.org/TR/REC-html40">

Body Weight (kg) | CdARho -- | -- 45.0 | 0.298 50.0 | 0.312 62.5 | 0.345 75.0 | 0.376 83.3 | 0.393 100.0 | 0.430 125.0 | 0.515

So yes, based on those numbers a CdARho = .51 (your current default) would be on the higher side relative to what I extracted from Zwift.

The only other anecdotal evidence that I can share... I did some CdA tests on my mountain bike by repeatedly riding up and down a paved hill at different speeds without using brakes and then solving for Crr and CdA. My measured CdA was about .4 and rho=1.0 (I live at high altitude). I weigh about 72 kg and am 180 cm tall. Since it was on a mountain bike, I was not in a very aerodynamic position - just resting hands on to the handlebars in a normal mountain bike posture. At sea-level, that would correspond to a CdARho=.51. But it was a mountain bike, not a road bike, so maybe a CdARho=.51 is on the higher side for simulating a road bike?

Drivetrain loss should be zero, in my opinion.

Cheers!

GCuser99 commented 2 years ago

Have you seen this?

http://www.trainingandracingwithapowermeter.com/2011/04/estimation-of-cda-from-anthropometric.html

Dr. Andrew Coggan mentions that Cd=.7 is a good number given the imprecision of study data.

So using round numbers so not to pretend there is a lot of precision:

CdA=.7*.4= .28
CdARho (sea-level) = .28*1.275= .36

That number would put you close to the middle of my Zwift table above...

mjunkmyjunk commented 2 years ago

@GCuser99 what is the % contribution of KE towards the equation? From the Pic, it seems like it's really v small % (2.3w).

Screenshot 2022-03-03 at 9 52 53 PM
GCuser99 commented 2 years ago

Yes @mjunkmyjunk , in that particular example from Martin et al, it is small because acceleration is small. If acceleration is zero (steady state), then powerKE=0. If magnitude of acceleration is large, then magnitude of powerKE is large, and can be negative in case where acceleration is negative (such as coasting).

Modelling the effects of inertia is very important for a realistic simulation. For example, without accounting for powerKE, if one stops pedaling (power=0), then speed=0 (a potential whiplash situation! 🙂).

Anyway, I don't know of any compelling reason to ignore it, especially if we have the maths to include it.

dvmarinoff commented 2 years ago

I’ve been working on a graphing solution for distance intervals and I realised that there are different conceptual ways in which this feature can be implemented.

I see it as 3 main options:

Having all 3 would be hard, given a choice which one would you pick?

NOTE: Option 1 is about 90% done, Option 3 is like 30% done, but you can't have both.

That's a preview of the requested K2C Cycle Tour, executed as one large set of distance intervals after a small ERG warmup interval: Screenshot 2022-03-07 at 21 52 01

GCuser99 commented 2 years ago

Wow - very cool! In your K2C Cycle Tour example graphic above...

I assume the profile shown after the warmup is essentially the elevation profile versus distance on the x-axis?

After the ERG warmup, is the rest of the workout performed in Sim (or Slope) mode, based on distance intervals which seems like it would be described by your option 2?

What would be the difference between option 1 and option 2 using K2C Cycle Tour as an example?

Thanks!

PhatWheZ commented 2 years ago

I’ve been working on a graphing solution for distance intervals and I realised that there are different conceptual ways in which this feature can be implemented.

I see it as 3 main options:

  • distance intervals inside an ERG workout,
  • Course simulation as a collection of distance intervals in a workout,
  • independent workout running (or not) while you ride a course simulation, (think Zwift + TR use case, but with a 'magical' way to switch control between the simulation and the workout)

Having all 3 would be hard, given a choice which one would you pick?

NOTE: Option 1 is about 90% done, Option 3 is like 30% done, but you can't have both.

That's a preview of the requested K2C Cycle Tour, executed as one large set of distance intervals after a small ERG warmup interval: Screenshot 2022-03-07 at 21 52 01

@dvmarinoff, that Distance model looks good! I would be happy to give it a go when its ready. I got a local 20km route that i would be keen to pop in and give a bash I do think that Option 1 would be the ideal choice from my initial request as it gives an "estimate" of the real thing. However Option 3 would be amazing IMO if one could just pull the GPX file of choice and ride that (similar to RGT magic roads)

dvmarinoff commented 2 years ago

@GCuser99 The main difference between 1 and 2 is that 1 would have both time and distance on the x-axis and altitude and power on the y-axis, which would require some sort of mapping so that both the ERG (power) and the distance (altitude) intervals are clearly visible at all possible screen sizes. With short distance intervals that's doable, but with long ones it's going to be complicated. Option 2 would separate the 2 use cases and have only distance OR time on the x-axis which would massively simplify things.

Option 3 would mean that there is always going to be some course maybe a flat default one, which the user would cover while doing an ERG workout. And the user would be able to switch the course for a different one. Here comes the quastion of what happens when the User is using the course and changes the slope with the +/-, this will affect speed so do we record the original elevation/distance or the new one.

The case of negative altitude also needs to be handled (Garmin .FIT files allow for down to -500m negative altitude). I think that negative altitude should not be allowed. Maybe the app should auto switch to 0% when the altitude is about to get negative.

On .gpx files, I don’t really like the format, it’s seriously inefficient at storing course data, and since the app is working 100% in the browser without any backend this can lead to problems resulting from high memory usage. The app already has support for .FIT file virtual activities, so I think it would be better to just extend support to .FIT courses and general .FIT files. I am currently using a custom extension of the zwift .zwo format to store the course data inside a <FreeRide> tag. As efficiency it sits in the middle, but it can store data for both ERG and Simulation.

<FreeRide Sim=“<slope> ,<duration>, …”>
<FreeRide Track=“<slope> ,<distance>, …”>

Where the Sim attribute stores a list of slope-duration pairs, and the Track attribute a list of slope-distance pairs.

This is the K2C workout/course from the screenshot:

<workout>
    <SteadyState Duration="300" Power="0.49"/>
    <FreeRide Altitude="727" Track="0,2000, 6.5,4800, -1.9,7800, 2.9,4200, -5.1,14000, 4.8,5800, -3.2,7400, 9.0,2600, 0,3400, 2.6,4800, -3.1,9800, 0.6,6600, 6.3,800, -4.2,9200, 4.6,2400, -4.8,4400"/>
</workout>
GCuser99 commented 2 years ago

Thanks for the detailed explanation @dvmarinoff.

It sounds like Option 3 is different than the others in that it is an ERG-only workout (while covering a course) whereas Option 1 allows for both ERG and sim modes in same workout, and Option 2 is distance-based sim only. If I understand that correctly, then I would vote for Options 1 or 2 so that user can choose between the uniquely different experiences of ERG and sim modes. I don't see much upside to Option 3, for all of the associated complication.

If implementing Option 1 is too complicated, then Option 2 (along with your original time-based workouts) exposes 90% of the desired functionality and satisfies the OP's request, I think.

On the negative altitude issue... If user does not specify the "Altitude" start value (or even if they do), then couldn't you calculate in the background a start value that insures positive altitude as in following pseudo code:

minAltitudeAllowed = 0
altitude = startAltitude //User supplied start altitude or zero if not specified
minAltitude = altitude

for each trackPair in workout
     altitude = altitude + Math.Sin(Math.Atn(trackPair.slope/100)) * trackPair.Distance
     if altitude < minAltitude then minAltitude = altitude
next trackPair 

If minAltitude < minAltitudeAllowed then startAltitude = startAltitude + (minAltitudeAllowed - minAltitude)

It would seem that you could eventually extend option 2 to take a .FIT file course as well, memory limitations allowing. And maybe there is a way to simplify the info from .FIT file in a pre-calculation phase using some sort of simplification algorithm (such as Douglas-Peucker) that finds a reduced trackPair representation satisfying a memory limitation objective?

I cannot wait to try this - sadly it will be weeks before I can get back on my trainer... 😒

GCuser99 commented 2 years ago

On the polyline simplification algorithm - I realized that I had already written code in VBA for simplifying (reducing) X-Y polygons - so I quickly ran a Strava distance-elevation profile through same algorithm (Douglas-Peucker) using an epsilon=1 meter. The data reduction using that epsilon was 27-to-1. Epsilon could be set to achieve more or less compression. Here is a graph of the original and simplified:

image

Here are some links if/when you are interested:

https://towardsdatascience.com/simplify-polylines-with-the-douglas-peucker-algorithm-ac8ed487a4a1 https://martinfleischmann.net/line-simplification-algorithms/ http://www.diva-portal.org/smash/get/diva2:444686/FULLTEXT01.pdf https://github.com/mourner/simplify-js http://mourner.github.io/simplify-js/ (nice demo of above!)

Cheers!

mjunkmyjunk commented 2 years ago

Thanks for the link, this makes more sense reading this after reading the dougles-peucker algo

https://www.gpsvisualizer.com/tutorials/track_filters.html

I think in terms of altitude data, it's gonna be a hit and miss. Many factors at play there. eg: https://ridewithgps.com/help/grade-and-elevation

Had used this https://www.potter.ca/Biking/smoother-beta/index.html prev for "magic road"

https://www.youtube.com/watch?v=u_SaHZJ7Bns

GCuser99 commented 2 years ago

Thanks for the links @mjunkmyjunk. Just to be clear - the intent of my post above was to offer ideas on data reduction in case memory is a problem. I would not use a simplification algorithm like Douglas-Peucker for the purpose of "conditioning" a noisy course track in prep for turning it into a workout. I agree that there are generally some problems to solve in terms of elevation accuracy (of which I do have some ideas on how to address!) but that is a different conversation....

dvmarinoff commented 2 years ago

GPX

My issue issue with .gpx is that the format is a text file that uses incredibly verbose xml syntax, which results in a very very long xml text file even for a simple course. .gpx is measured in MB and depending on the course profile, features, length, sampling rate, etc. it can reach 10s of MB, .FIT is binary and you can have a ton of data and still be in the kB range.

Elevation

I was wondering more about what to do with the free form slope mode now that the app records elevation and has speed from simulation parameters. Should the app enforce that the User first needs to climb before they can descent, which kinda defeats the purpose of the free form mode, or just limit altitude to never be negative, which may result in the User going at descent speeds while recording a flat course.

Options

I think I am going to go with option 2 for now, it can later be evolved into 1 or 3. Will wait for more opinions to see which way to go.

Thanks for the links they will come in handy, I am currently exploring WebGPU and trying to learn Blender for a future side project 😃

dvmarinoff commented 2 years ago

An update:

Now developer version 0.1.31 allows for .fit ride or run activities to be used as sources for course simulation. You can also import .fit course exports made with Garmin Connect or ridewithgps. Those are upload with the usual button at the bottom of the workouts tab.

The quality of the altitude data remains a big issue for this kind of functionality so the files are processed by a line simplification algorithm (I am using the one @GCuser99 linked here). This will result in a simplified version of the course with hopefully less frequent and erratic slope changes, but all this depends on the quality of the input file.

Here are 2 files of the same course, but the 1st one is exported from Garmin Connect and the second from ridewithgps:

Screenshot 2022-03-23 at 18 10 06

This is a local TRI course that is a 20km false flat downhill followed with 20km back up. The second course has clearly a lot of random slope oscillation that doesn't exists in reality and remains there even after processing. The first on is closer to reality, but still has an occasional moment of weirdness, where a really high grade pops up for a meter or two.

Another observation is that run activities might have much better altitude data recorded than ride activities.

Here are 2 virtual courses from RGT and one real world run activity:

Screenshot 2022-03-23 at 19 22 19

The flat Borrengo is translating well, the Hilly Iron Horse is mostly good too, but the 4 hills in the first 3rd of the course are loosing the down grades to a gentle false flat. The run activity though looks really well.

There are some things that don't work yet like:

So if you try out the feature know that page reload will result in a loss of data and position. Meanwhile I'll keep working on it. Thanks for the help, everyone!

GCuser99 commented 2 years ago

@dvmarinoff I'm back to being able to perform some short tests. Looks like you have made much progress!

On the setup page of Flux-Dev I see the new SIMULATION toggle at the bottom:

If I toggle SIMULATION to "Power", will that use trainer-calculated speed? And for "Speed" will it use the Flux calculated speed?

How does the SIMULATION toggle interact with the selectable Speed buttons under CONTROLLABLE and SPEED AND CADENCE? Should I be selecting which of those to use as well?

Does the SIMULATION toggle work in both Erg and Slope modes?

Thanks in advance - looking forward to doing a lot of tests and trying to slowly get back into shape!

dvmarinoff commented 2 years ago

Hey, @GCuser99, welcome back, nice to read you are getting better now!

Meanwhile I merged the work at the distance branch on the development branch.

The Simulation toggle will switch the source for speed between speed reported from a connected device on one hand and speed calculated from power on the other. Setting it to value speed will read and record speed from a chosen connected device, setting it to power will ignore the connected device speed and calculate speed from the reported power of a connected device.

It works independent of the modes, so it applies to all of them.

Here is an attempt at a diagram:


flowchart LR
    A[Device] --> B[Speed] & C[Power]
    B --> D([db.speed]) ----> H{{toggle}} --> J[UI / Recording]
    C --> E([db.power]) --> F[VirtualStateModel] --> G([db.virtualSpeed]) --> H
GCuser99 commented 2 years ago

Ok got it - nice graph! If I want to do a simple test where I ride in distance slope mode at 3% grade for the first 1000m, and 0% grade for the next 1000m, which of these syntaxes (if any) will get me there?

<workout>
    <FreeRide Track="3,1000, 0,1000"/>
</workout>
<workout>
    <FreeRide Track="3,1, 0,1"/>
</workout>
<workout>
    <FreeRide Sim="3,1000, 0,1000"/>
</workout>

I rode with the first of those above and did not get the expected change of slopes after 1000m - it just stayed at 3% grade the whole way. I will play around some more and then provide the usual browser log. Thx!

dvmarinoff commented 2 years ago

I decided to remove (at least temporary) the .zwo course syntax. The time based one with the Sim attribute still works, but the distance based Track is off. Distance + ERG intervals resulted in too much complexity.

The app supports now .fit file uploads for course simulation. The only way I can think of, right now, for scripting a course is to produce a .fit file. My fit-file project supports uploading a .fit file as .json and exporting it to binary .fit. So if you can edit the jsFIT format, which is just a fit file as JSON and import it to the site you will be able to download back a binary .fit which is importable into flux as a course simulation file.

This archive contains a minimal .fit file with 1000m at 3% and another 1000m at 0%, you can import this into flux-devel: test-course-files.zip

The .json file you can edit as a text file by adding more record data messages with altitude and distance values. To format the values you may use this formulas:

altitude = (altitude * 5) + (500 * 5), where altitude is in meters
distance = (distance * 100), where distance is in meters

and this is a data record message:

  {
    "type": "data",
    "message": "record",
    "local_number": 3,
    "fields": {
      "timestamp": 1017595716,
      "distance": 200000,
      "altitude": 2650
    }
  }

distance is 2000m and altitude is 30m.

GCuser99 commented 2 years ago

@dvmarinoff, I ran the course fit file that you prepared for me through Flux twice - thanks for the help. The simulations went very well and graphics/display are excellent - this is an awesome app!

Looking at the results, I noticed that the output .fit files included distance and altitude, but not speed. Is that intentional? I also ran some Erg tests (virtual speed calculated from power) and see that those results did include speed in the .fit files...

Will provide more feedback later... Again, great work!

dvmarinoff commented 2 years ago

Yes that's totally a bug. Thanks for catching it!

I've managed to reproduce the speed being recorded as 0 when the simulation toggle is set to speed. Will take a closer look later today.

GCuser99 commented 2 years ago

Ok. I ran the course in both Sim modes (speed and power) and did not get a speed channel in the .fit output file for either of those modes.

Also, for the Erg runs I did, an output speed channel was produced for Sim=power, but not for Sim=speed.

In all of the runs above the speed displayed in the Flux app just fine. Thx!

dvmarinoff commented 2 years ago

Pushed some quick fixes, recording should work correct now. There is still one issue with the toggle, it always displays in the view its' default value power on page reload despite the actual value in db. I am working to fix that.

BTW in my manual testing I am using a simulation script that broadcasts mock data as if a trainer is connected. It broadcasts the current power target as power or a default value + speed, heart rate and cadence. I keep the script in the repo, but at a place where it cannot be loaded and the code that initializes it is also commented.

If you are comfortable with running the app locally on your machine I can give you instructions on how to switch it on for your local version. You will need node.js and npm installed at a minimum and optionally git.