See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/222848337

Exact kinematics analysis of Car’s suspension mechanisms using symbolic

computation and interval analysis

Article  in  Mechanism and Machine Theory · April 2005

DOI: 10.1016/j.mechmachtheory.2003.07.003

CITATIONS READS

40 7,482

3 authors:

Yves A. Papegay Jean-Pierre Merlet

National Institute for Research in Computer Science and Control National Institute for Research in Computer Science and Control Sophia-Antipolis
43 PUBLICATIONS   471 CITATIONS    298 PUBLICATIONS   11,048 CITATIONS   

SEE PROFILE SEE PROFILE

David Daney
National Institute for Research in Computer Science and Control
96 PUBLICATIONS   2,058 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Personally Assisted Living INRIA View project

Assistance robotics View project

All content following this page was uploaded by Yves A. Papegay on 23 November 2017.

The user has requested enhancement of the downloaded file.

 Exact Kinematics Analysis of Car’s Suspension
Mechanisms
Using Symbolic Computation and Interval Analysis

Yves A. Papegay – Jean-Pierre Merlet – David Daney

INRIA Sophia Antipolis, SAGA team
2004 route des lucioles, 06902 Sophia-Antipolis, France
[email protected] -- [email protected] --
[email protected]

Abstract
This paper shows how some exact computational methods can be (and
have been) used for the kinematic analysis of mechanisms involved in sus-
pension of vehicles. It points out that, among the tools recently devel-
opped to solve exactly sets of algebraic equations, some are powerful and
efficient enough to analyse and solve the kinematics of such mechanisms.
This paper presents numerical results which are much more accurate and
complete than the one obtained through numerical iterative computa-
tions. Those results are illustrated on two mechanisms, the simple but
very popular Mac Pherson mechanism and one of the most complicated
suspension mechanism involving 5 rods with 5 spherical joints.

1 Introduction
Multi-link suspensions are yet very popular on a lot of vehicle and the numerical
simulation of their motions is a critical step of their design process. Hence the
kinematics of these mechanisms are the main topics of a large number of studies
([?],[?],[?]) even if mostly the simplest case of the well known Mac Pherson
mechanism has been thoroughly investigated ([?],[?]).
Among the papers addressing the problem of the kinematics analysis of such
mechanisms, not one (except [?]) considers the problem of the reliability of
the Newton-Raphson iterative scheme that is usually used to solve the highly
non-linear equations appearing in the kinematics. Indeed, being given a set of
numerical values for the parameters of the mechanism, such a numerical iterative
process leads to a single result for the position of the wheel while the system of
equations usually have several real solutions corresponding to several different
positions. So, when such a process is used at each step of the computation of a
discrete approximation of the trajectory of the wheel during a motion, there is
no guarantee that the whole set of points defining the positions are belonging to

1
 Greyed curves are the two branches of
the solution curve, which correspond to
two different trajectories, black curve is
the computed trajectory.
Figure 1: what could happen

the same branch of the corresponding curve, especially as the Newton-Raphson

scheme is not guaranteed to converge toward the solution which is the closest
to the initial guess (a fact which is unknown to many design engineers). Hence,
the reliability of the computed trajectory is done by a visual check which can be
faulty if two trajectories came close to each other at some point of the process
- see figure 1.
When small differences on camber or toe angle or on orientation of steering
axis may have a drastical influence on the handling performances of a car, it is
natural to investigate if existing computational methods can give more precise
and complete qualitative results and/or more accurate quantitative results for
the kinematics analysis of suspension mechanisms. For this challenging goal,
exact computational methods are good candidates.
Exact computational method denotes here any kind of method for solving
equations which returns the whole set of real solutions, given with its accuracy.
So, it is a general terminology for symbolical methods based on algebraic proper-
ties of the equations, certified numerical methods based on analytical properties
of the equations and hybrid methods combining a symbolic pre-processing of
the equations and a certified numerical solving of the resulting problem.
In the next section, we will describe how we suggest to build the geometrical
model of multi-link mechanisms for suspensions for their kinematics analysis
and what is the algebraic structure of this model and we will review some
recent methods for computing exactly the solutions of systems of equations. In
the next section, we will describe more precisely the method based on interval
analysis, that we have extensively experimented. In the two remaining sections,
we will present detailed studies of the kinematics of two typical mechanisms –
Mac Pherson Strut and 5-rod 5SK mechanism – which respectively corresponds
to one of the simplest and one of the most complicated mechanisms and leads
to typical systems of equations.

2 Modelling Suspension Mechanisms

At an early stage of the design process, a suspension mechanism is modeled as
a multibody rigid system forming several closed kinematic chain in the three-
dimensional space.
Most of the joints between the bodies are physically elastokinematic joints
but are modeled as spherical (ball) joints. Revolute joint may also appear, espe-

2
cially after simplifications due to an analysis of the geometry of the mechanism.
Some advanced mechanisms also include cylindrical joints and more seldomly
universal joints.

2.1 Geometrical analysis

In the geometrical analysis, we express the constraints induced by the presence
of a joint between two bodies by writing equations that imply that a point has
to remain respectively on a sphere, a circle or a cylinder – from this point of
view, it is not necessary to distinguish between spherical and universal joint.
These equations involve lengths and three-dimensionnal position of selected
points as parameters and they are algebraic with respect to them as instances
of one of the following formulas,
• for a spherical or universal joint:

(xi − xj )2 + (yi − yj )2 + (zi − zj )2 = lij (1.1)

• for a revolute joint:

(xi − xj )2 + (yi − yj )2 + (zi − zj )2


= lij
(1.2)
(xi − xj )ui + (yi − yj )vi + (zi − zj )wi = kij

• for a cylindrical joint:

lij = (vi2 + wi2 )(xj − xi )2 + (u2i + wi2 )(yj − yi )2

