Part 29: A Magical Trick

©Copyright February 2002

The Magic Formulae (*Physics of Racing* parts 21, 22, 24, 25) for grip *versus*
slip have some disadvantages for the simulation programmer. Chief among them are
the complicated mathematical structure and the large number of parameters. In
this instalment, we show a much simpler mathematical expression that mimics the
most important, overall features of the Magic Formulae with many fewer
parameters. This "Magic Trick" formula is much easier to code and
debug, requires much less computing horsepower, and differs by less than 10%
almost everywhere from Pajecka: probably sufficiently accurate for gaming
simulation. Just to tantalize you, let me offer the following two plots:

The one on the left is the full, longitudinal Magic Formula from Part 21. The
one on the right is a plot of , where *s* is longitudinal slip and *F _{z}*
is vertical load in Newtons, as before. This function is much easier to grasp
and remember than is the Pajecka formula. Also, it is plain to see that the two
plots are so similar that we should expect no gross handling differences were we
to use the right-hand in lieu of the left-hand. In the rest of this instalment,
I detail the methods for obtaining the numerical values of the three parameters.

The full Magic Formulae were developed for professional research in tyre
dynamics. That is why they have so many parameters: they must model everything
from large truck tyres to off-road tyres to racing slicks, plus every subtlety
of behaviour of every different type of tyre. For simulation and gaming,
however, the most important behaviour is simply breakaway instability: the
decrease in grip with increasing slip beyond the limits. Minor bumps and wiggles
in the curves may be critical to capture differences in brand or type or rubber
compound, but they are probably overkill even for detailed simulations like
Grand Prix Legends or NASCAR 4. Since the Magic Formulae are "the only game
in town," however, simulation programmers tend to use them whilst
purposefully wasting their subtle capabilities. In fact, my own presentations of
the formulae set most of the parameters to zero value, and mine was derived
directly from Genta's (*Motor Vehicle Dynamics*, pub. World Scientific,
1999). So the available modelling 'dynamic range,' as it were, of the Magic
Formulae, seems too large for our applications. We need a simpler formula. After
playing around for a little while, I found that the following bears an uncanny
resemblance to the Magic Formula, but with only three parameters, (*A* *B*, *P*):

This is just a generic form of the formula, in which *F*_{horizontal}
could be *F*_{lateral} or *F*_{longitudinal}, and could
be slip angle or slip ratio. In other words, the sketch above can account for
either the longitudinal (part 21) or the lateral Magic Formula (part 22). Let's
analyse the longitudinal formula in detail, using Mathematica (MMA, introduced
in *PhoRS* part 27). I make reference to the longitudinal formula from part
21 *without* copying it here, so you will need part 21 on hand to go
forward from this point. You will also need familiarity with Mathematica, and
possibly to have its manual on hand.

Encode the Pajecka *b* parameters twice as vectors, one symbolic, one as
numeric *rules* for Genta's supposed Ferrari:

`ClearAll [b] ` |

So, it's meaningful to do things like this

` b[symbolic] /. b[ferrari]`

Define some helper functions and more numeric rules (the nomenclature will be obvious if you have one eye on the equations in part 21 this is a straight transcription):

