Michael Puchowicz Models

.. updated to include all 3 variants of the 'Veloclinic' model
   linear, exponential and regeneration models for the second
   component
This commit is contained in:
Mark Liversedge
2014-05-03 09:54:07 +01:00
parent b035122107
commit b3f1ab1389
4 changed files with 102 additions and 25 deletions

View File

@@ -50,6 +50,9 @@
CPPlot::CPPlot(QWidget *parent, Context *context, bool rangemode) : QwtPlot(parent), parent(parent),
// model
model(0), modelVariant(0),
// state
context(context), rideCache(NULL), bestsCache(NULL), rideSeries(RideFile::watts), isFiltered(false), shadeMode(2),
shadeIntervals(true), rangemode(rangemode), showPercent(false), showHeat(false), showHeatByDate(false),
@@ -498,13 +501,7 @@ CPPlot::plotModel()
case 4:
{
cp = tau = t0 = 0;
deriveCPParameters();
// ooopsie no model for us!
if (cp == 0 && tau == 0 && t0 == 0) return;
// the Veloclinic model has the following formulation
// the Michael Puchowicz (aka @Veloclinic) model has the following formulation
//
// p(t) = pc1 + pc2
// Power at time t is the sum of;
@@ -518,37 +515,44 @@ CPPlot::plotModel()
// p1 - pmax - cp as derived from the CP2-20 model
// p2 - cp as derived from the CP2-20 model
// tau1 - W'1 / p1
// tau2 = 15,000
// w2 - A slow twitch W' derived from p2 * tau2
// alpha - 0.1
// beta
// tau2 - 15,000
// w2 - A slow twitch W' derived from p2 * tau2
// alpha- 0.1
// beta - 1.0
//
// Fast twitch component is:
// pc1(t) = W'1 / t * (1-exp(-t/tau1)) * ((1-exp(-t/10)) ^ (1/alpha))
//
// Slow twitch component has three formulations:
// a) pc2(t) = p2 * tau2 * (1-exp(-t/tau2))
// b) pc2(t) = p2 / (1 + t/tau2)
// c) pc2(t) = p2 / (1 + t/5400) ^ (1/beta)
// sprint capped linear) pc2(t) = p2 * tau2 * (1-exp(-t/tau2))
// sprint capped regeneration) pc2(t) = p2 / (1 + t/tau2)
// sprint capped exponential) pc2(t) = p2 / (1 + t/5400) ^ (1/beta)
//
// Currently deciding which of the three formulations to use
// as the base for GoldenCheetah (we have enough models already !)
//
// to keep things simple we just use formulation (a) for now.
cp = tau = t0 = 0;
deriveCPParameters();
// ooopsie no model for us!
if (cp == 0 && tau == 0 && t0 == 0) return;
double pmax = cp * (double(1.00f)+tau /(((double(1)/double(60))+t0)));
double w1 = cp*tau*60;
double p1 = pmax - cp;
double p2 = cp;
double tau1 = w1 / p1;
double tau2 = 15000;
double alpha = 0.1;
const double tau2 = 15000;
const double alpha = 0.1;
const double beta = 1.0;
//double w2 = p2 * tau2;
//qDebug()<<"model parameters: pmax="<<pmax<<"w1="<<w1<<"p1="<<p1
// <<"p2="<<p2<<"tau1="<<tau1<<"tau2="<<tau2<<"alpha="<<alpha;
// populate curve data with a CP curve of 100 points resolution
const int points = 3600 * 5; // 5 hours is enough
const int points = 3600 * 10; // 10 hours is enough
QVector<double> cp_curve_power(points);
QVector<double> cp_curve_time(points);
@@ -562,8 +566,27 @@ CPPlot::plotModel()
} else {
// two component model
double pc1 = w1 / t * (1.00f - exp(-t/tau1)) * pow(1-exp(-t/10), alpha); // ^ (1/alpha) missing XXX
double pc2 = p2 * tau2 / t * (1-exp(-t/tau2));
double pc1 = w1 / t * (1.00f - exp(-t/tau1)) * pow(1-exp(-t/10), alpha);
// which variant for pc2 ?
double pc2 = 0.0f;
switch (modelVariant) {
default:
case 0 : // exponential top and bottom
pc2 = p2 * tau2 / t * (1-exp(-t/tau2));
break;
case 1 : // linear feedback
pc2 = p2 / (1+t/tau2);
break;
case 2 : // regeneration
pc2 = pow(p2 / (1+t/5400),1/beta);
//pc2 = p2 / pow((1+t/5400),(1/beta));
break;
}
// p(t) as a sum of both component powers
cp_curve_power[i] = pc1 + pc2;
@@ -1293,7 +1316,7 @@ CPPlot::setShadeIntervals(int x)
// model parameters!
void
CPPlot::setModel(int sanI1, int sanI2, int anI1, int anI2, int aeI1, int aeI2, int laeI1, int laeI2, int model)
CPPlot::setModel(int sanI1, int sanI2, int anI1, int anI2, int aeI1, int aeI2, int laeI1, int laeI2, int model, int variant)
{
this->anI1 = double(anI1) / double(60.00f);
this->anI2 = double(anI2) / double(60.00f);
@@ -1306,6 +1329,7 @@ CPPlot::setModel(int sanI1, int sanI2, int anI1, int anI2, int aeI1, int aeI2, i
this->laeI2 = double(laeI2) / double(60.00f);
this->model = model;
this->modelVariant = variant;
clearCurves();
}