+(u2i + vi2 )(zj − zi )2 − 2ui vi (xj − xi)(yj − yi ) (1.3)
−2ui wi (zj − zi)(xj − xi ) − 2vi wi (yj − yi)(zj − zi ))

where (xi , yi , zi ) and (xj , yj , zj ) denotes the coordinates of two points belonging
to each of the linked bodies, and where (ui , vi , wi ) and kij , li,j denotes the
coordinates of an unitary vector and other constants which are only depending
on the geometry of the mechanism.
This set of non-linear equations expressing the mechanical constraints is
known as geometrical model using natural coordinates.
Depending of the choice of the points involved in the equations, the model
may also contains equations expressing that several points are belonging to the
same rigid body:
• for 2 points (x1 , y1 , z1 ) and (x2 , y2 , z2 ), it expresses that the distance is
constant:
(x1 − x2 )2 + (y1 − y2 )2 + (z1 − z2 )2 = l12 (2.1)

• for 3 points, it expresses that the three distances are constants:

 (x1 − x2 )2 + (y1 − y2 )2 + (z1 − z2 )2 = l12



(x1 − x3 )2 + (y1 − y3 )2 + (z1 − z3 )2 = l13 (2.2)

(x2 − x3 )2 + (y2 − y3 )2 + (z2 − z3 )2 = l23


3
 • for any additional point (xi , yi , zi ), we write down that their coordinates
is a linear combinations of the coordinates of the three first one:

 xi = k1,i x1 + k2,i x2 + k3,i x3
yi = k1,i y1 + k2,i y2 + k3,i y3 (2.3)
zi = k1,i z1 + k2,i z2 + k3,i z3


2.2 Kinematics Analysis

Several kinematics analysis can be performed with the same geometrical model
according to which parameters are considered as known. If we address a design
problem, a simulation or a calibration problem, some parameters will be variable
and be therefore the unknowns of the problem while the other will be fixed.
That is why, once the geometrical analysis has been done, we have to separate
parameters in three classes, physical constants, input variables, and output
variables. The kinematics analysis consist in studying how to express output
parameters in terms of the constants and the input parameters.
The maximum number of independent input parameters, say DoF , is gen-
erally known as number of degrees of freedom. It is also the minimum number
of variable parameters that are necessary to fully describe the kinematics of
the mechanism. To determine the number DoF of a mechanism, the following
general formula1 is used:
p
X
DoF = d(n − 1) − (d − fi ) (3)
i=1

where d is the dimension of the configuration space - 6 in our case -, n and p

are respectively the number of bodies and of joints and fi is the local number
of degrees of freedom corresponding to the i-th joint (1 for a revolute, 2 for a
cylindrical and 3 for a spherical joint). This formula is only valid when bodies
and joints are in generic position and some ”always singular” mechanisms such
as the Goldberg or Bennet mechanism may require further analysis.
As most of the equations of the geometrical model are non-linear, there are
usually, for a given set of values of the input variables, several sets of values of
the output which are solutions of the system of equations.
The kinematics analysis, which aims at describing all the possible motions
for the mechanism, has two parts:
• Solving i.e to find all the sets of values corresponding to real solutions for
a given set of values of the input parameters. It has to be done a lot of
time for a wide range of values of these parameters.
• Following trajectories, by describing how each solution varies when the
input is changing continuously, whith a special attention to pay at singular
points, where this trajectories are crossing.
1 known as Grübler-Kutzbach criterion [?, ?]

4
3 Exact Computational Methods
Three important exact computational methods have been reviewed by Raghavan
and Roth in [?] in the scope of kinematics analysis problems, Grbner basis,
continuation and dialytic elimination. They are respectively representive of
three family of methods.

Symbolical Methods may only be used for solving set of polynomial ex-
pressions with rational or symbolic coefficients. They are computing explicit
symbolic expressions of the output parameters (in our case, the parameters
defining the position of the wheel) in terms of the input and fixed parameters
(in our case, the length of the spring and the design parameters of the suspen-
sion mechanism). Using symbolical methods is very efficient for solving, the
whole set of numerical result being simply obtained by numerical evaluation of
the symbolical expressions but is also very efficient for following trajectories, as
the information about each branch of a trajectory is attached to the computed
expressions. Symbolic methods like Grbner basis or Triangular set general-
ize the Gauss elimination to set of non-linear equations. It gives a theoretical
framework for the human way of solving equations ”by hand” by successive
eliminations of variables, repeating:
Select one equation, select one variable appearing in the selected
equation, express the selected variable in term of the other using the
selected equation and then, substitute this expression in the others
equations.
Unfortunately the number of monomials involved in intermediate expressions is
huge even on academic examples and is very quickly increasing with the number
of variables and the number of equations, as well as the size of the rational
coefficients. In the worst cases, the time and space complexity bound is a twice
exponential function of the number of variables, the number of equations and
the degree of the equations. We will show that it is not practicable, even for
the simple Mac Pherson mechanism, to compute the symbolical expression of
the output parameters using such a method.

Symbolic-numeric methods have been mostly developped for solving set

of polynomial equations, to overcome the unefficiency of Grbner basis based
methods while taking benefit of the algebraic structure of equations. Typically,
a pre-processing of these equations is done, computing a matrix from the coef-
ficients of the polynomial equations. Using this matrix, one can transform the
solving problem into a linear eigenvalues computation problem or into an uni-
variate polynomial solving problem. In both cases, the information about the
number of solutions attached to the matrix is used by the numerical procedure
which computes the whole set of solutions. Even if the intermediate matrix is
by construction very large, numerous efficient hybrid methods have been devel-
opped and implemented recently ([?]). But those method provide a numerical

5
result for each given set of input parameters and hence have no capabilities for
following trajectories.

Certified Numerical Methods may be used even if the equations are not
algebraic. They are based on numerical scheme with adaptative steps or on
homotopy algorithms. Combined with interval arithmetics, or with arithmetics
with a controlled precision, they can guarantee qualitative and quantitative in-
formation on the results or report numerical problems. Like hybrid methods,
certified numerical methods are of no help for following trajectories and special-
ized numerical procedure must be design for this purpose.

