Synthesis Method for the Spherical 4R Mechanism with Minimum
Center of Mass Acceleration
a

O. Mendoza-Trejo∗,

a

Carlos A. Cruz-Villar†, b R. Peón-Escalante‡, b M. A. Zambrano-Arjona§, b F. Peñuñuri¶

Cinvestav-IPN, Departamento de Ingenierı́a Eléctrica, Av. IPN 2508, A. P. 14-740, 07300, México D.F., México.
b

arXiv:1412.6850v2 [] 5 Apr 2015

a

Facultad de Ingenierı́a, Universidad Autónoma de Yucatán, A.P. 150, Cordemex, Mérida, Yucatán, México.

Abstract
In the mechanisms area, minimization of the magnitude of the acceleration of the center
of mass (ACoM) implies shaking force balancing. This article shows an efficient optimum
synthesis method for minimum acceleration of the center of mass of a spherical 4R mechanism
by using dual functions, as well as, the counterweights balancing method. Once the dual
function for ACoM has been formulated, one can minimize the shaking forces from a kinematic
point of view. We present the synthesis of a spherical 4R mechanism for the case of a path
generation task. The synthesis process involves the optimization of two objective functions,
this multiobjective problem is solved by using the weighted sum method implemented in the
evolutionary algorithm known as Differential Evolution. Our results shows that the magnitude
of the ACoM can be reduced about 20 times. Which for the presented synthesis, is traduced
in a reduction of about 89% for the shacking forces.
Keywords— Center of Mass Acceleration, Shaking Force, Spherical 4R Mechanism, Dual
Numbers, Differential Evolution.

1

Introduction

A mechanism is said to be dynamically balanced if there is an elimination of dynamic reaction
forces and torques (moments) on the base [1–4]. The use of counterweights [5–8] and springs [9–11]
are the two most common methods for balancing a mechanism.
The balancing of mechanisms is a very active field in the mechanisms area, and there is a big
amount of bibliography on such a topic, we refer the reader to works [12,13] where the balancing of
mechanisms is reviewed. Specifically regarding to the concept of the center of mass, it is worthwhile
to mention [14–16]. In [14], the balancing condition of parallel robots is obtained by expressing the
position vector of the center of mass as a function of the position and orientation of the platform,
and of six actuated prismatic joints. In [15], a synthesis method for linkages with the center of
mass at an invariant link point is presented. In [16], the shaking force minimization of high-speed
robots is done via optimal control of the acceleration of the total mass center of moving links.
∗

omendoza@cinvestav.mx
cacruz@cinvestav.mx
‡
rpeon@uady.mx
§
miguel.zambrano@uady.mx
¶
francisco.pa@uady.mx
†

1

The problem of cancellation or minimization of the total external force can be addressed from
a pure kinematic point of view by minimizing the ACoM. The idea of an optimization approach
has been successfully applied in [17–21]. The smaller the ACoM, the smaller the shaking force in
the mechanism. In the present work we show how, by using dual numbers [22–24] and the counterweights balancing method, a simple and computationally efficient optimum synthesis method
with minimum center of mass acceleration can be achieved. First, we provide a method where the
acceleration of the center of mass can be precisely and efficiently calculated. Second, since in some
cases, a mass redistribution can become equal to zero the acceleration of the center of mass, we
use an optimization based methodology, to find such a mass redistribution.
Our study is focused on the path generation task for a spherical 4R mechanism. The problem
involves two objective functions to be optimized, namely the objective function to fulfill the desired
trajectory and the objective function related to the acceleration of the center of mass. This
multiobjective problem is solved by using the weighted sum method, which is implemented on
the evolutionary algorithm known as differential evolution (DE) [25]. The synthesized mechanism
was constructed, and actuated by a dc motor controlled by a classical PID law. This allow
us to compare the required power (and energy used) to operate the balanced and unbalanced
mechanisms. The experimental measurement of the energy consumption was contrasted with its
theoretical calculation which is conducted by using the work–energy theorem, providing in this
way an experimental validation for the reduction of the ACoM in the balanced mechanism.
The rest of the paper is organized as follows, section 2 presents the essentials of our method.
There, the balancing method, the dual number approach to obtain derivatives and the optimization
method are presented. In section 3 a Spherical 4R mechanism is synthesized for a prescribed path
generation task with minimum center of mass acceleration. This section shows all the parameters
and formulas that need to be written in its dual form in order to obtain the acceleration of the center
of mass vector. Also a theoretical calculation of the energy consumption is presented. Section 4
shows the results along some discussions. Finally, the conclusions are presented in section 6.

2

Synthesis Method

In general, a mechanism has an arbitrary time-varying ACoM. So our aim is to minimize the
magnitude of the ACoM. For that purpose we need a precise and efficient method to obtain first
and second derivatives, as well as an optimization method that allow us to cancel or minimize the
ACoM. This section shows the implementation of the dual numbers to obtain derivatives and the
method used for balancing.

2.1

Minimum Acceleration of Center of Mass and Counterweight Method

