Automatic CdA and Crr Estimation in Aerolab

Fixes #851.

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;
+}

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);
 
 };
 

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

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();