3.1 Solving with Interval Analysis

An alternative possibility to solve a system of equations is based on interval
analysis [?]. In this method each variable X is replaced by a range [x, x] and
each mathematical operator has a new definition. For example the ”+” operator
for two ranges X1 , X2 is defined as the range X12 such that:

X12 = [x1 + x2 , x1 + x2 ]

which means that the sum of any number in the ranges X1 , X2 is strictly included
in the range X12 . As the operators that have an interval definition are almost
all the classical mathematical functions, interval analysis enable to determine
lower and upper bounds for the value of any mathematical expression. Clearly
if 0 is not included between these bounds, then the value of the expression will
never be 0 whatever the values of the unknowns. A drawback is that usually
the bounds will be overestimated. For example if we consider the function
F = x2 + x when x ∈ [−1, 3] then:

F ([−1, 3] = [0, 9] + [−1, 3] = [−1, 12]

while the real bounds are [-1/4,12]. But an advantage of interval computation
is that it is possible to implement interval arithmetics in such way that round-
off errors are taken into account i.e. the interval evaluation of a function is
guaranteed to include the real bounds for the function (hence if 0 is not included
in the interval evaluation of an equation for given intervals for the unknowns
we may guarantee that there is no root of this equation within the intervals for
the unknowns).
A basic equation solver based on interval analysis will use a bisection process
that we will illustrate on the equation F in one unknown. If we are looking
for the solution of F (x) = 0 when x lie in the range [-1,3] we will start by
computing the interval evaluation of the function for this range. As this lead
to the range [-1,12] we cannot conclude if F has a real root in the range. The
bisection process consists in splitting the initial interval into two new intervals,
for example [-1,1],[1,3] and repeat the process on each new interval. In our case
the interval evaluation of F for the range [1,3] leads to [2,12] that does not
include 0. Hence this interval will be rejected and we will repeat the bisection

6
process only on the interval [-1,1]. The process will return as solution the range
whose width (the difference between the upper and lower bound) is lower than
a given threshold and for which the interval evaluation include 0.
One main drawback of interval analysis is that ranges for the unknowns
should be given and that only the roots within these ranges will be found.
For the kinematic analysis this is not a problem as the distance equations we
are dealing with enable one to determine easily ranges that are guaranteed to
contain all the roots.
There are numerous tricks to improve the efficiency of the basic solver.
Among them we may mention:
• using the monotonicity of the function. We may consider the derivatives
of the function and compute its interval evaluation: if the lower and upper
bound of a derivative with respect to a given variable have same sign
then the function is monotonous with respect to this variable and the
minimum and maximum of the function will be evaluated not by using
the range for the variable but by fixing its value according to the sign
of the derivative. Note that this must be a recursive process: fixing the
value of a variable may change the interval evaluation of the derivative
with respect to other variables. Using the monotonicity better interval
evaluation of the function will be obtained
• using Moore or Kantorovitch test [?], [?]. These methods give condi-
tions respectively on the gradient and on the gradient and Hessian of the
equations such that if they are fulfilled there will be a unique real root
within some ranges. Furthermore they guarantee that a numerical method
(respectively the Krawczyk operator and Newton scheme) will converge to-
ward this root. These tests enable one usually to avoid a large number of
bisection.
• using parallel computation: at some point of the process we have a list of
set of ranges that may contain a real root of the equations. Determining
if a set contain one root is not dependent upon the presence of roots in
any other set. Hence this set can be examined independently and interval
analysis is a good candidate for parallel computation
For solving the system of equation of the suspension mechanism we use two
others methods:
• the 3B method: let x1 , . . . xn be the unknowns in the equations. We
choose one unknown in this list, say xi whose range is [xi , xi ] and a small
number ǫ such that xi + ǫ < xi . We will substitute the range for xi by the
range [xi , xi +ǫ] while leaving the other ranges untouched and compute the
interval evaluation of the equations for this new set of ranges. If at least
one interval does not contain 0, then the range for xi may be reduced to
[xi + ǫ, xi ]. The process is then repeated until no improvement is obtained
for xi . We may also use this procedure for the other variables and for the
right side of the interval.

7
 • by using the simplex approach [?]: it may be seen that all the equations
that appears in the kinematics equations may be written as the sum of a
linear part in the unknowns and a non-linear part. Let substitute this non-
linear part by a new variable Yi . Using interval analysis we may compute
and interval evaluation of Yi = [Yi , Yi ]. Each equation is now linear in the
augmented set of unknowns and we have a set of inequalities constraints
for each unknowns. Hence we may use the simplex method to determine
if there is a feasible region for the unknowns and eventually to determine
the minimum and maximum values of the unknowns, thereby improving
their ranges. This method enable to deal with a major problem of interval
analysis: each equation is considered independently and while there may
be values for the unknowns in their range for which each equation is 0
there may not be values for which all the equations are 0 at the same
time
To test this solving method we use the ALIAS library, a C++ library 2 that
implement high-level solving procedures using the BIAS/Profil library3 for the
basic interval operations.

3.2 Following a trajectory

As we have seen there are various tools to determine all the possible configura-
tions of the mechanism for a given value of its degree of freedom l. But what
we are really interested in is to determine the configurations for a whole set of
value of the degree of freedom, say at l0 , l0 + ∆l, l0 + 2∆l, . . . l1 where ∆l is
arbitrary number, i.e. to get an approximation of the trajectories of the mecha-
nism as function of the degree of freedom. We may evidently use repeatedly the
solving procedures at each step but this will lead to a lengthy procedure and
furthermore problems may occur if two trajectories come close to each other as
we may have difficulties to determine to which branch belong a configuration.
Alternatively we may want to follow a given branch being given an approximate
value for the configuration of the mechanism.
The usual approach to both problems is to use the result of the solving
procedure at one step as the initial guess of a Newton scheme that will be used
to determine the configuration for the next step. But this is an unsafe method
as there is no guarantee that the Newton scheme will converge to the solution
that belong to the same branch. We propose to improve this method by using
Kantorovitch theorem. Let a system of n equations in n unknowns:
F = {Fi (x1 , . . . , xn ) = 0, i ∈ [1, n]}
Let x0 be a point, called P
the central point, and U = {x/||x − x0 || ≤ 2B0 }, the
norm being ||A|| = Maxi j |aij |. Assume that x0 is such that:
1. the Jacobian matrix of the system has an inverse Γ0 at x0 such that
||Γ0 || ≤ A0
2 www-sop.inria.fr/coprin/software/ALIAS/
3 www.ti3.tu-harburg.de/Software/PROFILEnglisch.html

8
 2. ||Γ0 F (x0 )|| ≤ 2B0
Pn ∂ 2 Fi (x
3. k=1 | ∂xj ∂xk | ≤ C for i, j = 1, . . . , n and x ∈ U

4. the constants A0 , B0 , C satisfy 2nA0 B0 C ≤ 1

Then there is an unique solution of F = 0 in U and Newton scheme used with
x0 as estimate of the solution will converge toward this solution [?].
Let us assume that we have computed the set of m solutions of the system
F , {X1 = (x11 , . . . x1n ), . . . Xm = (xm m
1 , . . . xn )}, for a given value l0 of the degree
of freedom l. Let dj be the minimal distance between the solutions set Xj and
all the other solutions. Assume now that the degree of freedom is changed to
l0 + ∆l. We use Kantorovitch theorem for each solution set Xj with as central
point the set itself. Two cases may occur when applying this theorem on the
set Xj :
1. the hypothesis of Kantorovitch theorem are satisfied and dj > 2B0
2. the hypothesis of Kantorovitch theorem are not satisfied or dj < 2B0
In the former case Newton scheme applied to the kinematic equations obtained
for l0 + ∆l with as initial guess Xj will converge to a solution that lie on the
same branch than Xj (otherwise there will two solutions within the ball, which
contradict Kantorovitch theorem). In the later case we are not sure that Newton
scheme will converge: hence we will consider the system of equations obtained
for l0 + ∆l/2 and repeat the procedure. Here also two cases may occur:
1. the process is successful for a given value of the degree of freedom
2. the process is not successful even for a very small change to l
In the former case we will apply the procedure to another branch and a syn-
chronisation mechanism will ensure that each branch is followed in sequence and
that the configuration solution of each branch are computed for each ∆l step. In
the later case Kantorovitch will fail because the central point of a branch is too
close to a configuration in which the jacobian of the system is degenerate, i.e.
we are close to a singularity of a mechanism and at least two branches collapse.
In that case the procedure will signal that a singularity has been encountered.
As we need a new starting point to follow the branches we will start again a
solving procedure on the system obtained when the degree of freedom is l0 + ∆l.
If the procedure is unable to find a Kantorovitch solution with this new value
we change it to l0 + 2∆l and so on until a valid set of solutions is found.
A drawback of this method is that we will follow only the branches that
have been detected when using the solving procedure: isolated branches that
may appear for a value of the degree of freedom strictly included between l0 and
l1 will not be detected. This is not a problem if we are interested in following
only one branch, as this is usually the case for the suspension mechanisms.

9
 Figure 2: MacPherson mechanism

4 Mac Pherson mechanism

The McPherson4 suspension mechanism - shown on figure 2 - is the most widely
employed front suspension mechanism for small cars. It is made of three main
components linking the vehicle body to the wheel spindle,
1. the lower arm, attached to the body through a revolute joint and to the
spindle through a spherical joint,
2. the tie rod, for steering control, which is a spherical-spherical crank,
3. and the strut of variable length rigidely fixed on the wheel side and fixed
with a ball joint on the body side – a spring-damper system connect the
ball joint and the wheel spindle.
From a kinematical point of view, this mechanism is hence made of 5 bodies and
has 3 kinematic chains with 2 loops. Using the Grübler-Kutzbach criterion, we
get that the number of degrees of freedom is 3 but this linkage has in fact only
one usefull degree of freedom for the wheel spindle as the 2 others are passive:
rod and struts links are rotating freely about their own axis. That is why it can
be denoted by RSCS-SS or equivalently by RS-SK-SP - in the Raghavan Atlas
[?].

4.1 Geometrical analysis

For the geometrical analysis, we assume the notations given in figure 3. We
define a point at each extremity of each link, and we denote by (xi , yi , zi ) the
4 from the name of its designer in the 1940s, Earl S. McPherson, General Motors

10
 Figure 3: MacPherson - notations

coordinates of point Ai (or Bi ) relative to a fixed reference frame attached to

the body of the car (O, Ox, Oy, Oz). Moreover we define a point B5 at the other
end of the coiled spring. We represent by lij the distance between point Ai (or
Bi ) and point Aj (or Bj ) and U0 is a unit vector of coordinates (u0 , v0 , w0 ) along
the rotation axis of the lower arm. We also need to introduce two constants, k32
and k42 describing respectively the angles (B5 B3 , B5 B2 ) and (B5 B4 , B5 B2 ).
As the mechanism has one degree of freedom, the purpose of the kinematic
analysis is to determine the pose of the wheel according to the value of one
variable parameter of the system. The pose of a body is fully defined if the
coordinates of three points on it are known in the reference frame – here, we
choose B3 , B4 and B5 . We select as the variable parameter the length of the
spring, l25 .
Using the geometrical model defined in section 2.1, expressing the constraints
for each link leads to the following system of equations:

(x3 − x0 )u0 + (y3 − y0 )v0 + (z3 − z0 )w0 = 0 (4.1)

(x3 − x0 )2 + (y3 − y0 )2 + (z3 − z0 )2 = l03 2 (4.2)
(x4 − x1 )2 + (y4 − y1 )2 + (z4 − z1 )2 = l14 2 (4.3)
2 2
(x4 − x2 )2 + (y4 − y2 )2 + (z4 − z2 )2 = l25 + 2k42 l25 l45 + l45 (4.4)
(x4 − x3 )2 + (y4 − y3 )2 + (z4 − z3 )2 = l34 2 (4.5)
2 2
(x3 − x2 )2 + (y3 − y2 )2 + (z3 − z2 )2 = l25 + 2k32 l25 l35 + l35 (4.6)
(x5 − x2 )2 + (y5 − y2 )2 + (z5 − z2 )2 = l25 2 (4.7)
(x5 − x3 )2 + (y5 − y3 )2 + (z5 − z3 )2 = l35 2 (4.8)
(x5 − x4 )2 + (y5 − y4 )2 + (z5 − z4 )2 = l45 2 (4.9)

As A’s points are fixed points, we only have as unknown variables the 9
coordinates of B3 , B4 and B5 . Other symbols correspond to fixed (but not
independent) parameters, except one of them taken as input, say for example

11
l25 expressing the degree of freedom.

4.2 Solving by Symbolical Methods

It is completely hopeless to enter this set of equations in a computer algebra
system like maple or mathematica and to try to solve it symbolically with the
help of the dedicated algebraic solving function to get the symbolic expressions of
the 9 variables in term of the parameters, especially without any pre-processing
of the equations.
Using a specialised software like gb or macaulay designed for solving very
efficiently algebraic systems of equations, we can have information on the set
of solutions. In this case, we get that the system is algebraic set of dimension
20 and degree 1024, which means that the values of 20 parameters may be
arbitrary fixed and once they are fixed there are at most 1024 set of solutions
for the 9 unknowns. Theoretically, we may also obtain a symbolic expression
of the solutions by computing a Grbner basis for the right elimination order
on the variables. But, due to the size of the problem and the high complexity
of the method – at least double exponential in the degree and the number of
the monomials – we did not succeed in getting the first result: macaulay fails
for memory size (exceeding 968 Megabytes) after having computed during 456
hours of cpu time on a dec alpha 1 with 633 Mb of ram.
The next challenge is to get the expressions of the 9 variables, not as a whole
symbolic formula, but as an evaluation program, i.e. a sequence of symbolic
formulas, each of them expressing a new parameter in terms of the parameters
previously expressed, for example

x3 = f1 (s, l25 )
y3 = f2 (s, l25 , x3 )
z3 = f3 (s, l25 , x3 , y3 )
z4 = f4 (t, l25 , x3 , y3 , z3 )
x4 = f5 (t, l25 , x3 , y3 , z3 , z4 )
y4 = f6 (t, l25 , x3 , y3 , z3 , z4 )
z5 = f7 (t, l25 , x3 , y3 , z3 , x4 , y4 , z4 )
x5 = f8 (t, l25 , x3 , y3 , z3 , x4 , y4 , z4 , z5 )
y5 = f9 (t, l25 , x3 , y3 , z3 , x4 , y4 , z4 , z5 )

where s and t denote respectively x0 , y0 , z0 , u0 , v0 , w0 , l03 , x2 , y2 , z2 , k32 , l35 and

s, x1 , y1 , z1 , l14 , k42 , l45 , l34 .
On our system, it can be achieved with a computer algebra system as maple
by using the following algorithm:
1. Solve equations (4.1) to get the expression of z3 in term of x3 and y3 ,
−u0 x3 + u0 x0 − v0 y3 + v0 y0 + z0 w0
z3 = (5.1)
w0

12
2. substitute the previous result into (4.2) to eliminate z3 from the equation
and solve it to get the expression of y3 in term of x3 , we get two expressions
√
2w0 2 y0 − 2u0 x3 v0 + 2v0 2 y0 + 2u0 x0 v0 ± 2 ∆
y3 = 1/2 (5.2)
v0 2 + w0 2
with

∆ = 2v0 2 w0 2 x3 x0 + 2w0 2 u0 2 x3 x0 − w0 4 x3 2 − w0 4 x0 2 + w0 4 l03 2

+2w0 4 x3 x0 − w0 2 u0 2 x3 2 − w0 2 u0 2 x0 2 − v0 2 w0 2 x3 2 − v0 2 w0 2 x0 2
+v0 2 l03 2 w0 2

Note that these expresions correspond to real values of y3 only when x3

is such as ∆ is non-negative.
3. substitute (5.1) and one (resp. the other) of the two previous results (5.2)
into (4.6) to eliminate both y3 and z3 and obtain a quadratic equation
governing x3 , solve it (in both cases) to get (in each case) two results
expressing x3 in terms of the fixed parameters and l25 . We obtain four
results on the same pattern:
p
4 Pi (t) + Qi (t)
x3 =
Ri (t)

where Ri (t) is a polynomial with 30 monomials of degree 4, Qi (t) is a

polynomial with 106 monomials of degree 5 or 6 and Pi (t) is a polynomial
with 2047 monomials of degree 10, 11 or 12. They are as many candidates
for the f1 expression of the program, but only two at most are valid one.
4. Compute the differences between equation (4.5) and equation (4.4) and
between equation (4.5) and equation (4.3). Quadratic terms in x4 , y4 or
z4 are vanishing and we get two linear equations in those variables:

(−2x3 + 2x2 )x4 + (−2y3 + 2y2 )y4 + (−2z3 + 2z2 )z4

−y2 2 + x3 2 − x2 2 + y3 2 + z3 2 − z2 2
= l34 2 − l25 2 − 2k42 l25 l45 − l45 2
(−2x3 + 2x1 )x4 + (−2y3 + 2y1 )y4 + (−2z3 + 2z1 )z4
−y1 2 + x3 2 − x1 2 + y3 2 + z3 2 − z1 2
= l34 2 − l14 2

By solving them, express x4 and y4 in terms of z4 :

−2Dy4 = (−2x1 z3 + 2x3 z1 − 2x3 z2 + 2x1 z2 − 2x2 z1 + 2x2 z3 )z4 − x3 l25 2

−x3 l45 2 − 2x3 k42 l25 l45 + 2x1 k42 l25 l45 + x3 l14 2 + x3 z2 2 + x3 x2 2
+x3 y2 2 − x1 y2 2 − x1 z2 2 − x1 l34 2 + x1 l25 2 + x1 l45 2 − x3 x1 2
+x1 z3 2 − x1 x2 2 + x2 y1 2 + x2 x1 2 − x2 z3 2 − x2 y3 2 − x2 x3 2

13
 −x2 l14 2 − x3 z1 2 − x3 y1 2 + x1 y3 2 + x2 l34 2 + x2 z1 2 + x1 x3 2
2Dx4 = (−2z1 y2 + 2z3 y2 − 2y3 z2 + 2y1 z2 + 2z1 y3 − 2y1 z3 )z4
−2y3 k42 l25 l45 + 2y1 k42 l25 l45 + y3 z2 2 + y3 x2 2 − z3 2 y2 + l34 2 y2
−y1 2 y3 + y1 2 y2 − z1 2 y3 − l14 2 y2 − x3 2 y2 + y1 y3 2 + y1 x3 2
+l14 2 y3 − y1 x2 2 − y1 y2 2 − y1 z2 2 − y1 l34 2 + y1 l25 2 + y1 l45 2
+y1 z3 2 + y3 y2 2 + x1 2 y2 − x1 2 y3 + z1 2 y2 − y3 2 y2 − y3 l45 2 − y3 l25 2

with
D = −x1 y3 − x2 y1 + x2 y3 + x3 y1 − x3 y2 + x1 y2

5. substitute the two previous expressions in equations (4.3) to eliminate

both x4 and y4 , solve it to get the expression of z4 in terms the fixed
parameters, l2 5 and x3 , y3 , z3 . We obtain two results on the same pattern:
p
4 Pi (t, x3 , y3 , z3 ) + Qi (t, x3 , y3 , z3 )
z4 =
Ri (t, x3 , y3 , z3 )
where Ri (t, x3 , y3 , z3 ) is a polynomial with 63 monomials of degree 4,
where Qi (t, x3 , y3 , z3 ) is a polynomial with 305 monomials of degree 5
or 6 and Pi (t, x3 , y3 , z3 ) is a polynomial with 14275 monomials of degree
10, 11 or 12.
6. execute analoguous steps to step 4 and 5, replacing respectively equations
(4.3), (4.4) and (4.5) by (4.7), (4.8) and (4.9) and the variables x4 , y4 , z4
by x5 , y5 , z5 .
The implementation of this algorithm compute all the expressions appearing
in the complete program in less then 2 minutes.
And, if we give numerical values to the parameters – a realistic set of data is
described below – we can use this evaluation program for computing numerically
the solutions. This can be done in around 1 minute by maple, but such a
performance may be largely enhanced by generating a C program and using it
for the numerical computations.

 l03 = 317.4905466 
 u0 = 0.9901559  x1 = 106


  l14 = 306.8979231
 
v = 0.09522448  y1 = 290
  
0
  

  l34 = 180.2625676
 
w0 = −0.1025845607 z1 = 90
  
l35 = 126.3992495

 x 0 = 33.69279  k42 = −0.02052180122 
  2 = 10
 x
y = 342.8049  y2 = 510
 
0
  

  k43 = 0.86574862
 
z0 = −0.03617776374  z2 = 583
 
 
l45 = 134.0027263

The figure 4 and 5 presents respectively a 3-dimensional view of the set of

acceptable positions for B4 and for B5 when the parameters l25 varies between
811 and 231. The figure 6 shows the two curves obtained by connecting the
points of the set of acceptable positions for B5 , which are two separate supports
of any trajectory of B5 .

14
Figure 4: acceptable positions for B4

15
 Figure 5: acceptable positions for B5

Figure 6: acceptable positions for B5 – the two branches

16
 A5

B5 B4
A4

B3 chassis
A3

B2

B1 A2

wheel
A1

Figure 7: 5-rod 5SK mechanism

4.3 Solving with ALIAS

Using the numerical data presented in the previous section and the algorithm
presented in section 3.1 it is possible to obtain similar plots. The computation
time for determining all the branches is about 15 minutes on a Pentium III at
600 Mhz. If only one branch has to be followed the computation time is about
1mn 40s. Note that this time may be drastically reduced if we use the specific
solver presented in section 5.2.1: in that case the computation time for finding
the 8 initial starting points of the branches is about 22 seconds while following
the 8 branches takes about 5 seconds (0.6 second if only one branch is followed).

5 5-rod 5SK
The mechanism, shown in figure 7, which is denoted by 5SK in the standard
classification has been used for the design of the suspension of the rear wheels
for more then 15 years on the Mercedes 190 serie. The wheel support is linked
to the body of the car through 5 rigid rods of fixed length. On both sides,
the links between the rods and the other parts are realised by elasto-kinematic
joints, modeled as ball joints. The spring-damper system is connected to one of
the 5 rigid rods and on the other side to the body of the car through ball joints.
From a kinematical point of view, this mechanism is highly redundant, it is
made of 8 bodies with 12 joints building 6 kinematics chains and 11 loops.

5.1 Geometrical analysis

The geometrical analysis is performed automatically as described in the section
2.1. The points on the body of the car are denoted by Ai = (xi , yi , zi ) for
i = 1 . . . 6 and the points on the wheel support by Bi = (ui , vi , wi ) for i = 1 . . . 5

17
assuming that for each i, Ai is connected to Bi through the i-th rod. A6 is
the point where is attached on the body the spring-damper system and B6 =
(u6 , v6 , w6 ) is at the other end of this system, on the first rod. Let l6 denotes
the length of the spring damper system while the other parameters are physical
constants describing the geometry of the mechanical system.
With these notations, constraints on the length of the rods can be writ-
ten:

(u1 − x1 )2 + (v1 − y1 )2 + (w1 − z1 )2 = l1 2 (6.1)

2 2 2 2
(u2 − x2 ) + (v2 − y2 ) + (w2 − z2 ) = l2 (6.2)
(u3 − x3 )2 + (v3 − y3 )2 + (w3 − z3 )2 = l3 2 (6.3)
(u4 − x4 )2 + (v4 − y4 )2 + (w4 − z4 )2 = l4 2 (6.4)
2 2 2 2
(u5 − x5 ) + (v5 − y5 ) + (w5 − z5 ) = l5 (6.5)
(u6 − x6 )2 + (v6 − y6 )2 + (w6 − z6 )2 = l6 2 (6.6)

We express the fact that the spring-damper system is attached to the first rod
by:

u6 = k6 u1 + (1 − k6 )x1 ) (6.7)
v6 = k6 v1 + (1 − k6 )y1 ) (6.8)
w6 = k6 w1 + (1 − k6 )z1 ) (6.9)

