Last time I proposed the following functional form for the aerobic component:
f = 1 / [ 1 + (t / τ)]α
For times much less than the threshold τ this saturates to 1. For times much longer, than τ it has a power-law dependence with exponent α.
Differentiating this with respect to time yields the following:
∂ f / ∂ t = -(α / τ) / [1 + (t / τ)]α + 1
In the limit of small t this has the value -α / τ, so no singularity.
This would add an additional fitting parameter to the model. A reasonable initial guess for α might be 0.5.
One issue here is that the time derivative, to first order, is is proportional to α / τ, while the value for small times depends on neither to zero order. So if the fitting algorithm is working separately on α and τ, it will have a difficult time separating the contribution of α from the contribution of τ. It would be tempted to increase or decrease both together, keeping the ratio the same. Actually this isn't based on personal insight, but rather experience: I tried fitting for this formula and this is what I observed in particularly difficult cases.
So it might be better to recast the function to solve for the ratio of these values directly, rather than solving for them separately. That leaves an additional value to solve for. I can choose one or the other, but I'll take α. So I'll call this new adjusted time constant τ' = τ / α, or solving for the old term, τ = α τ'. Then the equation becomes:
f = 1 / [ 1 + (t / ατ')]α
I differentiate this with respect to time and get:
∂ f / ∂ t = -(1 / τ') / [1 + (t / α τ')]α + 1
which no longer depends on α to first order in t.
For the Jacobian the derivative with respect to τ' is needed. That's the following:
τ' ∂ f / ∂ τ' = (t / τ') / [1 + (t / α τ)]α + 1
If I want to fit for α, I need that derivative, so for that I rewrite the formula:
f = exp( ln[ 1 + (t / α τ')] × -α )
which makes it easier to solve the derivative:
α ∂ f / ∂ α =
α [ (t / α τ') / (1 + t / α τ') - ln ( 1 + [ t / α τ' ] ) [1 + (t / α τ')]-α
This worked fairly well. I played around with the dynamic damping and even set the damping limit to sqrt(2), a fairly high value, with a transition set with η = 0.95, relatively gradual. It converged with all 29 of my datasets. With the change of variables, I changed the initial guess for τ1 from AWC / P1 to AWC / 2 P1, simply because that worked better.
Note, however, I've gone from 4 fitting parameters to 5 with the addition of α2 to what I had before: P1, τ1, P2, and τ2 (In this model, I'm using α2 to represent what I had above as α, and τ2 to represent what I had as τ').
There is the persistent issue that it almost always overpredicts power for under approximately 20 seconds. I may try to address that next.