another attempt at modeling neuromuscular power in the power-duration curve (ugh)

I decided to give one more try at modeling the maximal power curve in the neuromuscular power region, using the low-ceiling approach I described last time, sort of a post-hoc adjustment to calculated power which assures that what was originally a non-increasing function would remain non-increasing.

First I calculate an initial power as before:

P_init = P₁ (τ₁ / t) (1 - exp[-t / τ₁]) + P₂ / (1 + t / α₂τ₂)^α₂

I like this model because each parameters has an easy-to-interpret meaning. It also seems to fit fairly nicely using my iteratively weighted nonlinear least square fitting scheme, where I use weights of approximately 1 for points below the envelope curve and 10 thousand for points above the envelope curve.

On weighting I say "approximately 1" because I use 1 if the jump to surrounding points is 2.5%, but if the jump is more than 2.5% I scale the weight by ln( t_n+1/ t_n-1 ) / 2 ln 1.025. Note I don't apply this to the points falling above the envelope curve since this would unbalance the result. Note the principal region affected by this weighting scheme is the region < 40 seconds. For example, at 1 second the jump to the next time points is 100%, not 2.5%, so this point gets a weight ln 2 / ln 1.025 = 28.4 (I didn't have a preceding point for the 1 second point so I just used the jump to the following point).

Anyway, the post-hoc adjustment to this is the following, which is related to the norm function:

P = [ P_init^-k₃ + P₃^-k₃ ]^-1/k₃

where k₃ is a parameter which describes how sharply the power is truncated (1 is very gradually, 10 would be very abruptly) and P₃ is the power ceiling.

At first I thought this would wreak havoc on my derivatives but then I remembered the chainrule. If function f is a function of g and g is a function of x then the derivative of f on x can be calculated: ∂f / ∂ x = ( ∂f / ∂ g ) ( ∂ g / ∂ x ). This makes the differentiation easy.

P₃ ∂ P / ∂ P₃ = P₃^-k₃ / (P_init^-k₃ + P₃^-k₃)^{-(1/k₃ + 1)}

This expressed the derivative of P with respect to P₃ in terms of the initally defined P_init, the power without the application of this model.

Next I need to recalculate the derivative of all of the other terms. But here's where the chainrule comes in. I only need to calculate one derivative, which is:

∂ P / ∂ P_init = P_init^{-(k₃ + 1)} / (P_init^-k₃ + P₃^-k₃)^{-(1/k₃ + 1)}

Then I simply multiply the derivatives calculated for the other terms, which now apply to P_init, by this factor, yielding the derivative for P.

As always, I check the derivatives numerically to eliminate the inevitable initial mistakes, and as written above, the numerical calculation and the analytic calculation are in excellent agreement.

I could treat k₅ as yet another fitting parameter but I was afraid this would be a bit too much, so instead I chose a reasonable value and sticked with this. For example, 5 seemed to work fairly well in early tests. So I just need to worry about fitting P₃. And worry I did; note that if the power is not overestimated at short times, the optimal value of P₃ is infinite. And if I set P₃ too small, then the anaerobic power component could get squeezed out. But I figured I'd fire it up and give it a try.

And it worked.... in some cases. But in other cases it crashed and burned, not even attaining an envelope fit, but rather burying the model curve lower than the data. And in still other cases, it did absolutely nothing, because P_init was already close to the power peak.

Any model needs to not work in just theory, or in just a single case, but rather to robustly fit data of a wide quality range, so this experiment was not a success. I revert to the 5-parameter model without the power ceiling.

Out of interest, I did the following plots for this model, showing the derivative of the power on the natural log of each parameter. Note this derivative for total power is total power, which is equal to the sum of the derivatives with respect to the natural logarithms of P₁, P₂, and P₃. You can see from the plot that each power coefficient has a certain domain over which it is dominent, from P₃ up to around 10 seconds, then P₁, up to around 50 seconds, then P₂. τ₁ has a little peak of influence around 60 seconds while τ₂ and α₂ show some strength for very long times.

In contrast, here's the result for the 5-parameter model, without the attempt to model neuromuscular power. It's cleaner at short times, with less compatition for control of the power curve:

It's not obvious from these plots that the 5-parameter version should converge so much better than the 6-parameter version, but it did. In general there's an issue when multiple parameters are fighting for dominance over the same regime of the plot. Ideally each parameter has its zone of control. In this example, the P₃ term is strong for times under 10 seconds. But this is a strong function of the value of P₃. I think this is the core problem: the form of the model has P₃ lose control quickly when the value increases much higher than P₁, and similarly P₁ loses significance when P₃ becomes significantly smaller.