There are several methods for balancing mechanisms, each one having advantages and disadvantages [3, 12, 13]. Since a mathematical point of view, more free parameters are introduced to the
synthesis problem to balance a mechanism. We implement this approach for balancing a spherical 4R (Fig. 1) mechanism by adding counterweights [5–7, 14] to the links as shown in Fig. 2,
therefore we have more independent variables in order to minimize the ACoM of the Spherical 4R
mechanism.

2

2.2

Dual numbers approach

Usually, the process of calculating a derivative is not difficult. However, for the case of a spherical
mechanism, obtaining first and second derivatives of the position vector for the center of mass is
not simple. Even when such derivatives can be explicitly obtained, the resulting expressions could
be of great complexity and useless for practical purposes. An alternative solution is to numerically
calculate such derivatives. Nevertheless, traditional methods for calculating numerical derivatives
(finite-differences) are not efficient enough to be used in the optimum synthesis of mechanisms.
Moreover, they are subject to both truncation and subtractive cancellation errors. This is an issue
to be dealt with in the optimum synthesis of mechanisms, since in some cases the errors introduced
by the numerical derivative method are of the same order than the quantity to be optimized.
A different approach that is not subject to the above mentioned errors can be constructed
by using dual numbers [26–28]. In order to make this paper self contained we briefly review the
essential ideas following [28] and bearing in mind a numerical implementation.
Dual numbers and derivatives
A dual number r̂ is a number of the form
r̂ = a +  b,

(1)

where a (the real coefficient) and b (the dual coefficient) are real numbers and 2 = 0. As in the
case of complex numbers, a dual number can also be defined as ordered pairs
r̂ = [a, b].

(2)

The algebraic rules for dual numbers can be found elsewhere in the literature, see for example
[22–24].
For a function f : R → R which accepts a Taylor expansion we have:
f 00 (x) 2
h + . . . + O(h3 ).
(3)
2
Now, considering the particular dual number x̂ = x+, i.e., a dual number where the coefficient
of the nilpotent  is equal to one we have fˆ(x̂) = f (x) + f 0 (x). In the notation of (2) we may write
f (x + h) = f (x) + f 0 (x)h +



fˆ(x̂) = f (x), f 0 (x) .

(4)

So, if instead of calculating over the reals, we calculate over the duals, we end up with a dual
function where the real part is the function itself and the dual part its derivative.
The extension to obtain the second derivative is straightforward. For this purpose we promote
the f 0 (x) function to be a dual function, that is:


fˆ0 (x̂) = f 0 (x), f 00 (x) .
(5)
The information of the function itself and its first and second derivative can be stored in a
vector of three components. We will use the notation


f˜(x) = f (x), f 0 (x), f 00 (x)
(6)
to represent the extended dual version of the original function f .
The following pseudo-code implements the extended dual version f˜ = [f0 , f1 , f2 ] of the function
f , as a function of the dual variable g̃ = [g0 , g1 , g2 ]:
3

start the definition of the dual function
declare the type of variables (if necessary)
dualf = [f_0(g_0),f_1(g_0)*g_1, f_2(g_0)*g_1^2 + f_1(g_0)*g_2 ]
end the definition of the dual function.
Following the above discussion, one can obtain the extended dual version of the center of
mass vector from where its velocity and acceleration are directly obtained. The derivatives are
calculated to the same precision of the implemented standard functions of the used programming
language, and its calculation is reduced to a simple function evaluation without the need of a
limit calculation process. As an added value to our work, some useful functions common in the
rotational kinematics are dualized for the first time and provided as additional material to this
article1 . So they will be of great help to any one interested in applying the dual numbers to obtain
derivatives. Since we are interested in the optimal synthesis of mechanisms (fast and efficient
programs are desired), we have coded the dual functions in Fortran but a translation to a more
user-friendly language is almost direct.

2.3

Optimization method

The optimal dimensional synthesis is conducted by using the evolutionary DE method [25]. This
method has been successfully applied to the optimal dimensional synthesis of mechanisms [29–32]
and mainly consists of three operators: mutation, crossover and selection; applied to a possible
solution or individual x. Below an outline of the DE method, following [25] is presented.
2.3.1

Differential Evolution Algorithm

1. The population is described by:
Px,g = (xi,g ), i = 1, ..., m; g = 0, ..., gmax
xi;g = (xji;g ), j = 1, ..., D;

(7)

where D represents the dimension of x, m represents the number of individuals, and g is the
generation.
2. Initialization of population.
xji;0 = randj (0, 1) · (bjU − bjL ) + bjL .
Vectors bU and bL are the limits of parameters. randj (0, 1) represents uniformly distributed
random number in [0, 1). Superscript j means that a random value for each parameter is
generated.
3. Mutation.
vi;g = xr0 ;g + F · (xr1 ;g − xr2 ;g ).

(8)

The main difference between DE and other evolutionary algorithms like genetic algorithms
is due to this mutation operator. xr0 ;g is called the base vector which is perturbed by the
difference of other two vectors. r0 , r1 , r2 ∈ {1, 2, ...m}, r1 6= r2 6= r3 6= i. F , the mutation
scale factor, is a real number greater than zero.
1

