The Vehicle Road Load Problem - An Approach by
Non-Linear Modeling
Glenn D. Thompson
United States Environmental
Protection Agency
-------
ABSTRACT
When vehicle exhaust emission tests or vehicle
fuel consumption measurements are performed on a chassis
dynamometer, the dynamometer is usually adjusted to simulate
the road experience of the vehicle. Therefore, a method
to accurately characterize the road experience of the vehicle
is a prerequisite for accurate dynamometer testing of this
nature. Commonly used methods for determining the road
characteristics of the vehicle are manifold vacuum measure-
ments, drive line torque measurements, and the coast down
technique.
In the usual coast down technique, the vehicle road load
force is calculated from a numerical derivative of the speed
versus time measurements. This study significantly improves
the coast down technique by eliminating the differentiation
of the speed-time data while maintaining the ability to treat
the forces on the vehicle as a general quadratic function of
vehicle speed. A model equation for the forces on a freely
decelerating vehicle is constructed; this equation is then
analytically transformed into an expression for the speed of
the freely decelerating vehicle as a function of time. This
intrinsically non-linear equation is then computer fitted to
the speed versus time data recorded from the vehicle.
The results of road load measurements by this method are
compared with road load measurements on the same vehicle by
the drive shaft torque method.
-------
Introduction
Vehicle exhaust emission measurments are usually
performed on a dynamometer so that the vehicle engine may
be exercised in the laboratory at diverse power outputs.
The power demand of the engine is usually determined by
first adjusting the dynamometer power absorber to simulate
the power versus speed or force versus speed characteristics
of the vehicle; and then the vehicle is operated at the
desired simulated speeds. A method to accurately char-
acterize the road experience of the vehicle is therefore
a prerequisite for accurate dynamometer testing of this
nature.
Common methods for determining the road characteristics
of a vehicle are manifold pressure measurements, drive line
torque measurements and the coast down technique. Manifold
pressure measurements are the easiest to conduct, but various
engine emission control techniques, notably exhaust gas
recirculation and fuel injection have made this approach
unsatisfactory for many vehicles. Drive line torque measure-
ments are quite satisfactory, but require placing strain
gages or torque meters in the drive line of each test vehicle.
The deceleration or coast down technique is well suited for
measurements on diverse vehicles since only instrumentation
to record the vehicle speed as a function of time is required.
-------
-2-
The Deceleration Technique
The deceleration technique is perhaps the oldest
method of determining the forces acting on an automobile
as a function of the vehicle speed.*•*•* The concept
is to determine the rate of deceleration of the vehicle
and then, knowing the mass of the vehicle, the road
load force may be calculated by Newton's, second law.
Experimentally it is not practical to measure the
deceleration as a function of velocity directly, however,
it is easy to record speed as a function of time. It
is possible to differentiate the recorded velocity versus
time data to obtain acceleration versus time, and then
to use the velocity versus time data to map the accelera-
tion as a function of velocity. This is the most common
approach but it is theoretically undesirable for two
reasons. The non-analytical differentiation process is
inherently noise sensitive. This can be a problem when
attempting a least squares fit to the differentiated data,
and can lead to a falsely large linear term in the least
squares fitted curve. Also since the acceleration must
be derived from the velocity, the .initially random errors
may no longer have a normal distribution. This makes many
-------
-3-
of the usual statistical tests invalid or of question-
able validity.
Theoretically it is a better approach to assume a
model or form of the acceleration versus speed equation
and then perform analytical operations on this equation
to convert it to the form of a speed versus time equation.
This expression may then be directly fitted to the
velocity versus time data. This approach has been used
to investigate road load force equations of the usual
constant plus a velocity squared term form.[2] However,
tire datat3} indicate that terms linear in speed should be
included in the force equation and this method may be
generalized to include a linear term.
The System Energy
An expression for the system energy is constructed
to develop a generalized expression for the force acting
on the vehicle. The total energy of the decelerating
vehicle is the sum of the translational kinetic energy
and the rotational kinetic energy of any vehicle compo-
nents in rotational motion. For all mechanical components
of the wheels and drive train, the rotational velocity is
proportional to the vehicle velocity, arid the energy of
-------
-4-
the system is
E = 1/2 mv + 1/2 vli «i (1)
where:
m = the vehicle mass
v = the vehicle speed
Ii = the rotational inertia of the ith rotating
component
oti = the proportionality constant between the
rotational velocity of the i*1*1 rotating
component and the vehicle speed
Differentiating equation (1) with respect to time; and
comparing the resulting power expression with the similar
time derivative of a purely translational system; the
generalized force on the system may be expressed as:
f = (m + Ameq) v (2)
dt
The term Ameg is equal to irfl.«i2 and is the mass equiva-
lence of the inertia of the rotating components.
The Model Equation for the Vehicle Deceleration
For ideal wind free conditions the acceleration of
-------
-5-
a freely coasting vehicle is assumed to be a polynomial
function of the velocity:
dv ~
— = A(v) = -(a0 + aiv + a2v^) (3)
The ag term is usually identified as representing™]
tire losses, hpwever, both tire losses and mechanical drive
train losses probably appear in both the ag and a-^v term.
The a2y2 term of (3) is identified as the aerodynamic
term and the velocity in this term must be the vehicle
air speed. Therefore, if any ambient wind exists,
equation (3) must be written as:
A = [aQ t a^v + a2(v - W-S)2] (4)
where:
S = a unit vector in the direction of vehicle travel
w= the wind velocity vector.
This analysis only considers the effect of the component
of wind in the direction of vehicle travel. Also equation
-------
-6-
(4) is only valid when (v - w-8) £0. If the component of
wind speed in the direction of vehicle travel exceeds the
vehicle speed, then the aerodynamic drag term must be
replaced by an aerodynamic driving term by changing the
sign of a2- The effect of cross winds can be treated in
this analysis if sufficient aerodynamic informatipn is
available about the vehicle to allow the cross wind
effect to be approximated as a simple polynomial
expression of the vehicle speed.
A term must be added to equation (4) to describe
the effects of grade. The grade term is equal to g SinV
where 7 is the angle between the test surface and the
horizontal and g is the gravitational acceleration.
Since y is a very small angle the approximation of using
the grade is quite accurate. Inserting this grade term,
equation (4) becomes:
A = -[an + gs + anv + a~(v - w-S)2] (5)
where:
g = a vector in the track direction of increasing
grade with a magnitude equal to the product of
the test surface grade and the gravitational
acceleration.
-------
-7-
Expanding and regrouping, equation (5) becomes:
A = - ) [a0 + a2(ws )2 + g-s] +
[a-L - 2a2(W-S )]v + a2v2 } (6)
Integration and Inversion
Using A = dv/dt arid separating variables of
equation (6):
dv
_dt =
[a0 + a2(W-S)2 + g-S] + [ai - 2a2(w-S)]v + a2v2 (7)
Equation (7) can be integrated in closed form.
V
The integrals of equation (2) will depend on the
relative magnitude and sign of ag, a^, a2, w-S and g-s.
The appropriate choice for the form of the integrals is
more easily seen if the variable transformation:
2u - [a-, - 2a0( w-s )]
v= — ± _ - (8)
2a2
and
4a2[ a0 + a2(ws)2 + (g-s)] - [a-^ - 2a2(W-S)]2 (9)
B2 =
-------
-8-
are applied to equation (7). This yields:
dt = du/(B2 + u2) , (10)
Because of the constraint [v - (W-S) ] > 0, u will always
be real and non-negative. The B is not so strongly
restrained, B2 will always be real but it may be either
negative, positive or zero.
Case A; B2< 0
If B2 is negative the transformation:
B2 = -D2 (11)
yields:
dt = -du/(u2 - D2) , (12)
Now both u2 and D2 are real and positive, hence both u
and D are real.
The integral of equation D-12 can take either of
two forms:
-------
-9-
1 ,
t = -[tanh -"-(u/D)]* C, (13)
D -1
or
_i
t = -tcoth •L(u/D)]+ C2 (14)
D
The terms GI and C2 are constants of integration, deter-
mined by the initial conditions.
The choice of equation D-13 or D-14 to represent the
coast down data must depend on the physical possibilities
of the motion of the system. If equation D-13 is used,
the magnitude of u cannot exceed the magnitude of D since
the argument of the tanh "•*- function is bounded between
+1 and -1. But, if u is less than D, then equation D-12
shows dt/du is positive. Hence, considering the transfor-
mations 8 and 9, dv/dt is positive and the vehicle must be
accelerating. Consequently equation D-14 describes a
"coast-up" which can only occur if there is a grade or
wind driving term of sufficient magnitude that the vehicle
freely accelerates from its initial speed to some higher
steady state speed.
If equation 14 is to be used to describe the coast
down data,u must be greater than D since the magnitude
-------
-10-
of the -argument of the coth "•*• function must exceed 1.
The physical interpretation of this can be seen by re-
writing equation D-12 as an expression for du/dt.
du/dt = -(u2 - D2) (15)
If the vehrcle is initially at some velocity:where u > D,
the rate of change of u is negative hence the vehicle
will decelerate until u = 0 or u = D. If the u = D
condition "occurs, du/dt vanishes and the vehicle will
remain at the velocity where u = D. Physically this means
there must be a grade or wind .effect acting as a driving
term. The vehicle will decelerate until the drag forces
are balanced by this driving term,, and the vehicle will
remain in motion at a constant speed.
Case B; B = 0
The case B = 0 is a special case that can theore-
tically occur, however in practice it will probably .never
be observed because of .experimental error. It is pr»e- '
sented since it is the transition between two classes of
motion, and numerical analysis difficulties may Occur in
this area.
If: B = 0
-------
-11-
dt = -du/u2 (16)
and:
t = 1/u + C3 (17)
where C-j is the constant of integration to be determined
by the initial conditions of the system.
Case C; B2 > 0
In the case B2 > 0 equation D-9 may be integrated
in the form:
-t = (1/B)[tan -1 (u/B)] + C4 (18)
or:
t = (1/B)[cot ~1 (u/B)] + C5 (19)
Again C4 and C$ are constants of integration.
Equations 18 and 19 are equivalent. This may be shown
by inverting equation 18 to yield:
-------
-12-
u = B tan t-B(C4 + t) ] (20)
Using the trignometric identity relating the tangent
and cotangent functions, equation 20 may be written as:
u = B cot (Bt + V2 + BC4) (21)
Inverting equation 19 yields:
u = B cot [B (t-C5)] (22)
Comparing equations 21 and 22 these equations differ only
by the sign of the constants of integration, and by the
phase angle of »/2. Either equation may be used to
describe the vehicle coast down. However, since the B2< 0
case requires the coth -1 solution, equation 22 will be
used for the case B^ > 0 because of the convenient parallel
nature of the cot and coth functions.
»
Applying the original transformations 8 to 22 and
solving for v yields:
1 ' '•'•'-
v = - [B cot (Bt - C5) - ax/2 + a2(w-5)] (23)
3. o • . *•-."•
-------
-13-
where:
C5 = - cot -1 [ (2a2v0 +
The velocity VQ is the initial vehicle speed at t = 0.
Applying the same operations to equation D-13 yields:
1
v = — [D coth (Dt - C2) - a,/2 + a2(ws)] (24)
•a *• J-
a2
and
C2 = - coth 1 [(2a2vQ + a.|J
Either equation 23 or 24 must be fitted to the speed
time data by the method of least squares. The values of
the coefficients a^, a-., and a~ will determine which
expression, 23 or 24 is the appropriate form to describe
the data.
*
The^fcitting Technique
The least squares method requires the sum of squares:
S = £[V;L ~ V(ti)] (25)
i
-------
-14-
be minimized with respect to the fitted parameters. To
simplify the fitting process, the change of variables:
b0 = a0 + a2(w-S )2 + g
b2 = a2 (26)
b3 = the constant of integration
is made in equations 23 and 24. The equation for the
vehicle velocity may now be written as
v = v(t,b) (27)
where b designates the four parameters of the fitting
process, bQ, b, , b2, and b.,.
Equations 23 and 24 are transcendental and intrin-
sically non-linear in b; hence, an iterative approxima-
tion method must be used somewhere in the fitting process,
The equations expressing the necessary conditions for a
minimum of equation 25 may be constructed analytically,
but then this set of equations would have to be solved
numerically by some iterative approximation method. An
alternate approach, which will be used, is to approximate
-------
-15-
v(ti,b) by a Taylor series in b, about some trial point,
and then iterate the entire fitting process until the
system converges on a value for b.
Expanding v about the trial point b^; and retaining
only the terms through the first derivative:
v(tfb) = v(t,bu) +
k=0
avitd
(28)
Substituting equation 28 into equation 25
S(b) =
v(ti|b°) -E
1 k=0
$
(29)
If S is to be a minimum, the partial derivative of S with
respect to each of the fitting parameter must vanish.
Applying this condition to equation 29;
K—U
dv(t.
= 0 (30)
b°
Equation 30 represents a system of 4 equations linear
in the four unknowns b . This system can be solved for
the b, which are then used as a new trial band the process
repeated until convergence is obtained.
-------
-16-
In the analysis of vehicle coast down data, the
non-linear cotangent and hyperbolic cotangent models
gave rise to a system which exhibited convergence
difficulties. Convergence sometimes occured slowly,
requiring many iterations, each one decreasing S(b), the
sum of squares, but the solution vector would not stabi-
lize. Sometimes the system would pscillate, each new
correction complementing the correction of the previous
step, while S(b) both increased and decreased. Finally,
divergence sometimes occurred with S(b) increasing
with each iteration.
Several techniques have been devised to circumvent
these problems. A "relaxation technique,"which amends
the correction vector by halving it whenever the
correction increases S(b) and doubling the correction
vector when it reduces S(b) is well known . Another
technique computes a weighted average of the correction
vector found from the linearization technique and the
correction vector found from the steepest descent.
Reparameterization of the model is yet another technique
which often helps, if a convenient reparameterization can
be found. Unfortunately, however, there is as yet no
general convenient way to determine a priori whether a
-------
-17-
proposed reparameterization will improve the situation.
Although any of these techniques may provide the additional
help needed to make a balky system converge, if they fail
to help; no further insight about the problem or its
possible solution is provided.
The approach used in this analysis attacks one of
the fundamental problems and provides insight into the
fitting of non-linear models. Since the theoretical con-
vergence of these systems can be proven, failure to con-
verge in practice must partially arise from a deficiency
in the tools used to compute the solution. Specifically,
'computers use only a small subset of the rational numbers,
and this must be considered when any algorithm-computer
combination is used to approximate the solution of a
difficult mathematical problem.
In virtually all cases where convergence difficulties
are encountered, the linear system of equations which must
be solved at each step in the iteration is "poorly" conditioned,
Simply stated, a poorly conditioned linear system is one
in which propagated round-off error dominates the computed
solution.
A precise definition of the condition of a system of
linear equations and an explanation of how this directly
affects the convergence of the system can be expressed
-------
-18-
in terms of the norms of the system vectors and matrices.
The norm of a vector is defined as
(31)
for p = 2 this norm is merely the length in the Euclidian
sense, when p = oo equation 31 designates the infinity
norm, which is equivalent to the magnitude of the largest
component of the vector. Analogous to Lite norm of a
vector, the norm of a matrix is define! as:
UAH = max || A x II
where (32)
llxll = 1
This norm can be viewed as the longest image of the unit
sphere under the mapping associated with A. Norms defined
in this manner" have the properties.
llcMH = | cNlMil
||M N || < IIMII HN|| (33)
Consider a system of linear equations represented
byAx=r where A is a non-singular matrix. If the components
of rare not known precisely, we would like to know the
effect of such uncertainty on the solution vector X .
Let 5 r represent the uncertainty in r. Since;
-------
-19-
(34)
Ax = r
then:
A (x + »x) = r + sr (35)
and consequently:
*x = A^r (36)
To establish a bound of the relative errors 2X/X, the norms
of ix and r are constructed from equations 34, 35, 36, and
the property ,of norms 33. After rearranging terms:
nsx n ,, , , r1 "&r"
-i. < HAlMiA || (37)
nx ii I, r n
The quantity II A II-HAH is defined as the condition
of the matrix A and is designated by COND(A).
A system of equations with a large condition number
for the coefficient matrix is said to be poorly conditioned.
Such poorly conditioned systems are difficult to solve.
The propagation of round-off error in the many multiplications
and additions involved in computing the solution vector can
lead to large inaccuracies. These problems are minimized
by using Gaussian elimination with pivoting and scaling
[ 8 ] , augmented by iterative improvement. *• *
-------
20-
These algorithms are capable of solving systems of linear
equations to machine precision accuracy. However,
equation 37 demonstrates such techniques are' not suffi-
cient to produce accurate solutions if the right hand
side vector of the system is uncertain.
Considering equation 30, the effect of the system
condition on the analysis of the coast down data may be
directly understood by defining:
r . =
o
bk " bk
a ., —
dv(ti/ b)
dv(t±,b
(38)
the system of equations depicted by 30 may then be written
in the same form as equation 34:
Ax= r <39>
where the .components of r,x and A are the r. x , .and a.,
3 K JK
respectively. As the iterative solution proceeds,
the components of r tend .to vanish since X must tend to
zero if the system is converging. However the indivi-
dual elements which comprise the sums of each component of
r remain about the same size. Thus the precision of
each component r. is given by the accuracy of the largest
term of the sum. If the calculations are performed on
-------
-21-
a t-digit machine using a base number system N, the
precision of r- may be expressed as:
~ -t
*r. = N max
-------
-22-
makes excellent initial progress toward convergence, then
suddenly slows down- and ceases to make additional progress in
reducing the sum of the squares of the errors.
Considering the right aide of equation 41; the
condition of.A is controlled by the choice of the model and
its parameterization, as is ||D||. The "noise" of the data
is represented by IIv. - v(t.,t>.) ll and can only be reduced
by improving the data collection instrumentation. The only
computational hope for improving convergence is to decrease
N by increasing the precision in computing r.
In analyzing the vehicle coast down data obtained from
a Ford station wagon, the system of equations resulted in a
+ 6
coefficient matrix with a condition number of the order 10
This poor condition caused severe convergence difficulties.
Since the solution has been programmed in FORTRAN on an
IBM 370/168 computer the precision of r could be increased
by using double precision variables in the computation.
This immediately improved the percentage of data sets which
converged, however, results were still unsatisfactory. This
computer had an additional mode of arithmetic^ called extended
precision, which was then used to compute r to a precision
of 28 hexadecimal digits. This resulted in a dramatic
improvement, and successful convergence was obtained for
almost all data sets.
-------
-23-
Data
A 1971 Ford station wagon was used as the test
vehicle. This vehicle was equipped with an incremental
digital magnetic tape recorder to record vehicle speed,
drive shaft torque and time, all at one second intervals.
The vehicle speed signal was generated by a rotary shaft
encoder in the vehicle speedometer cable. The torque
transducer was a rotary transformer torque meter installed
in the vehicle drive shaft. The time signal was generated
by a quartz crystal clock.
The skid pad of the Transportation Research Center
of Ohio, in East Liberty Ohio, USA, was used for the test
track. This facility is a multilane concrete straight
track with large turnaround loops at each end. Approx-
imately 1 kilometer of this straight track has a constant
grade of 0.5% and this section was used for all measure-
ments.
The test vehicle and tires were warmed up for 30
minutes at a steady speed of 80 km/hr. Steady state
speed and torque data were then collected for 1 kilometer
in each direction of travel on the track at approximately
9, 13, 18, 22 and 26 m/sec. Following these measurements,
speed versus time data were recorded for ten coast downs,
five in each direction of the track. The initial speed
-------
-24-
of the coast down was approximately 26 m/sec., and the
terminal speed was 9 m/sec. The vehicle was then weighed
with all instrumentation and test personnel still on board.
Wind velocity and direction were recorded manually
during the test period. A mean value for the component of
wind velocity in the direction of the test track was
calculated for the time the steady state torque data were
collected and for the time required for the coast downs.
The mean values were 0.64 m/sec. during the steady state
torque measurements and 0.85 m/sec. during the coast
downs.
Each coast down was analyzed by the non-linear curve
fitting process to yield the four b coefficients of
equations 26 . An estimate of the standard error of each
coefficient was also calculated. Equations 26 were then
used to calculate the a coefficients from the b , the mean
wind velocity and test track grade. A weighted mean
was then calculated for aQ, a, and a. using the reciprocal
of the square of the estimate of the standard error as a
weighting factor. The resulting set of acceleration coef-
ficients were then multiplied by the total vehicle mass to
yield a set of coefficients for predicting the force on
the vehicle by the equation:
f = f0 + fxv + f2v2 (45)
-------
-25-
The force coefficients calculated by the coast down method
were:
ffl = 4.4 x 102 nt
f, = -9.3 nt/(m/sec.)
f2 = 9.9 x ICf1 ,nt/(m/sec.)2
The data from each set of steady state torque measure-
ments were converted to vehicle driving force measurements
using the rear axle ratio of the vehicle and the measured
rolling radius of the tire. A regression analysis was
then used to remove any acceleration effects and to yield
an improved estimate of the vehicle driving force at the
mean speed of each data set. The estimate of the mean
force and the mean speeds were then fitted to the equation:
f = fQ + h S + fjv + f2(v - W-S) (46)
where
h = a vector in the track direction of increasing
grade with a magnitude equal to the product
• pf the test track grade, the gravitational
acceleration, and the mass of the vehicle
Equation 46 is analogous to equation 5 . The f
coefficients describe the vehicle road load on a level
track with no wind. The values of these coefficients
calculated from the drive shaft torque meter data are:
-------
-26-
fQ = 3.6 xlO2 nt
f- = -2.0 nt/ (m/sec.)
nt/ (m/sec* ) 2
The force versus speed, curves using the coefficients
calculated by each method are plotted in figure 1.
Conclusions
The drive shaft torque measurements indicated slightly
higher calculated road load forces, than did the coast down
method^although the shape of the curves are in very good agree-
ment. The difference in the predicted road load forces near the
center of the measured speed range is less than 4%. There are
several reasons why the road load calculated by the steady
state drive shaft torque method might be expected to be
slightly higher than the road load calculated by the coast
down method. For example, the energy dissipated in tires
transmitting power is greater than energy dissipation in
freely rolling tires. Because of fuel consumption, the
measured vehicle mass used in the coast down calculation was
slightly less than the true mass during each coast down. Also
the vehicle mass during the drive shaft torque measurements
would be greater than the vehicle mass during the coast
downs. However, the standard estimate of error for the
coefficients calculated by each method would indicate that the
4% difference is not statistically significant.
-------
-27-
The coast down technique can yield road load results
very similar to drive shaft torque meter measurements
when the added flexibility of a linear term is included.
The results cannot be expected to be identical since
the deceleration technique is intrinsically a transient
experiment, however good precision may be obtained with
either method. Since coast down measurements can be con-
ducted with less instrumentation of the test vehicle than
drive shaft torque meter measurements, the coast down
technique is particularly useful for measurements on a
diverse class of vehicles.
-------
VEHICLE ROAD LOAD
1971 FORD STATION WAGON
Force
(nt)
900
800
700
600
500
400
coast down method
drive shaft torque meter method
10
12
i
14
16 18 20 22
24
26
SPEED (m/sec.)
FIGURE 1
-------
Bibliography
1. Hoerner, Sighard, F. Aerodynamic Drag, The Otterbein
Press, Dayton, Ohio, 1951.
2. White, R.A. and Korst, H.H., "The Determination of
Vehicle Drag Contribution from Coast-Down Tests",
Society of Automobile Engineers 720099, New York,
N.Y., 1972.
3. Curtiss, W.W., "Low Power Loss Tire", Society of
Automobile Engineers 690108, New York, N.Y., 1969.
4. Marks, L.S. Mechanical Engineers' Handbook, McGraw-
Hill Book Company, New York, N.Y., 1941.
5. Draper, N.R. and H. Smith Applied Regression Analysis,
John Wiley and Sons, New York, N.Y., 1966.
6. Marquardt, D.W. "An Algorithm for Least Squares
Estimation of Non-Linear Parameters" Journal of the
Society for Industrial and Applied Mathematics, Vol.2,
1963, pp. 431-441.
7. Guttman, I. and Meeter, D.A. "Use of Transformations
on Parameters in Non-linear Theory, I. Transformation
to Accelerate Convergence in Non-Linear Least Squares",
Technical Report No. 37, Department of Statistics,
University of Wisconsin, Madison, Wisconsin, 1964.
8. Forsyth, G.E. and C.B. Moler Computer Solution of
Linear Algebraic Systems, Prentice Hall, Englewood Cliffs,
New Jersey, 1967.
9. Ibid.
10. Scarborough, I.E. Numerical Mathematical Analysis
The Johns Hopkins Press, Baltimore Maryland, 1966.
------- |