And we express that B1 , B2 , B3 , B4 and B5 belongs to the same rigid body by:

(u2 − u1 )2 + (v2 − v1 )2 + (w2 − w1 )2 = l12 (6.10)

2 2 2
(u3 − u1 ) + (v3 − v1 ) + (w3 − w1 ) = l13 (6.11)
(u3 − u2 )2 + (v3 − v2 )2 + (w3 − w2 )2 = l23 (6.12)

and by:

x4 = k14 x1 + k24 x2 + k34 x3 (6.13)

y4 = k14 y1 + k24 y2 + k34 y3 (6.14)
z4 = k14 z1 + k24 z2 + k34 z3 (6.15)
x5 = k15 x1 + k25 x2 + k35 x3 (6.16)
y5 = k15 y1 + k25 y2 + k35 y3 (6.17)
z5 = k15 z1 + k25 z2 + k35 z3 (6.18)

Alltogether, these 18 equations build a redondant geometrical model of this

mechanism governing the 18 coordinates of the points Bi .
As only one rigid body may be moving in this mechanism, we know that
theoretically, only 6 parameters are required to fully express its geometry. So,
it is natural to try to reduce this model to a simpler one before going to the
solving step.

18
 By substitution of the expressions of the coordinates of B4 , B5 and B6 we
