Ability to create workouts off Distance

[PhatWheZ]

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>

[PhatWheZ]

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]

@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:

• calculating speed from power
• figuring out the coexistence with ERG, and the other Slope modes
• extending the .zwo format and its parser

[PhatWheZ]

@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]

• 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]

@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]

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]

@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]

@GCuser99 Brilliant!

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

[dvmarinoff]

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

[GCuser99]

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:

• dynamic rolling resistance (http://www.kreuzotter.de/english/espeed.htm#forml)
• cadence adjustment for CdA (http://www.kreuzotter.de/english/espeed.htm#forml)
• power loss due to bearing friction
• spoke drag for air resitance
• wheel inertia for Pke
• power adjustment due to altitude hypoxia
• air density as a function of altitude

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]

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

[dvmarinoff]

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]

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]

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]

@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]

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]

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]

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]

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]

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:

• You need to pass to powerToMaxSpeed in function "test" the last state of acceleration (which is not zero)
• I completed coding and testing of the other modeling options in virtualSpeed
• For more accurate virtualSpeed and really tight convergence with powerToMaxSpeed use smaller dt

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:

• Probably don't need smallAngleApprox as an option - I put it in there so that one could see the impact of this approximation on high road slopes, but I really do not know why anyone would make this approximation with modern fast computers.
• My opinion is that dtloss should default to zero because the power is being measured at the fly wheel, not at pedal. I did check that it is working properly though...
• I assume that you are going to structure this in some way that you input to virtualSpeed all of the raw power measurements and associated dt's, and then integrate and sample to a preferred sample rate (such as 1s)? I can help do some simulations for this...

Here is the modified code: physics.txt

Thanks for letting me help - cheers!

[dvmarinoff]

@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]

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]

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

Cheers!

[GCuser99]

@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]

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]

@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.

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]

@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]

@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:

• the switch between trainer and virtual speed feature, clicking the new button in the simulation settings section labeled power/speed will switch the source for virtual state (speed, distance, altitude) between the selected power or speed source. The Speed field on the home page will now display the active one.
• persistence for source selection settings between workout sessions.
• the new virtualSpeedAlt function is added as virtualSpeedCF

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]

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:

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]

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...