# Calculating Cd and Crr coefficient experimentally

1. Feb 25, 2014

### jumpjack

On this site I found how to calculate Cd and Crr coefficient from experimantal data:
http://www.instructables.com/id/Measure-the-drag-coefficient-of-your-car/

I am trying to implement it in a web page and in a smartphone application, so I must "replace" Excel Solver by an algorithm, but I'm very in trouble.
I figured out how to calculate coefficients of fitting polynomial, so I now have the equation which approximates the speed curve, in form of y = ax^2+bx+c

But now how do I correlate this equation to real data to obtain Cd and Crr?

The excel sheet relies on this formula:
v1 = v0 - a*t
where:
a = F/m
and
F = air drag + roll friction
air drag = 0.5 * r *A * Cd * v^2
roll friction = m * g * Crr

so I think it is:

v1 = v0 - 0.5 * r * A * Cd * v0^2 / m - g * Crr

if this is correct, I can "zip" to:

v1 = v0 + K v0^2 + H

and

K v0^2 + H = v1-v0

My idea is to use the next sample to build an equation couple which I can solve for K and H, and get Cd and Crr from them... but it looks like it does not work, what am I doing wrong?

K v0^2 + H = v1-v0
K v1^2 + H = v2-v1

2. Feb 26, 2014

### voko

The equation of motion is $\dot v = a (1 + \alpha^2 v^2)$, where the dot means differentiation w.r.t. time, and $a$ and $\alpha$ are obtained from the constants of the problem. This equation can be integrated: $$\int\limits_{v_i}^{v_f} {dv \over 1 + \alpha^2 v^2} = aT = \frac 1 \alpha \int\limits_{\alpha v_i}^{\alpha v_f} {d u \over 1 + u^2} = \frac 1 \alpha (\arctan \alpha v_f - \arctan \alpha v_i)$$ where the i and f subscripts mean "initial" and "final", respectively, and $T$ is the time between them. We can further rewrite that as $$a \alpha T - \arctan \alpha v_f + \arctan \alpha v_i = 0$$ Now assume you have made $N$ measurements of $v_{in}, v_{fn}, T_n$, where $n$ goes from 1 to $N$, inclusive. For all those measurements you should have, in theory, the above equation exactly right. But it won't be exactly right, because your measurements are not infinitely precise, and the equation of motion is not exactly exactly right to begin with (it is a model). So instead of the zero on the right hand side you will get some error: $$a \alpha T_n - \arctan \alpha v_{nf} + \arctan \alpha v_{ni} = \epsilon_n$$ Now you can square all those errors, sum them, and find $a$ and $\alpha$ that minimise the sum: $$\sum\limits_{n = 1}^N \left[ a \alpha T_n - \arctan \alpha v_{nf} + \arctan \alpha v_{ni}\right]^2 \to \min$$ Once they are found, you will obtain your coefficients from them.

3. Feb 26, 2014

### jumpjack

I found this "suggestion" in several sources, but no examples at all, so I don't understand how to proceed!
Is there a "minimisation algorithm" somewhere? Or should I use "try&fail" method?

4. Feb 26, 2014

### voko

5. Feb 26, 2014

### Staff: Mentor

There are many minimization algorithms with different strengths and weaknesses. It's a standard topic in textbooks on "numerical methods" or "numerical analysis".

6. Feb 26, 2014

### jumpjack

7. Feb 26, 2014

### voko

I do not know whether there is such a thing. Anyway, for the mobile application there are C/C#/Java libraries, so you should be all set. For a web page, you can run the computation server side and use whatever library that works for you.

8. Feb 26, 2014

### jumpjack

I found this version for VBA:

Code (Text):
Attribute VB_Name = "Lebenber_Marquart"
Option Explicit
Option Base 1