get a reduced model with 9 equations and 9 unknowns:

(u1 − x1 )2 + (v1 − y1 )2 + (w1 − z1 )2 = l1 2 (7.1)

(u2 − x2 )2 + (v2 − y2 )2 + (w2 − z2 )2 = l2 2 (7.2)
(u3 − x3 )2 + (v3 − y3 )2 + (w3 − z3 )2 = l3 2 (7.3)
(u4 − (k14 x1 + k24 x2 + k34 x3 )) + (v4 − (k14 y1 + k24 y2 + k34 y3 ))2 +
2

(w4 − (k14 z1 + k24 z2 + k34 z3 ))2 = l4 2 (7.4)

(u5 − (k15 x1 + k25 x2 + k35 x3 ))2 + (v5 − (k15 y1 + k25 y2 + k35 y3 ))2 +
(w5 − (k15 z1 + k25 z2 + k35 z3 ))2 = l5 2 (7.5)
(u6 − (k6 u1 + (1 − k6 )x1 )))2 + (v6 − (k6 v1 + (1 − k6 )y1 )))2 +
(w6 − (k6 w1 + (1 − k6 )z1 )))2 = l6 2 (7.6)
2 2 2 2
(u2 − u1 ) + (v2 − v1 ) + (w2 − w1 ) = l12 (7.7)
(u3 − u1 )2 + (v3 − v1 )2 + (w3 − w1 )2 = 2
l13 (7.8)
(u3 − u2 )2 + (v3 − v2 )2 + (w3 − w2 )2 = 2
l23 (7.9)