** b[n_Integer] := b[symbolic] [[n+1]]
muPeak[fz_] := b[1] fz + b[2]
bigD[fz_] := muPeak[fz] fz
bb0d[fz_] := (b[3] fz**

**
** numRulez = {MNMegaNewton, KNKiloNewton,
K1000,

Kilo1000, Mega1000000,

0.0, 1.1, -1.1,
100.100}

In our numerical exposition of the formula in part 21, we found the constant 0.0822203 popping up over and over. So, test the current forms to see if the same constant appears by substituting the Ferrari's data into them:

` `

**Simplify[bigB[aa KN] /. b[ferrari]] /. numRulez**

` 0.0822203`

Add more helpers:

` `

**bigS[sigma_, fz_] := (100 sigma) + b[9] fz + b[10]**

` bigE[fz_] := b[6] fz`

^{2} + b[7] fz + b[8]

and, finally, the whole formula:

** fx[sigma_, fz_ ] := Module[
{SB = bigS[sigma, fz] bigB[fz],
e = bigE[fz],
bigD[fz] Sin[b[0] ArcTan[SB + e (ArcTan[SB] - SB)]]]**

Make a version of the formula with the Ferrari substituted in:

** fxNum[s_, fz_] := fx[s, fz] / Newton /. b[ferrari] /.
numRulez**

Evaluate this final formula symbolically to check against part 21, then numerically, with a big plot:

` `

**fxNum[sig, zed]
(1.688 zed Sin[
1.65 ArcTan[8.22203 sig - 10. (-8.22203 sig +
ArcTan[8.22203 sig])]])**

** Plot3D[fxNum[sigma, fz Newton],
{sigma, -.20, .20}, {fz, 1, 8000}, ViewPoint{-1.3,
2.4, 1.5}] **

We've successfully reproduced the Magic Formula, and this is the plot we want to emulate. So far, there has been nothing new, just recasting of familiar territory in Mathematica, where we can manipulate the data and formulas. The new function is vastly easier to write down (that's the whole point!), and the plot looks very similar:

** fx2Num[s_, fzN_, A_, B_, P_ ] := B fzN s / (1 + Abs[A s]^ P)
Plot3D[fx2Num[sig, fzn, 9.625, 31, 2.375],
{sig, - .20, .20}, {fzn, 1, 8000}, ViewPoint {-1.3,
2.4, 1.5}] **

It remains only to show how to obtain the parameter values

(*A, B, P*) = (9.625, 31.0, 2.375).

The short answer is by the ** least squares method**. The slightly
longer answer is to subtract the two plots, sum up the squares of the
differences, and then find those values of (

First, plot the errors, that is, the difference between the two forms.

** errdiff1[sig_, fzn_, A_, B_, P_] :=
(fxNum[sig, fzn Newton] - fx2Num[sig, fzn, A, B, P])**

** Timing[Plot3D[errdiff1[sig, fzn, 9.625, 31, 2.375],
{sig, -.20, .20}, {fzn, 1, 8000}, ViewPoint{-1.3,
2.4, 1.5}]]**

`{0.812 second, - SurfaceGraphics -}`

Note that the maximum magnitude of the errors is about 800 N, or less than 10% of the maximum force. Given that the shape of the plot is very similar to that of Pajecka, hence, that the behavioural characteristics will be similar, we should feel encouraged to go on, with the reservation that the errors are large at the origin, where the force is small. In fact, the errors do reach around 32% at the origin, as shown in the next plots

** errpct[sig_, fzn_, A_, B_, P_] :=
**

` `

**TimingPlot[Plot3D[errpct[sig, fzn, 9.625, 31, 2.375],
{sig, -.20, .20}, {fzn, 1, 8000}, ViewPoint ]]**

However, this does not concern me very much, for a couple of reasons. First,
consider the fact that most authors are satisfied to use a *linear*
approximation near the origin. Witness the ever-present *cornering stiffness*
in the lateral case. Neither of our forms is particularly linear and would
probably deviate from a linear approximation at least as much as they deviate
from one another. If a linear approximation to Pajecka is good enough, why
wouldn't a handy nonlinear approximation, which happens to have much better
overall behaviour, be good enough? Secondly, in racing simulation, we are seldom
in the linear region, and we should be much more concerned with the differences
near the force maxima, around *s* = 0.15. At these values, the errors are
much smaller, well under 10%.

Turning attention back to the plot of errors, notice that some errors are
negative and some are positive. That's the reason that we square them: we want
the overall * sum* to be positive for any values of the parameters.
We don't want to add negative errors to positive errors and risk getting a low
estimate of variance. Note that we might just as well add the absolute values
rather than the squares, but squares are much easier to manage symbolically
since they are smooth (differentiable) at the origin. Also, squaring errors
amounts to an assumption of independence: it treats the errors just like
components of a vector in a large-dimensional Euclidean space. Historically, a
considerable number of tools were developed to solve such problems analytically
(symbolically), so squares it is. In fact, we

For now, in the interest of getting results quickly, we do a brute-force
search for the best values of (*A, B, P*). I like to call this
search method ** archaeological**, because it's the same one
archaeologists use to search a site for artefacts: grid off the search space,
then exhaustively examine each grid cell. Since (

First, zoom in on some two-dimensional slices of the search space and see if
we're near a minimum. Look at , the sum over from, say, -0.2 to
0.2 by 0.02 and *F _{z}* from, say, 1 to 8000 by 500 of the squares
of the errors. Keep the exponent

** s2[A_, B_, P_] := Module[
{v = 0},
Do[v += errdiff2[s, f, A, B, P],
{s, -0.2, 0.2, 0.02}, {f, 1, 8000, 500}]
v]
Timing[Plot3D[s2[i, j, 2.375],
{i, 8, 12}, {j, 27, 37}, ViewPoint{-1.3,
2.4, 1.5}]]**

`{247.496 Second , - SurfaceGraphics -}`

Notice the timing: this plot took some time to produce, so we're having to wait a little for results. More importantly, notice the broad trough at the bottom, suggesting that a minimum can be found in the trough somewhere. A contour plot of the same region gives us even more confidence that we're on the correct archaeological site and can begin digging:

** Timing[ContourPlot[s2[i, j, 2.375], {i, 6, 14}, {j, 17,
47}]]**

`{256.328 Second , - ContourGraphics -}`

The contour plot shows that the trough is roughly elliptical and most likely
holds a minimum of . It is not *certain* to do, and this business of
parameter searching can be quite hazardous--it's at least as tricky as
integration (part 28). However, for this application, we're almost splitting
hairs. The "Magic Trick" function, by eyeball, looks like an effective
substitute for full Pajecka, for a wide range of (*A, B, P*)
parameters. We can probably get "close enough" just by eyeball. In
fact, the trough in the plot is quite broad, suggesting that a large range of
values produce almost equally good results. In any event, we should perform a
modicum of due diligence just in the name of professionalism of practice.

Look at another slice, this time fixing *B*.

** Timing[Plot3D[s2[i, 31, k],
{i, 9.5, 9.75}, {k, 2.25, 2.5}, ViewPoint{-1.3,
2.4, 1.5}]]**

`{262.708 Second , - SurfaceGraphics -}`

** Timing[ContourPlot[s2[i, 31, k], {i, 9.5, 9.75}, {k, 2.25,
2.5}]]**

`{263.458 Second , - ContourGraphics -}`

This slice has almost the same character as the previous one: a broad, narrow trough that likely contains the minimum of . We are now fortified with confidence and proceed to search:

** run[step_, {alo_, ahi_}, {blo_, bhi_}, {plo_, phi_}] :=
**

Do[

If[etest < ebest, ebest = etest

abest = ax bbest = bx pbest = px

Print[{ebest, abest, bbest, pbest}], Null],

{ax, alo, ahi, step},

{bx, blo, bhi, step},

{px, plo, phi, step}]

{ebest, abest, bbest, pbest}]]

run[1 / 8., {9.5, 9.75}, {30.50, 31.50}, {2.250, 2.5}]

` {4.08339 x 10`

^{7}, 9.5, 30.5, 2.25}

{2.90981 x 10^{7}, 9.5, 30.5, 2.375}

{2.90128 x 10^{7}, 9.625, 30.875, 2.375}

{2.89942 x 10^{7}, 9.625, 31., 2.375}

{82.208 Second, {2.89942 x 10^{7}, 9.625, 31., 2.375}}

Voila! We have found the promised values (*A, B, P*) = (9.625,
31.0, 2.375).

Notice how trivial is the searching code: just a triple loop over the
parameter space. However, the embedded function,

, is,
itself a double loop over the slip **s2***s* and vertical load *F _{z}*,
so we have a quintuple loop. I invite you to experiment with longer-running,
finer-grained searches by varying the first input to

**run**

.
The inputs are the search step size and the search boundaries for the three
parameters. All three parameters are stepped by the same increment. I warn you
that as the search grain is made finer, the minimum becomes slippery and the
boundaries of the search have to be adjusted. This is the trial-and-error that I
mentioned at the beginning of the article. It doesn't matter much that we have
an