Sub LMNoLinearFit(x() As Double, y() As Double, c() As Double, ch2 As Double, Niter As Long, Optional Itmax As Long = 10000)
' Rutina para calcular los valores de ajuste por minimos cuadrados a un modelo Fun(x)
' mediante el algoritmo de Lebenberg-Marquart
' Los valores de arranque son c
' ver. 14.04.2006 by Luis Isaac Ramos Garcia and Foxes Team
' Bibliografia
' Numerical Recipies in Fortran77; W.H. Press, et al.; Cambridge U. Press
'
Dim i As Integer, k As Integer, kmax As Integer
Dim Mcoef() As Double, det As Double, tpc() As Double
Dim Diffc() As Double       ' Vector de diferecias entre antiguos c y nuevos
Dim delta As Double         ' Valor delta de algoritmo Levenberg-Marquardt
Dim tdelta As Double        ' relative increment
Dim mu As Double            ' reduction/amplification quadratic error
Dim ch20 As Double          ' quadratic error previous step
Dim delta0 As Double        ' increment previous step
Dim tol As Double           ' machine tolerance
Dim fact As Double          ' step factor

ReDim Diffc(UBound(c) - LBound(c) + 1), tpc(LBound(c) To UBound(c))
delta = 0.001  '0.1%
k = 0
Niter = 0
fact = 10
tol = 2 * 10 ^ -16
kmax = Itmax / 4
ch20 = chi2(x, y, c)

Do While k < kmax And Niter < Itmax
Call CalculateMcoef(x, y, c, Mcoef, delta)
Call SolveLS(Mcoef, Diffc, det)
For i = LBound(c) To UBound(c)
tpc(i) = c(i) + Diffc(i - LBound(c) + 1)
Next
ch2 = chi2(x, y, tpc)

If ch2 > tol Then mu = ch2 / ch20

tdelta = (ch2 - ch20)
If ch2 > tol Then tdelta = tdelta / ch2

'DEBUG <<<<<<<<<<<<<<
'        Debug.Print Niter, delta, fact
'DEBUG <<<<<<<<<<<<<<<<<<

If mu > 10 Then
delta = delta * fact
Else
If Abs(tdelta) < 0.001 Then k = k + 1
delta = delta / fact
For i = LBound(c) To UBound(c)
c(i) = tpc(i)   'update new point
Next i
ch20 = ch2
End If

'check delta oscillation
If Niter Mod 4 = 0 Then
If delta0 = delta Then
fact = 2  'relaxed step
Else
fact = 10 'fast step
End If
delta0 = delta
End If

Niter = Niter + 1
Loop

ch2 = chi2(x, y, c)  'best approximation

End Sub

Private Function chi2(x() As Double, y() As Double, c() As Double) As Double
Dim i As Integer, f() As Double
Funct f, c, x       'take all function values
chi2 = 0
For i = LBound(x) To UBound(x)
chi2 = chi2 + (y(i) - f(i)) ^ 2
Next
End Function

Sub CalculateMcoef(x() As Double, y() As Double, c() As Double, Mcoef() As Double, delta As Double)
' un modelo Fun(x) mediante el algoritmo de Lebenberg-Marquart
' Dfun es la rutina que calcula el vector gradiente
Dim i As Integer, l As Integer, k As Integer
Dim dv() As Double                             ' Vector gradiente de la funcion de ajuste
Dim temp As Double
Dim f() As Double
ReDim Mcoef(UBound(c) - LBound(c) + 1, UBound(c) - LBound(c) + 2)
Funct f, c, x       'take all function values
DFunct dv, c, x     'take all derivatives values
For i = 1 To UBound(x) - LBound(x) + 1
For k = 1 To UBound(c) - LBound(c) + 1
Mcoef(k, UBound(c) - LBound(c) + 2) = (y(i) - f(i)) * dv(i, k) + Mcoef(k, UBound(c) - LBound(c) + 2)
For l = k To UBound(c) - LBound(c) + 1
Mcoef(k, l) = dv(i, k) * dv(i, l) + Mcoef(k, l)
Next
Next
Next
For k = 2 To UBound(c) - LBound(c) + 1
For l = 1 To k - 1
Mcoef(k, l) = Mcoef(l, k)
Next
Next

For i = 1 To UBound(c) - LBound(c) + 1              'Algoritmo Marquart
Mcoef(i, i) = Mcoef(i, i) * (1 + delta)
Next
End Sub

'------------------------------------------------------------------------------------------

http://digilander.libero.it/foxes/optimiz/Optimiz1.htm

Unfortunately it's poorly commented.
As far as i can understand, it requires a function in explicit form (Y=ax^2+bx+c), which i do not have. Is this source suitable for my task? if so, how do i put y1=y0 + K y0^2 + H in this code?!?

9. Feb 26, 2014

### jumpjack