5.2 Solving the kinematics analysis

We have made many attempt to get the symbolical expression of the Bi points
in term of the constants, the positions of the Ai and the length l6 of the spring-
damper system. But, due to the huge size of the expression generated by the
resolution, we didn’t succeed even allowing a maple process to use 10 Gigabytes
of memory. We failed also in trying to get an evaluation program for computing
the position of the Bi , as we did with the Mac Pherson mechanism.
In order to reduce the complexity of the computation, we gave explicit nu-
merical values to some of the parameters using the following realistic set of
data.

19
  

 x1 = −102.6 
 u1 = −48.1
y1 = 220.1 v1 = 660.5

  

l1 = 444.3865097

 

z1 = 238.1 w1 = 214.5

 
 

l2 = 313.9128701

 
 

x2 = 309.8 u2 = 34.6

 
 

l3 = 245.2243258

 
 

y2 = 506 v2 = 631.5

 
 

l4 = 246.4982353

 
 

z2 = 253.5 w2 = 169.5

 
 

l5 = 297.6022345

 
 

x3 = 201.9 u3 = 137.9

 
 

l12 = 98.51543026

 
 

y3 = 426.1 v3 = 660.5

 
 

 l13 = 191.7244116

 
 

z3 = 294.1 w3 = 261
 