A file (Dual ACoM files.zip) containing the code of such functions can be downloaded from http://www.meca.
cinvestav.mx/personal/cacruz/archivos-ccv/

4

4. Crossover.
It uses a dual recombination of vectors to generate the trial vector:
 j
vi;g if randj (0, 1) 6 Cr or j = jrand
j
ui;g = ui;g =
 j
xi;g otherwise.

(9)

The crossover probability, Cr ∈ [0, 1], is a user-defined value, jrand ∈ [1, D].
5. Selection.
The selection is made according to

xi;g+1 =
2.3.2

ui;g if f (ui;g ) 6 f (xi;g )
xi;g otherwise

(10)

Objective function for the prescribed path generation task

In this study we are interested in the prescribed path generation task. So, besides an objective
function related to the minimization of the ACoM, we also consider an objective function intended
to fulfill the path generation task.
In the path generation problem, n points are given and a mechanism passing through these
desired points is required. In order to synthesize such a mechanism (in this study a spherical 4R
mechanism), we minimize the mean structural error (root mean square error) which is defined as:
n

fpath (θ1 , φ1 , η1 , ψ, β, γ, α1 , α2 , α3 , α4 ) =

1X
krdi − rgeni k2
n

!1/2
,

(11)

i=1

where rdi − rgeni is the difference between the i-th desired point rdi and the i-th generated point
rgeni . The θ1 parameter, is the angle for the input link corresponding to the first point. In the
prescribed path generation task the angle corresponding to the other points are given.
2.3.3

Objective function for minimum center of mass acceleration

In order to obtain a mechanism with minimum ACoM we minimize the objective function
n

fcm (θ1 , φ1 , η1 , ψ, β, γ, α1 , α2 , α3 , α4 , e1 , e2 , e3 , e4 ) =

1X
kr̈cm;i k2
n

!1/2
(12)

i=1

where r̈cm;i is the acceleration of the center of mass of the mechanism when it is passing through
the i-th point.

3

Synthesis of a Spherical 4R Mechanism

This section presents the application of the method above explained, to the synthesis of a spherical
4R mechanism for a path generation task, with minimum ACoM. Table 1 (as presented in [33])
shows the desired points. This problem was originally studied in [34] but there, the desired points
are not on a unit sphere.

5

Figure 1: Spherical 4R Mechanism and variables used to obtain rgen .

3.1

Kinematics

A spherical four-bar mechanism (Fig. 1) is a closed chain, consisting of four links and four revolute
joints. These linkages have the property that every link in the system rotates about the same fixed
point [35]. Thus, in a spherical mechanism, any point in a moving body is confined to move within
a spherical surface, and all spherical surfaces of motion are concentric. In order to find the coupler
point curve rgen , we will use the parameters shown in Fig. 1. Such equation has been already
obtained in [36], however, we want to change to a more concise description by using the lengths
of the links as parameters instead of the coordinates for the joints. Thus, the coupler-point curve
can be written as:
rgen (θ, φ1 , η1 , ψ, β, γ, α1 , α2 , α3 , α4 ) = R (π/2, rcp (β)) rcp (β + γ),

(13)

where φ1 and η1 are the azimuth and polar angle for the point x1 , R(ν, ŵ) is the active rotation
matrix of angle ν in direction of the unit vector ŵ.
rcp (ν) = R(ν, n̂23 ) r2 (θ, ψ, φ1 , η1 , α1 , α4 ),

n̂23 =

r2 (θ, ψ, φ1 , η1 , α1 , α4 ) × r3 (φ(θ), ψ, φ1 , η1 , α1 , α2 , α3 , α4 )
,
kr2 (θ, ψ, φ1 , η1 , α1 , α4 ) × r3 (φ(θ), ψ, φ1 , η1 , α1 , α2 , α3 , α4 )k
r2 (θ, ψ, φ1 , η1 , α1 , α4 ) = R(θ, x1 ) R(α1 , n̂14 ) x1 .

For η1 6= 0, π, we obtain2

(14)

(15)
(16)

x1 × T14
kx1 × T14 k

(17)

T14 = R(ψ, x1 ) T1

(18)

n̂14 =

T1 =


∂ 
R(u, k̂) x1
∂u

2

(19)
u=0

Similar formulas can be obtained when η1 = 0 or η1 = π although a simple solution is to chose another basis,
for instance {k̂, î, ĵ}.

6

Table 1: Desired points for the prescribed timing path generation task.
Point number
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