But... do I actually need to implement the algorithm??
I already found the best-fit equation for my data in form y=ax^2+bx+c , so the least square minimisation is already performed... or not?
Maybe I've just to figure out how to express (1) y1 = y0 + H y0^2 + H in (2) y=ax^2+bx+c form? Or I could get y0 and y1 from y=ax^2+bx+c equation which I already found, and obtain (1) with known values for y1 and y0. Then I could do the same for y2 and y1, and so I would have to equations and two unknown terms...

10. Feb 26, 2014

### voko

You need to make an effort and understand what I wrote in #2. The function to minimise, in terms of the page you referenced (see example 2), is $$f(x_1, x_2, c_1, c_2) = \frac 1 {c_1 c_2} \left[ \arctan c_1 x_1 - \arctan c_1 x_2 \right]$$ where $x_1 = v_i, x_2 = v_2, c_1 = \alpha, c_2 = a, y = T$.

11. Feb 26, 2014

### voko

If your objective is to produce some numbers, however bogus, maybe yes.

If your objective is to produce numbers that approximate the coefficients, then definitely not.

The problem is that when you take measurements seconds apart - and especially 10 seconds apart as in the spreadsheet - the difference scheme approximation employed there fails miserably, because your time steps (leaps!) are just too wide. You need subsecond time resolution to get anywhere near decent approximation. And you cannot get that, plain and simple.

So you have to use a more accurate approximation of motion between the rare moments when you get a velocity fix, and you need more than two measurements because mobile phone GPS/accelerometer readings are not known for exceeding accuracy.

12. Feb 26, 2014

### jumpjack

I thought this was the function to minimize:
$$\sum\limits_{n = 1}^N \left[ a \alpha T_n - \arctan \alpha v_{nf} + \arctan \alpha v_{ni}\right]^2 \to \min$$

Anyway I don't know which are the steps to minimize a function of two variables: I figured out how to determine A, B and C coefficients in the y=ax^2+bx+c equation which approximates my data, so I think I already minimized the errors... but how do I get Cd and Crr out from these results?!?

Additionally, does it exist any tool allowing symbolic calculations, so I'm sure my steps in multiplying and adding polynomials are correct?

13. Feb 26, 2014

### voko