l23 = 141.0111343
x4 = 206.7 u4 = 75.5
k14 = 1.706216400

 
 

y4 = 443 v4 = 651.5

 
 

k24 = −2.641106193

 
 

z4 = 398.2 w4 = 389.5

 
 

k34 = 1.805302996

 
 

x5 = −2.1 u5 = 5.1

 
 

k15 = 2.559324080

 
 

y5 = 363 v5 = 660.5

 
 

k25 = −3.529570293

 
 

z5 = 421.5 w5 = 424.5

 
 

k35 = 1.815276435

 
 

x6 = −112.5 u6 = −75.35

 
 

k6 = 0.5

  

y6 = 438.0 v6 = 440.3

 


 

z6 = 417.5 w6 = 226.3
 

Still with this numerical data we were not able to solve the problem using only
symbolic computation. In that case however we may still solve the problem
using the interval analysis approach. As mentioned in section 3.1 the first step
of the approach is to determine the starting points of the branches by solving
the system of equations for the initial value of the parameter. For the 5SK
mechanism this is a quite difficult problem, which is equivalent to solve the
direct kinematics of a parallel robot, a problem that is recognized to be one of
the most challenging in modern kinematics.

5.2.1 The distance solver

To solve this problem we use a specific solver of ALIAS called the distance solver.
If Xi is the 3D point of coordinates [ui , vi , wi ] we notice that all the equations
appearing in the system express that the distance from one Xi to either a fixed
point or another Xj is constant. We may hence specialize the methods proposed
in section 3.1 to improve their efficiency.
Let us give three examples of specific methods that improve the efficiency of
the distance solver, compared to the generic solver implemented in ALIAS. The
equations we have to solve may be written as:
j=3
X
Fi = (Xkj − Xm
j 2
) − L2i = 0 (8)
j=1