Point
(0.85737,
(0.82985,
(0.80241,
(0.77567,
(0.75011,
(0.72607,
(0.70381,
(0.68352,
(0.66533,
(0.64933,
(0.63559,
(0.62415,
(0.61504,
(0.60833,
(0.60396,
(0.60196,
(0.60230,
(0.60485,
(0.60959,
(0.61637,
(0.62500,
(0.63530,
(0.64700,
(0.65989,
(0.67370,
(0.68811,
(0.70293,
(0.71789,
(0.73274,
(0.74737,
(0.76167,
(0.77544,

-0.18481,
-0.20167,
-0.21996,
-0.23967,
-0.26056,
-0.28244,
-0.30515,
-0.32833,
-0.35185,
-0.37537,
-0.39867,
-0.42159,
-0.44389,
-0.46541,
-0.48596,
-0.50548,
-0.52381,
-0.54085,
-0.55656,
-0.57081,
-0.58363,
-0.59489,
-0.60456,
-0.61259,
-0.61896,
-0.62363,
-0.62652,
-0.62759,
-0.62678,
-0.62404,
-0.61933,
-0.61256,

Point number

Point

33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
52
54
55
56
57
58
59
60
61
62
63
64

(0.7887, -0.60370, 0.11578)
(0.80152, -0.59270, 0.07900)
(0.81378, -0.57959, 0.04311)
(0.82552, -0.56433, 0.00841)
(0.83678, -0.54700, -0.02478)
(0.84759, -0.52763, -0.05611)
(0.85807, -0.50641, -0.08530)
(0.86819, -0.48344, -0.11200)
(0.87804, -0.45889, -0.13596)
(0.88763, -0.43304, -0.15689)
(0.89704, -0.40611, -0.17448)
(0.90626, -0.37837, -0.18848)
(0.91537, -0.35022, -0.19867)
(0.92433, -0.32193, -0.20481)
(0.93322, -0.29396, -0.20667)
(0.94196, -0.26667, -0.20400)
(0.95052, -0.24048, -0.19667)
(0.95885, -0.21581, -0.18441)
(0.96685, -0.19315, -0.16707)
(0.97430, -0.17289, -0.14452)
(0.98096, -0.15541, -0.11656)
(0.98652, -0.14104, -0.08315)
(0.99052, -0.13007, -0.04430)
(0.99244, -0.12263, -0.00019)
(0.99174, -0.11870, 0.04874)
(0.98774, -0.11822, 0.10185)
(0.98000, -0.12085, 0.15807)
(0.96819, -0.12626, 0.21604)
(0.95226, -0.13415, 0.27426)
(0.93252, -0.14411, 0.33115)
(0.90956, -0.15600, 0.38515)
(0.88422, -0.16959, 0.43519)

0.48037)
0.52030)
0.55478)
0.58389)
0.60785)
0.62693)
0.64152)
0.65193)
0.65844)
0.66144)
0.66115)
0.65781)
0.65167)
0.64293)
0.63170)
0.61819)
0.60241)
0.58448)
0.56448)
0.54244)
0.51844)
0.49248)
0.46467)
0.43507)
0.40378)
0.37093)
0.33674)
0.30133)
0.26500)
0.22800)
0.19059)
0.15307)

r3 (φ(θ), ψ, φ1 , η1 , α1 , α2 , α3 , α4 ) = R(φ(θ), −x4 ) R(α3 , −n̂14 ) x4 ,

(20)

x4 = R(α4 , n̂14 ) x1 .

(21)

Finally, the output angle φ(θ) is given by:
φ(θ) = 2 tan

−1

−A ±

√

A2 + B 2 − C 2
C −B

!
where

(22)

A = sin α1 sin α3 sin θ
B = cos α1 sin α3 sin α4 − sin α1 sin α3 cos α4 cos θ
C = cos α1 cos α3 cos α4 + sin α1 cos α3 sin α4 cos θ − cos α2 .
It is worthwhile to mention that, since points x2 and x3 are not given as input parameters, we
arbitrarily can choose any of the two branches. Notice however, that if points x2 and x3 are given
7

as input parameters, then a careful choice of the branch is needed (see for example appendix C
of [37] for details).

3.2

Center of Mass Acceleration

Figure 2: Extensions of the Input and Output Links.

Let us consider that a link has the geometry of a geodesic joining two points, p1 and p2 on
the unit sphere, whose associated vectors subtends an arc α. Assuming that the link is made of
an homogeneous material, the center of mass vector of the link is
rcm;l =

2
sin(α/2) R (α/2, p̂12 ) p1 ,
α

(23)

where p̂12 = p1 × p2 /kp1 × p2 k.
Since the density of mass is a constant, the center of mass vector of the input link (without
extensions) can be calculated as:
rcm;l1 =
with
n̂12 =

2
sin(α1 /2) R (α1 /2, n̂12 ) x1
α1

x1 × r2 (θ, ψ, φ1 , η1 , α1 , α4 )
.
kx1 × r2 (θ, ψ, φ1 , η1 , α1 , α4 )k

(24)

(25)

The center of mass for the other links and extensions (Fig. 2) can be similarly calculated. Once
all the center of mass for the links and extensions have been calculated, the center of mass of the
complete mechanism is:
.X
X
rcm =
sk rcm;k
sk ,
(26)
k

k

where sk is the length of the k − th link (or extension) for k = e1 , e2 , e3 , e4 , l1 , etc.
Now, the acceleration of the center of mass will be
r̈cm = θ̈

∂ 2 rcm
∂rcm
+ θ̇2
.
∂θ
∂θ2

(27)

Since it is possible to control the angular velocity of the input link, the calculation of the acceleration of the center of mass is reduced to calculate ∂rcm /∂θ and ∂ 2 rcm /∂θ2 . Both first and second
8