You are correct, I should have said "the function to find the best fit for" (which results in the function to minimize from #2). The function I wrote about in #10 is formulated in the terms of http://digilander.libero.it/foxes/optimiz/Optimiz1.htm, example 2, so that you could adapt that example to your needs.

Since you have found an implementation of the Levenberg-Marquardt algorithm, you do not need to know all its details. You do need to understand what it expects from you.

You cannot. As explained in #2, the real curve is not quadratic. It is given by $$v = \frac 1 \alpha \tan (a \alpha t + \arctan \alpha v_i)$$ which looks like this: http://www.wolframalpha.com/input/?i=plot:+y+=+tan(2+-+x)+,+.45+<+x+<+2 The best you could do is fit the tangent-function curve to your quadratic curve, but that means you lose accuracy two times.

The site I just linked does. There is also commercial software doing that. Probably some freeware, too, but I cannot advise on that.

14. Feb 26, 2014

### voko

The simplest way to get your coefficients would be via two measurements of equal duration but different initial speeds. Then we have $T_1 = T_2$ and so $\arctan \alpha v_{1f} - \arctan \alpha v_{1i} = \arctan \alpha v_{2f} - \arctan \alpha v_{2i}$. Using Newton's method, $$\alpha_{k + 1} = \alpha_{k} + {\arctan \alpha v_{2f} - \arctan \alpha v_{2i} - \arctan \alpha v_{1f} + \arctan \alpha v_{1i} \over {v_{1f} \over 1 + \alpha_k^2 v_{1f}^2} - {v_{1i} \over 1 + \alpha_k^2 v_{1i}^2} - {v_{2f} \over 1 + \alpha_k^2 v_{2f}^2} + {v_{2i} \over 1 + \alpha_k^2 v_{2i}^2}}$$ where $\alpha_k$ are successive approximations to $\alpha$. Start with $\alpha_0 = 0.001$ and go till approximations become very close to one another.

The difficulty with this method is that the durations must be very accurately equal, which is probably going to be very hard to achieve in a mobile application. But this is something you could try. Repeated measurements and recalculations should show how successful your attempts are. If you see significant variation, then you will have to adopt the more complex method we discussed earlier.

15. Feb 26, 2014

### jumpjack

Before I switch to better-fit equation you provided, which will require a lot more additional study by me :-( and also additional data recording, could you please give some hints about the quadratic approximation ax^2+bx+c ? I plotted the function obtained from least squares calculation, and it fits quite good to the experimental data:

Samples:
19,44
16,76
14,31
12,08
10,35
8,92
7,55
6,18

Approximation:
19,45
16,72
14,30
12,17
10,34
8,82
7,60
6,68

Considering I already manually performed all algebric and matrix calculations which lead to A, B and C coefficients starting from experimental data, do I already have needed data/equations needed to determine Cd and Crr? Or do I need something else?

Alternatively, how could I use the formula you provided in #12 using current data, which all start from 70 kph?

16. Feb 27, 2014

### voko

You have samples of velocity at various times. Which you can interpret as a number of pairs $(v_i, v_f)$.

17. Feb 27, 2014

### jumpjack

I meant #14: you talk about 2 different start speeds so I thought you meant two series of samples starting at different speeds, but actually I now think you meant two different samples of same series.
Can you please check spelling of your formulas? I think there are some typos here and there (for example in #13; and in #10 i don't see any Y to replace)
Finally I don't understand how a and alfa relates to "my" K and H or to Cd and Crr .

18. Feb 27, 2014

### voko

If you have four different samples in a series taken at times $t_1 < t_2 < t_3 < t_4$ such that $t_2 - t_1 = t_4 - t_3$, then this is all you need. Or, if $t_3 - t_2 = t_2 - t_1$, then you need just these three samples, and the speed measured at $t_2$ serves both as $v_{1f}$ and $v_{2i}$ (and you can simplify the formula a bit in this case).

I cannot see any problem in #13. In #10, there are typos. I can no longer edit that message, so I rewrite here correctly (I hope): $$f(x_1, x_2, c_1, c_2) = \frac 1 {c_1 c_2} \left[ \arctan c_1 x_1 - \arctan c_1 x_2 \right]$$ where $x_1 = v_f, x_2 = v_i, c_1 = \alpha, c_2 = a, y = T$

There is no $y$ in the function, but for the purpose of the VBA algorithm it will be the time difference between the initial and final sample in a pair of samples.

From Newton's second law we have $m \dot v = -R - D$ where $m$ is the mass, $R$ is the magnitude of the force of resistance to rolling, and $D$ is magnitude of the force of drag. We know further that $D = d v^2$, so we have $m \dot v = -R - dv^2$ and so we have $$\dot v = -{R \over m} - {d \over m}v^2$$ So let $$a = - \frac R m$$ and $$\alpha = \sqrt {\frac d R }$$ then $$\dot v = a(1 + \alpha^2 v^2)$$ $d$ and $R$ are constructed from $c_d$ and $c_{rr}$ and other constants similarly to your K and H, but mind the signs.

19. Feb 27, 2014

### jumpjack

you'll drive me crazy with these typos... ;-)
Did you mean $$\alpha = \sqrt {\frac d m }$$ ?

20. Feb 27, 2014

### jumpjack

Wait, no there's no typo.
But you write
V'= d/m v^2 -R/m

I write:
v1-v0= K v0^2+H

How do they relate?
I don't think it's just
d/m = K
-R/m = H

even if I have constant time intervals. Or is that so?

21. Feb 27, 2014

### voko

I wrote $$\dot v = - \frac d m v^2 - \frac R m$$ Note two minus signs. You have basically $$\Delta v = K v^2 + H$$ That means that $$K = -\frac d m \Delta t$$ and $$H = - \frac R m \Delta t$$

But I have a feeling that you are just getting yourself confused where it isn't needed. Just work out $a$ and $\alpha$ in the original constants, and be done with that.

22. Feb 27, 2014

### jumpjack

Cutting, copying, pasting, adding, fixing, summarizing, shaking, adding a small icecube, and here it is what I get:

From Newton's second law we have

$m \dot v = -R - D$

where

$m$ is the mass,
$R$ is the magnitude of the force of resistance to rolling $( = mgC_r )$
$D$ is magnitude of the force of drag $( = \frac 1 2 \rho A C_d v^2 )$

We set $d = \frac 1 2 \rho A C_d$ so we have:

$D = d v^2$

so we have

$m \dot v = -R - dv^2$

and so we have

$$\dot v = -{R \over m} - {d \over m}v^2$$

So let

$$a = - \frac R m = - \frac {mgC_r}{m} = -gC_r$$

and

$$\alpha = \sqrt {\frac d R } = \sqrt {\frac {d}{mgC_r}} = \sqrt {\frac { \frac 1 2 \rho A C_d }{mgC_r}} = \sqrt {\frac { \rho A C_d }{2mgC_r}} = \sqrt {\frac { \rho A }{2mg}} \sqrt {\frac { C_d }{C_r}}$$

then

$$\dot v = a(1 + \alpha^2 v^2)$$

$$\dot v = a(1 + \alpha^2 v^2)$$
$$\frac {dv} {dt} = a(1 + \alpha^2 v^2)$$
$$\frac {dv} {1 + \alpha^2 v^2} = a * dt$$

Integrating:

$$\int\limits_{v_i}^{v_f} {dv \over 1 + \alpha^2 v^2} = a \int\limits_{v_i}^{v_f} {dt}$$
where the i and f subscripts mean "initial" and "final", respectively.

$$\int\limits_{v_i}^{v_f} {dv \over 1 + \alpha^2 v^2} = a * (t_f-t_i)$$

$$u = \alpha v \to v = \frac u \alpha \to dv = \frac {du} \alpha$$
$$v = v_i .. v_f \to u = \alpha v_i .. \alpha v_f$$

$$\int\limits_{v_i}^{v_f} {dv \over 1 + \alpha^2 v^2} = \int\limits_{\alpha v_i}^{\alpha v_f} { \frac {du} \alpha \over 1 + u^2} = \frac 1 \alpha \int\limits_{\alpha v_i}^{\alpha v_f} { du \over 1 + u^2} = \frac 1 \alpha (\arctan \alpha v_f - \arctan \alpha v_i)$$

Hence:

$$\frac 1 \alpha (\arctan \alpha v_f - \arctan \alpha v_i) = a * (t_f-t_i)$$
$$\frac 1 \alpha (\arctan \alpha v_f - \arctan \alpha v_i) = a * T$$
$$\arctan \alpha v_f - \arctan \alpha v_i = a \alpha T$$

To determine Cd and Crr, which are contained in $a$ and $\alpha$ , I have to find $a$ and $\alpha$ that minimise the sum:

$$\sum\limits_{n = 1}^N \left[ a \alpha T_n - \arctan \alpha v_{nf} + \arctan \alpha v_{ni}\right]^2$$

because this is the sum of squared errors (squared error= (experimental - model)^2 ), and I want to minimize discrepancy (=errors) between experimental data and model.

The function to minimise, in terms of example 2 in this page, which says:

is

(1) $f(x_1, x_2, c_1, c_2) = \frac 1 {c_1 c_2} \left[ \arctan c_1 x_1 - \arctan c_1 x_2 \right]$ ($c_1 x_1$ and $c_1 x_2$ or $c_1 x_1$ and $c_2 x_2$ ??? What about c3 and c4?)

where

$x_1 = v_f$
$x_2 = v_i$

(2) $c_1 = \alpha = \sqrt {\frac { \rho A }{2mg}} \sqrt {\frac { C_d }{C_r}}$

(3) $c_2 = a = -gC_r$

$y = T$

There is no $y$ in the function, but for the purpose of the VBA algorithm it will be the time difference between the initial and final sample in a pair of samples.

In red my last doubts, to which I add:
is VBA source for "Example 2. Bi-dimensional gaussian regression" suitable in this case?

why do I have to minimise $\sum\limits_{n = 1}^N \left[ a \alpha T_n - \arctan \alpha v_{nf} + \arctan \alpha v_{ni}\right]^2$ ?

As far as I can see, $a \alpha T_n - \arctan \alpha v_{nf} + \arctan \alpha v_{ni}$ is not the difference between measurement and model, it IS the model! Or not?

23. Feb 27, 2014

### voko

There are just two parameters in this model. I chose the second example so that you could see how to deal with two variables (x) and multiple parameters (c).

In an exact solution of the equation of motion, this should be zero (go back and compare this with the equation). So whatever non-zero value we have here is an error we want minimized.

24. Feb 27, 2014