where Xkj is the unknown k-th coordinate of the point Xk and Xm j

is the k-th
coordinate of point Xm which may be either a fixed point or an unknown point.

20
A first example of improvment of the general solving method is to show that
it is not necessary to compute the gradient (which may be computer intensive)
to improve the interval evaluation of the Fi . Indeed, as the unknowns in this
equation appears only once in the Fi , a fundamental theorem of interval analysis
states that the bounds provided by the interval evaluation of Fi are exact, i.e.
there are values of the unknowns in the interval such that the value of Fi is
exactly one of the bounds given by the interval evaluation.
Compared to the general interval solving algorithm we may use here a second
improvement by using the so-called 2B method. Let us assume that Xm is a
fixed point with coordinates am , bm , cm . Equation (8) way be written as:
(Xk1 − am )2 = L2i − (Xk2 − bm )2 − (Xk3 − cm )2 (9)
Let U be the interval evaluation of the right side of this equation for the current
interval values of the unknowns Xk2 , Xk3 and let U be the upper bound of U . If
U is negative, then there are no root for this equation as the left side term of
the equation should be positive; hence we may state that there is no solution to
the system of equations for the current values of the intervals for the unknowns.
If we assume now √ that U is positive
√ then we get that Xk1 must belong to the
interval V1 = [am − U , am + U]. As Xk1 is an unknown it has already an
interval value V2 and Xk1 must be compatible with both the intervals V1 , V2 , i.e.
it must belong to V1 ∩ V2 . This operation may lead to a decrease in the width of
the interval for Xk1 . Evidently this procedure has to be repeated for Xk2 and Xk3
and for each equation in the system. Furthermore it has to be restarted as soon
as the interval for one unknown is modified as this change may have an effect for
an equation that has already been processed (however the convergence of this
method is asymptotically slow and it is better, in general, to stop this procedure
if the decrease of width of the interval is small as performing a bissection will
be more efficient). Note also that this procedure may be used to determine
efficiently initial values for the intervals for each unknowns.
The third improvement is related to Kantorovitch: for the specific case of
distance equations it is possible to show that a ”stronger” version of this theorem
may be demonstrated in the sense that a larger convergence domain will be
obtained. Using this version in the distance solver allows to determine that a
unique solution exists within larger interval than with the general version of
Kantorovitch test. This has also a drastic effect on the trajectory following
algorithm, as may be seen on the numerical example in the next section.
Note also that the distance solver may also be used to deal with the McPher-
son mechanism providing that the two first equations of the system are sub-
stituted by two equations expressing the point B3 (x3 , y3 , z3 ) is at a constant
distance from any two arbitrary points lying on the rotation axis of the joint of
the lower arm.

5.2.2 Numerical example

Using the distance solver the computation time for determining the initial start-
ing points for the branches is low. For example if we use the numerical data

21
 –50
0
50

u1
100
150
200
250
600
500
v1
400
300
600

500

400

300

200

100
w1
Figure 8: The 6 trajectories of the point of coordinates (u1 , v1 , w1 ) when l6 is
in the interval [173,280]

0
20
40
60
80
120
u2
160
200
240
700
600
500
v2
400
300
500

400

300

200

100
w2

Figure 9: The 6 trajectories of the point of coordinates (u2 , v2 , w2 ) when l6 is

in the interval [173,280]

presented in the previous section the 6 initial points for l6 = 280 are determined
in about 50 seconds.
Then, using the trajectory following algorithm, we are able to draw in an
exact way the 6 trajectories obtained when l6 varies between 174 and 280 in
about 10 seconds. These trajectories are represented in figure 8, 9, 10.
Clearly if we are interested in following only one branch the computation
time will be lower. Using the same data the computation time for following the
first branch is only 1.5 seconds.

6 Conclusion
We have presented two methods to solve exactly the kinematics of suspension
mechanisms
• a symbolic approach that has the advantages to provide the kinematics
solution directly as symbolic functions of the design parameters and is
hence appropriate for design purposes but has the drawback that not all
suspension mechanisms may be solved using this method,

22
 0
50
100
u3
150
200
250
650
600
550
500
v3
450
400
500

400

w3300

200

100
350
Figure 10: The 6 trajectories of the point of coordinates (u3 , v3 , w3 ) when l6 is
in the interval [173,280]

• an interval-based approach that allows to solve the kinematics of all the

suspension mechanisms but only for fixed values of their design parame-
ters.
Both methods allow to detect singularities in the mechanism but the symbolic
approach has the advantage that it allows to determine what parameters has to
be changed to avoidthe singularity.
It must be noted that both methods are exact in the sense that the computed
trajectories are guaranteed. This is a large improvement as undetected jumps
between branches may occur when using the usual method based on an iterative
use of Newton scheme.

23

View publication stats