derivatives are obtained when we write Eq. (26) in its dual form (see the rcm_dual_mod.f90 file
of the additional material). In order to clarify, let us find the dual version for the center of mass of
e θ̃, ũ) be the dual version of the active rotation matrix of dual angle θ̃ about
the input link. Let R(
the unit dual vector ũ. Then
r̃cm;l1 =

2
e 1 /2, ñ12 ) ˜ x̃1
sin(α1 /2) R(α̃
α1

(28)

where ˜ represents the dual multiplication of a dual matrix with a dual vector, α̃1 = [α1 , 0, 0], ñ12
the dual version of Eq. (25), which, in order to be constructed all the involved operations need to
be promoted to dual operations. Needless to say, such operations are coded in the downloadable
module. It is interesting to note that neither the sum nor the scalar multiplication, need to be
dualized, this is, of course due to the linearity of the derivative operator.

3.3

Theoretical calculation of the energy consumption.

From the theoretical point of view, the energy consumption of the spherical 4R mechanism can be
calculated as follows.
The external forces to the spherical 4R mechanism are,
• R1 : The reaction force of the frame on the x1 point.
• R2 : The reaction force of the frame on the x4 point.
• Fg : The weight of the mechanism.
• FM : The external force due to the motor (the agent who generates the rotation of the input
link).
In our description, the work done by R1 and R2 is zero, as the points x1 and x4 are fixed.
Since Fg is a conservative force, the work done by this force will be WFg = −∆U, where, choosing
the acceleration of gravity as g = −g k̂, the gravitational potential energy for the mechanism of
mass m is,
U (θ) = m g zcm (θ),
(29)
with zcm (θ) the z–coordinate of the center of mass as a function of the input angle θ.
From the work-energy theorem, a differential amount of work done by the force FM can be
written as
δWFM = (dU + dK) + δEothers ,
(30)
where K is the kinetic energy and the positive quantity δEothers represents all the other kinds of
energies we have not considered, for instance the heat dissipated by friction, etc.
The energy expended by the external agent (the motor) will be
δE = |(dU + dK)| + δEothers .

(31)

Notice that the absolute value is used, this is because the quantity (dU + dK) can be negative
and it could happen that δWFM to be zero3 , which in fact could be true. However, this fact does
3

From a theoretical point of view, if the quantity δEothers is known, one can obtain an angular velocity profile
–eigenmotion– for the input link in such a way that δWFM = 0.

9

not mean that the motor moves the mechanism without consuming energy. In our case, with
a conventional motor operating the mechanism, there is no difference between a negative work
and a positive work since with a conventional motor we cannot reuse the energy taken from the
mechanism. In the best scenario the energy consumption will be δEothers .
Considering U and K as functions of θ, using a Taylor expansion we have
U (θ + dθ) − U (θ) = U 0 (θ)dθ

(32)

dU = U 0 (θ)dθ

(33)

dK = K 0 (θ)dθ.

(34)

from where
and similarly

Thus in a cycle the energy consumption is given by
Z 2π
E=
|U 0 (θ) + K 0 (θ)|dθ + Eothers .

(35)

0

Applying this equation to the balanced (B) and unbalanced (N B) mechanisms, we have,
Z
EN B − EB =

2π


|U 0 (θ) + K 0 (θ)|N B − |U 0 (θ) + K 0 (θ)|B dθ + (Eothers;N B − Eothers;B ).

(36)

0

A quantification of Eothers is not a trivial task. However if both mechanisms are constructed
as similar as possible (except of course by the extensions), the frictional forces are also as reduced
as possible and if the deformation energy for both mechanisms are either similar or small, then
a plausible approximation is that Eothers;N B ≈ Eothers;B . Of course this approximation is not in
general true, but in our case it can be justified a posteriori as our theoretical results are in good
agreement with the experimental measure for EN B − EB .

4
4.1

Results
Optimum synthesis with minimum acceleration of the center of mass

The multiobjective optimization problem (simultaneous optimization of Eqs. (11) and (12)) was
solved by using the weighted sum method and the evolutionary algorithm known as Differential
Evolution [25]. Specifically, we use the DE/rand/1/bin strategy with the dither variant. For the
numerical results we have used a population of 50 individuals and a number of generations of
10 000 and a crossover probability of 0.9.
Assuming a constant angular velocity of 1 rad/s for the input link (we use a PID controller to regulate this velocity), the DE method gives as result the following values: x1 =
[1.0000, 0.0000, 0.0000], x4 = [0.54462, 0.80817, 0.22413], α1 = 0.40144 rad, α2 = 0.82035 rad,
α3 = 0.92505 rad, β = 0.23066 rad, γ = 0.47437 rad, e1 = 4.46818 rad, e2 = 0.94312 rad,
e3 = 0.24879 rad and e4 = 4.30884 rad. With these design parameters we obtain fpath = 0.00002
m and fcm = 0.01351 m/s2 . Unfortunately it is not possible to build this mechanism, as the links
(with extensions) would collide with the supports for the frame points (x1 and x4 ). A feasible
mechanism can be synthesized demanding that the z-coordinate for the points x1 and x4 to be zero,
10

Table 2: Parameters of the constructed spherical 4R mechanism.
1.00

x1
0.00

0.00

0.39

x4
0.92

0.00

α1
0.40

α2
0.88

α3
1.39

β
0.34

γ
0.40

e1
1.92

e2
1.06

e3
0.00

e4
2.79

Figure 3: Constructed spherical 4R mechanism.

and choosing horizontal supports for the frame points, thus avoiding collisions with the extensions.
With these considerations, we obtain the mechanism whose parameters are shown in Table 2. Fig.
3 shows a first prototype of the constructed mechanism. For this mechanism fpath = 0.00053 m,
while fcm = 0.014 m/s2 . It is worthwhile to note that if no extensions are considered for the
mechanism; i.e. ei = 0 for i = 1, ..., 4; the multiobjective optimization method yields fcm = 0.28
m/s2 which is twenty times bigger than the corresponding value for the mechanism with extensions. As a result, the shaking forces are minimized. From Fig. 4 we can see that the reduction of
the shaking force amplitude for the balanced mechanism, reaches values of 95% and in the worst
scenario of 78%, having a mean value of 89%. Notice that this reduction is independent (from a
theoretical point of view) of the used angular velocity for the input link (θ̇). Because if θ̇ in Eq.
(27) is constant, the ratio between ACoM for the balanced and unbalanced mechanisms, does only
depend on the ratio of their second partial derivatives with respect to θ.

4.2

Experimental validation

In order to validate our theoretical results (shown in Fig 4), either a measurement of the ACoM
or of the external forces is required. Such experimental measurements are not easy to be performed. Fortunately an experimental validation of theoretical results can be made by measuring
the difference in the energy consumption EN B − EB , of the balanced and unbalanced mechanisms.
The fact that the quantity under the integral sign in Eq. (36) depends on the velocity and the
acceleration of the center of mass allows to compare such theoretical quantity with its experimental
counterpart. To this end, we have manufactured both mechanisms with and without extensions,
considering a sphere of mean radius4 equal to 0.23 m. Specifically, the radii for the input link the
coupler link and the output links were of 0.250 m, 0.235 m and 0.210 m respectively. The radius
for the links forming the triangle (Fig. 3) for the coupler point was of 0.255 m. The links were
4

The links of the spherical 4R mechanism where constructed in concentric spheres of different radius. The
deviation in the radius for the links depends on the thickness of the links, in our case 6.35 mm.

11

Reduction of the
shaking force amplitude (%)

0.96
0.94
0.92
0.9
0.88
0.86
0.84
0.82
0.8
0.78

0

1

2

3

4

5

6 7 8
Time (s)

9

10 11 12 13

Figure 4: Reduction of the shaking force amplitude for the mechanism with minimum center of
mass acceleration. The mean value is of 89%.
1.6
1.4
Power (W)

1.2
1
0.8
0.6
0.4
0.2
0

0

5

10

15

20 25 30
Time (s)

35

40

45

50

Figure 5: Experimental power consumption comparison between the ACoM balanced (solid line)
and unbalanced (dashed line) mechanisms.

constructed with aluminum slab of 6.35 mm thickness and 2.54 cm width. For both mechanisms a
dc motor is used as actuator to the input link, which is controlled by a PID law with a reference of
1 rad/s. The energy required to operate the mechanism is considered equal to the electrical energy
required to actuate the motor, which was obtained from direct measurement of voltage and current
at its terminals. In Fig 5, we show the power consumption for both mechanisms. Fig 6, presents
the time behavior of the angular velocity of the input link of both mechanisms (which is obtained
as the output of the filter s/(1 + 0.01 s), whose input comes from a digital encoder attached to
the joint of the input link), we can observe that the regulation objective is not absolutely reached,
however the oscillations around the reference are smaller for the case of the balanced mechanism.
After an integration of the experimental data showed in Fig. 5, we have obtained a deviation
of only 5% with respect to a theoretical calculation of EN B − EB using Eq. (36) and neglecting
Eothers;N B − Eothers;B .

12

1.03

1.02

1.02

Angular velocity HradsL

Angular velocity HradsL

1.03

1.01
1.00
0.99
0.98
0.97

10

20

Time HsL

30

1.01
1.00
0.99
0.98
0.97

40

(a)

10

20

Time HsL

30

40

(b)

Figure 6: Controlled velocity for the unbalanced mechanism (a) and the ACoM balanced mechanism (b).

5

Discussion

Before concluding, some comments on the used derivative method and the energy consumption
are worthwhile. Below we present a brief discussion.

Numerical derivatives: Dual numbers vs Finite Differences
Using finite differences (FD), the second derivative of a function can be calculated as,
f 00 (x) =

f (x + h) − 2 f (x) + f (x − h)
+ O(h2 ).
h2

(37)

In order to obtain the second derivative using this expression, three evaluations of the function
f are required. Since in the dual number approach we only need one evaluation, it is reasonable
to think that this method will be more efficient. However, we cannot say that the dual number
approach will be three times faster than the finite difference method, as the involved evaluations
do not correspond to the same functions. In our case, for 5000 points equally spaced in (0, 2π),
the matrix [ṙcm , r̈cm ] was calculated 2.5 times faster by using the dual numbers approach instead
of the FD method.
Now, from Eq. (37) a truncation error of order h2 is present. One could think that this problem
can be solved by taking a small value for h, but then subtractive cancellation errors will appear.
For example taking h = 1 × 10−6 and using double precision for the variables, the error in the
calculation of the ACoM is not of the order h2 but of the order of 10−2 . This situation can be
alleviated by writing a more elaborated formula for f 00 but there is a cost to pay, we will need more
evaluations of the function f . These problems are not present when the dual numbers approach is
used to calculate the derivatives. Moreover, this fact is not the only advantage of the dual number
approach. For instance, one could be interested in obtaining the derivative of a function that is the
solution of some equation which can only be known numerically. In this case, we can implement
the numerical solution method in the context of the dual number approach and then both, the
solution and its derivative will be directly obtained.

13

Energy consumption
Even when minimization of ACoM does not imply a constant gravitational potential energy –
excepting if ACoM is zero and for a mechanism operating in cycles5 –, we have found (integrating
the absolute value of Eq. (33)) that the variations in the potential energy are reduced about 92%
for the balanced mechanism.
It is worthwhile to mention that not only the shacking forces were reduced by balancing the
mechanism but also the energy consumption (obtained by integrating the experimental data of
Fig. 5) was reduced by a 21 %. Of course this is not a rule since we are not optimizing Eq. (35).
Moreover, unlike the percentile reduction of the total external force, the percentile reduction of
the energy consumption depends on the angular velocity of the input link.

6

Conclusions

By implementing the center of mass vector in its dual version, a simple yet effective method for
balancing a spherical 4R mechanism using the counterweights method has been presented. The
minimization of the magnitude of the center of mass acceleration vector has been conducted with
the Differential Evolution algorithm. The synthesis process for the prescribed path generation task
turn out to be highly efficient giving excellent results. The reduction of the force shaking amplitude,
has reached a mean value of 89% while the variations in the gravitational potential energy have been
reduced about 92%. The synthesized mechanisms, the balanced and its unbalanced counterpart,
were constructed in order to experimentally measure the power consumption and hence the energy
needed for both mechanisms to perform the requested path. By doing this, an experimental
validation of theoretical results has been conducted. As a side effect, it is interesting to note
that the energy used to operate the balanced mechanism was reduced by 21% with respect to its
unbalanced counterpart. Nevertheless this cannot be considered as a rule as a minimization of
ACoM does not necessarily implies a minimization of Eq (35).

Acknowledgments
Authors thank financial support from the mexican PROMEP and the National System for Researchers. Authors also acknowledge the anonymous reviewers whose comments and suggestions
greatly improved this work.

5

Let us assume that ACoM is zero, this means that the velocity of the center of mass is a constant. Then, in
a cycle the mechanism must return to its initial point, therefore the only allowed constant for the velocity of the
center of mass is zero. Thus GPE is constant. Notice that the converse is not true.

14

References
[1] Y. Wu, C. Gosselin, Synthesis of reactionless spatial 3-dof and 6-dof mechanisms without
separate counter-rotations, The International Journal of Robotics Research 23 (6) (2004)
625–642.
[2] Q. Jiang, C. Gosselin, Dynamic optimization of reactionless four-bar linkages, Journal of
Dynamic Systems, Measurement, and Control 132 (2010) 041006.
[3] V. van der Wijk, J. L. Herder, B. Demeulenaere, Comparison of various dynamic balancing principles regarding additional mass and additional inertia, Journal of Mechanisms and
Robotics 1 (4) (2009) 041006–1–041006–9.
[4] V. van der Wijk, S. Krut, F. Pierrot, J. L. Herder, Design and experimental evaluation of a
dynamically balanced redundant planar 4-rrr parallel manipulator, The International Journal
of Robotics Research 32 (6) (2013) 744–759.
[5] C. Bagci, Complete shaking force and shaking moment balancing of link mechanisms using
balancing idler loops, Journal of Mechanical Design, ASME 104 (1982) 482–493.
[6] G. Feng, Complete shaking force and shaking moment balancing of 26 types of four-, five- and
six-bar linkages, Mechanisms and Machine Theory 25 (1990) 183–192.
[7] V. Arakelian, M. Smith, Complete shaking force and shaking moment balancing of linkages,
Mechanisms and Machine Theory 34 (6) (1999) 1141–1153.
[8] B. Demeulenaere, E. Aertbeliën, M. Verschuure, J. Swevers, J. D. Schutter, Ultimate limits
for counterweight balancing of crank-rocker four-bar linkages, Journal of Mechanical Design,
ASME 128 (6) (2006) 1272 (13 pages).
[9] J. Wang, C. M. Gosselin, Static balancing of spatial three–degree–of–freedom parallel mechanisms, Mechanism and Machine Theory 34 (1999) 437–452.
[10] A. Agrawal, S. K. Agrawal, Design of gravity balancing leg orthosis using non–zero free length
springs, Mechanism and Machine Theory 40 (2005) 693–709.
[11] S. R. Deepak, G. Ananthasuresh, Static balancing of a four-bar linkage and its cognates,
Mechanism and Machine Theory 48 (2012) 62–80.
[12] G. G. Lowen, F. R. Tepper, R. S. Berkof, Balancing of linkages– an update, Mechanisms and
Machine Theory 18 (1983) 213–230.
[13] V. Arakelian, M. R. Smith, Shaking force and shaking moment balancing of mechanisms:
A historical review with new examples, Journal of Mechanical Design, ASME 127 (2005)
334–339.
[14] A. Russo, R. Sinatra, F. Xi, Static balancing of parallel robots, Mechanism and Machine
Theory 40 (2005) 191–202.
[15] V. van der Wijk, J. L. Herder, Synthesis method for linkages with center of mass at invariant
link point — pantograph based mechanisms, Mechanism and Machine Theory 48 (2012) 15–28.
15

[16] S. Briot, V. Arakelian, J.-P. L. Baron, Shaking force minimization of high-speed robots via
centre of mass acceleration control, Mechanism and Machine Theory 57 (2012) 1–12.
[17] N. M. Qi, E. Pennestrı́, Optimum balancing of four-bar linkages, Mechanism and Machine
Theory 26 (1991) 337–348.
[18] L. Zhe, B. Shixian, Optimum balancing of linkages with clearances, Mechanism and Machine
Theory 27 (1992) 535–541.
[19] H. Chaudhary, S. K. Saha, Balancing of four-bar linkages using maximum recursive dynamic
algorithm, Mechanism and Machine Theory 42 (2007) 216–232.
[20] H. Chaudhary, S. K. Saha, An optimization technique for the balancing of spatial mechanisms,
Mechanism and Machine Theory 43 (2008) 506–522.
[21] R. Saravanana, S. Ramabalanb, P. D. Babub, Optimum static balancing of an industrial robot
mechanism, Engineering Applications of Artificial Intelligence 21 (2008) 824–834.
[22] W. Clifford, Preliminary sketch of biquaternions, Proc. London Mathematical Society 1 (s1-4)
(1873) 381–395.
[23] V. Brodsky, M. Shoham, Dual numbers representation of rigid body dinamics, Mechanism
and Machine Theory 34 (1999) 975–991.
[24] E. Pennestrı̀, P. Valentini, Linear dual algebra algorithms and their application to kinematics,
Multibody Dynamics Computational Methods and Applications 12 (2008) 207–229.
[25] K. V. Price, R. M. Storn, J. A. Lampinen, Differential Evolution: A Practical Approach to
Global Optimization, Springer, Germany, 2005.
[26] Y.-L. Gu, J. Y. S. Luh, Dual-number transformation and its applications to robotics, IEEE
Journal of Robotics and Automation 3 (6) (1987) 615–623.
[27] H. H. Cheng, Programming with dual numbers and its applications in mechanisms design,
Engineering with Computers 10 (4) (1994) 212–229.
[28] F. Peñuñuri, R. Peón-Escalante, C. Villanueva, C. A. Cruz-Villar, A dual number approach for numerical calculation of derivatives and its use in the spherical 4R mechanism,
arXiv:1301.1409 [].
[29] J. Cabrera, F. Nadal, J. P. M. noz, A. Simon, Multiobjetive constrained optima synthesis
of planar mechanisms using a new evolutionary algorithm, Mechanism and Machine Theory
42 (7) (2007) 791–806.
[30] R. R. Bulatović, S. R. Dordević, On the optimum synthesis of a four-bar linkage using differential evolution and method of variable controlled deviations, Mechanism and Machine Theory
44 (1) (2009) 235–246.
[31] F. Peñuñuri, R. Peón-Escalante, C. Villanueva, D. Pech-Oy, Synthesis of mechanisms for
single and hybrid tasks using differential evolution, Mechanism and Machine Theory 46 (2011)
1335–1349.
16

[32] J. Cabrera, A. Ortiz, F. Nadal, J. Castillo, An evolutionary algorithm for path synthesis of
mechanisms, Mechanism and Machine Theory 46 (2011) 127–141.
[33] G. Mullineux, Atlas of spherical four-bar mechanisms, Mechanism and Machine Theory 46
(2011) 1811–1823.
[34] J. Chu, J. Sun, Numerical atlas method for path generation of spherical four-bar mechanism,
Mechanism and Machine Theory 45 (2010) 867–879.
[35] J. M. McCarthy., Geometric Design of Linkages, Springer, New York, 2000.
[36] F. Peñuñuri, R. Peón-Escalante, C. Villanueva, C. A. Cruz-Villar, Synthesis of spherical 4R
mechanism for path generation using differential evolution, Mechanisms and Machine Theory
57 (2012) 62–70.
[37] J. J. Cervantes-Sánchez, H. I. Medellı́n-Castillo, J. M. Rico-Martı́nez, E. J. González-Galván,
Some improvements on the exact kinematic synthesis of spherical 4R function generators,
Mechanism and Machine Theory 44 (2009) 103–121.

17

