From b3da5c7636a45d072f4efee68de8de1e20c89945 Mon Sep 17 00:00:00 2001 From: Alejandro Martinez Date: Fri, 24 Aug 2012 12:34:30 -0300 Subject: [PATCH] Automatic CdA and Crr Estimation in Aerolab Fixes #851. --- src/Aerolab.cpp | 115 ++++++++++++++++++++++++++++++++++++++++++ src/Aerolab.h | 2 +- src/AerolabWindow.cpp | 22 +++++++- src/AerolabWindow.h | 2 +- 4 files changed, 138 insertions(+), 3 deletions(-) diff --git a/src/Aerolab.cpp b/src/Aerolab.cpp index f097f60642..e00ea2363d 100644 --- a/src/Aerolab.cpp +++ b/src/Aerolab.cpp @@ -769,3 +769,118 @@ void Aerolab::refreshIntervalMarkers() } } } + +/* + * Estimate CdA and Crr usign energy balance in segments defined by + * non-zero altitude. + * Author: Alejandro Martinez + * Date: 23-aug-2012 + */ +bool Aerolab::estimateCdACrr(RideItem *rideItem) +{ + // HARD-CODED DATA: p1->kph + const double vfactor = 3.600; + const double g = 9.80665; + RideFile *ride = rideItem->ride(); + + if(ride) { + const RideFileDataPresent *dataPresent = ride->areDataPresent(); + if(dataPresent->alt && dataPresent->alt) { + double dt = ride->recIntSecs(); + int npoints = ride->dataPoints().size(); + double X1[npoints], X2[npoints], Egain[npoints]; + int nSeg = -1; + double altInit = 0, vInit = 0; + /* For each segment, defined between points with alt != 0, + * this loop computes X1, X2 and Egain to verify: + * Aero-Loss + RR-Loss = Egain + * where + * Aero-Loss = X1[nSgeg] * CdA + * RR-Loss = X2[nSgeg] * Crr + * are the aero and rr components of the energy loss with + * X1[nSeg] = sum(0.5 * rho * headwind*headwind * distance) + * X2[nSeg] = sum(totalMass * g * distance) + * and the energy gain sums power in the segment with + * potential and kinetic variations: + * Egain = sum(eta * power * dt) + + * totalMass * (g * (altInit - alt) + + * 0.5 * (vInit*vInit - v*v)) + */ + foreach(const RideFilePoint *p1, ride->dataPoints()) { + // Unpack: + double power = max(0, p1->watts); + double v = p1->kph/vfactor; + double distance = v * dt; + double headwind = v; + if( dataPresent->headwind ) { + headwind = p1->headwind/vfactor; + } + double alt = p1->alt; + // start initial segment + if (nSeg < 0 && alt != 0) { + nSeg = 0; + X1[nSeg] = X2[nSeg] = Egain[nSeg] = 0.0; + altInit = alt; + vInit = v; + } + // accumulate segment data + if (nSeg >= 0) { + // X1[nSgeg] * CdA == Aero-Loss + X1[nSeg] += 0.5 * rho * headwind*headwind * distance; + // X2[nSgeg] * Crr == RR-Loss + X2[nSeg] += totalMass * g * distance; + // Energy supplied + Egain[nSeg] += eta * power * dt; + } + // close current segment and start a new one + if (nSeg >= 0 && alt != 0) { + // Add change in potential and kinetic energy + Egain[nSeg] += totalMass * (g * (altInit - alt) + 0.5 * (vInit*vInit - v*v)); + // Start a new segment + nSeg++; + X1[nSeg] = X2[nSeg] = Egain[nSeg] = 0.0; + altInit = alt; + vInit = v; + } + } + /* At least two segmentes needed to approximate: + * X1 * CdA + X2 * Crr = Egain + * which, in matrix form, is: + * X * [ CdA ; Crr ] = Egain + * and pre-multiplying by X transpose (X'): + * X'* X [ CdA ; Crr ] = X' * Egain + * which is a system with two equations and two unknowns + * A * [ CdA ; Crr ] = B + */ + if (nSeg >= 2) { + // compute the normal equation: A * [ CdA ; Crr ] = B + // A = X'*X + // B = X'*Egain + double A11 = 0, A12 = 0, A21 = 0, A22 = 0, B1 = 0, B2 = 0; + for (int i = 0; i < nSeg; i++) { + A11 += X1[i] * X1[i]; + A12 += X1[i] * X2[i]; + A21 += X2[i] * X1[i]; + A22 += X2[i] * X2[i]; + B1 += X1[i] * Egain[i]; + B2 += X2[i] * Egain[i]; + } + // Solve the normal equation + // A11 * CdA + A12 * Crr = B1 + // A21 * CdA + A22 * Crr = B2 + double det = A11 * A22 - A12 * A21; + if (abs(det) > 0) { + // round and update if the values are in Aerolab's range + double cda = floor(10000 * (A22 * B1 - A12 * B2) / det + 0.5) / 10000; + double crr = floor(1000000 * (A11 * B2 - A21 * B1) / det + 0.5) / 1000000; + if (cda >= 0.001 and cda <= 1.0 and crr >= 0.0001 and crr <= 0.1) { + this->cda = cda; + this->crr = crr; + return true; + } + } + } + } + } + return false; +} diff --git a/src/Aerolab.h b/src/Aerolab.h index 3aca0f9782..fa6a14e6d3 100644 --- a/src/Aerolab.h +++ b/src/Aerolab.h @@ -135,7 +135,7 @@ class Aerolab : public QwtPlot { int intRho() const { return (int)( rho * 10000); } int intEta() const { return (int)( eta * 10000); } int intEoffset() const { return (int)( eoffset * 100); } - + bool estimateCdACrr(RideItem* rideItem); }; diff --git a/src/AerolabWindow.cpp b/src/AerolabWindow.cpp index ab485d902b..9c954f7373 100644 --- a/src/AerolabWindow.cpp +++ b/src/AerolabWindow.cpp @@ -200,6 +200,9 @@ AerolabWindow::AerolabWindow(MainWindow *mainWindow) : comboDistance->setCurrentIndex(1); smoothLayout->addWidget(comboDistance); + QPushButton *btnEstCdACrr = new QPushButton(tr("&Estimate Cda and Crr"), this); + smoothLayout->addWidget(btnEstCdACrr); + // Add to leftControls: rightControls->addLayout( mLayout ); rightControls->addLayout( rhoLayout ); @@ -236,6 +239,7 @@ AerolabWindow::AerolabWindow(MainWindow *mainWindow) : connect(eoffsetLineEdit, SIGNAL(textChanged(const QString)), this, SLOT(setEoffsetFromText(const QString))); connect(eoffsetAuto, SIGNAL(stateChanged(int)), this, SLOT(setAutoEoffset(int))); connect(comboDistance, SIGNAL(currentIndexChanged(int)), this, SLOT(setByDistance(int))); + connect(btnEstCdACrr, SIGNAL(clicked()), this, SLOT(doEstCdACrr())); connect(mainWindow, SIGNAL(configChanged()), aerolab, SLOT(configChanged())); connect(mainWindow, SIGNAL(configChanged()), this, SLOT(configChanged())); connect(mainWindow, SIGNAL(intervalSelected() ), this, SLOT(intervalSelected())); @@ -466,7 +470,23 @@ AerolabWindow::setByDistance(int value) aerolab->setData(ride, false); } - +void +AerolabWindow::doEstCdACrr() +{ + RideItem *ride = mainWindow->rideItem(); + /* Estimate Crr&Cda */ + if (aerolab->estimateCdACrr(ride)) { + /* Update Crr/Cda values values in UI */ + crrLineEdit->setText(QString("%1").arg(aerolab->getCrr()) ); + crrSlider->setValue(aerolab->intCrr()); + cdaLineEdit->setText(QString("%1").arg(aerolab->getCda()) ); + cdaSlider->setValue(aerolab->intCda()); + /* Refresh */ + aerolab->setData(ride, false); + } else { + /* report error: insufficient data to estimate Cda&Crr */ + } +} void diff --git a/src/AerolabWindow.h b/src/AerolabWindow.h index dea6268b32..f8b4960219 100644 --- a/src/AerolabWindow.h +++ b/src/AerolabWindow.h @@ -59,7 +59,7 @@ class AerolabWindow : public GcWindow { void setEtaFromText(const QString text); void setEoffsetFromSlider(); void setEoffsetFromText(const QString text); - + void doEstCdACrr(); void setAutoEoffset(int value); void setByDistance(int value); void rideSelected();