Optimum Synthesis of Mechanism for single- and hybrid-tasks using
Differential Evolution
F. Peñuñuri,∗ R. Peón-Escalante,† C. Villanueva,‡ and D. Pech-Oy§

arXiv:1102.2017v2 [] 18 Jun 2011

Facultad de Ingenierı́a, Universidad Autónoma de Yucatán,
A.P. 150, Cordemex, Mérida, Yucatán, México.
The optimal dimensional synthesis for planar mechanisms using differential evolution (DE) is
demonstrated. Four examples are included: in the first case, the synthesis of a mechanism for
hybrid-tasks, considering path generation, function generation, and motion generation, is carried
out. The second and third cases pertain to path generation, with and without prescribed timing.
Finally, the synthesis of an Ackerman mechanism is reported. Order defect problem is solved by
manipulating individuals instead of penalizing or discretizing the search space for the parameters.
A technique that consists in applying a transformation in order to satisfy the Grashof and crank
conditions to generate an initial elitist population is introduced. As a result, the evolutionary
algorithm increases its efficiency.

I.

INTRODUCTION

Dimensional synthesis of mechanisms comprises the problems of path, function and motion generation. There are three types of methods for this purpose: graphical, analytical, and those involving
optimization [1].
Graphical methods offer a quick solution by sacrificing accuracy, and are rarely used since computers
can do the same work faster and better.
Analytical methods are based on algebraic expressions [1, 2], displacement matrix [3], complex
numbers [4], or continuation methods [5] resulting in mechanisms whose error will be zero at the
precision points.
The problem of motion generation, in the case of a planar four-bar mechanism, can be designed
based on the Burmester curve. This is one of the first proposed analytical methods for the dimensional
synthesis of mechanisms. In [6] an algorithm for the robust computation of the solution of the fiveposed Burmester problem is introduced. In [7] a Matlab-based graphical user interface to the algorithm
of [6] is done.
Also, the general equation of the coupler curve of a four-bar linkage has attracted the attention of
researchers. For a given set of points on the coupler curve, Blechschmidt and Uicker [8] have used the
equation of the coupler curve to synthesize a four-bar linkage by determining the coefficients of the
curve.
The main disadvantage of the analytical methods lies in the maximum number of points of accuracy
that can be set. The mechanisms are restricted to move exactly in a number of points equal to the
number of independent parameters that define them [9, 10]. Even though the mechanisms obtained
can reach the precision points, they may have other problems, known as design defects, that are not
taken into account during the synthesis process, thereby preventing the mechanisms from fulfilling the
task for which they were designed [9].
Optimization methods are based on numerical methods and allow a large number of design points
tolerating a loss of accuracy. These are formulated in terms of nonlinear programming problems.
The optimal solution is found by optimizing an objective function within an iterative procedure. The
objective function can be defined as a difference between the generated and the specified movement,
known as the structural error [3]. In general, it can be defined as the design error, i.e., the error
∗
†
‡
§

francisco.pa@uady.mx
rpeon@uady.mx
cesar.villanueva@uady.mx
davidj@hotmail.com

2
that arises when we are trying to satisfy a design equation [11] (which could be the Freudenstein
equation). An interesting definition for the objective function is presented in [12] where it is defined as
a kind of entropy that is maximized. The use of optimization methods is inevitable when the number
of positions to be covered during the duty cycle exceeds a certain number (in the case of motion
generation synthesis the classical analytical approach is limited to five specified points for a four-bar
mechanism).
The interest in optimum synthesis of mechanisms is not new. There have been a large number
of studies on this topic using a variety of methods. For example, some local search methods have
been described in references [13–20]. The main disadvantages of these methods are that the objective
function must be differentiable. Also, they are very sensitive to the initial search point.
Within the global search methods some of the techniques that have been used are Simulated Annealing (SA) [21], Neural Network [23, 24], Genetic Algorithm (GA) [25–30], Particle Swarm Optimization
Technique (PSO) [30], and Differential Evolution (DE) [30–34]. There are works that use a combination of two optimization methods such as SA-Powell’s Method [35], GA-FL [36], Tabu-Gradient [37],
Ant Colony Optimization-Gradient (AG) [38], and GA-DE [39].
The use of evolutionary algorithms has been of significant interest in recent years. For instance,
Ullah and Kota solved the path generation problem by presenting an objective function based on
Fourier descriptors that evaluates only the shape differences between two curves [35]. This function
is first minimized using a simulated annealing followed by Powell’s method. The size, orientation
and position of the desired curve are addressed at a later stage by determining analogous points on
the desired and candidate curves. Similarly, Vasiliu and Yannou [23] synthesized the dimensions of a
planar mechanism whose purpose is to generate a trajectory shape by using a neural network.
Laribi et al. presented the combined GA-FL method to solve the problem of path generation in
mechanism synthesis [36]. The FL-controller monitors the variation of the design variables during the
first run of the GA and modifies the initial bounding intervals to restart a second run of the GA.
Smaili and Diab (2007) apply AG to the mechanism synthesis problem for single- and hybrid-tasks
[38]. Shiakolas introduced a technique called the Geometric Centroid of Precision Points for defining
initial bounds for the design variables combining with DE [31].
Acharyya and Mandal carry out the path synthesis of a four-bar linkage using three different methods
[30]. They found that the DE with /rand/1/exp method performs better than the two others; one
being a binary-coded genetic algorithm (BGA) with multipoint crossover, and the other a PSO with
the constriction factor approach.
In [39] they used a GA-DE hybrid algorithm to make a path synthesis of a four-bar linkage. A
real-valued genetic algorithm, where the crossover operation of GA is replaced by differential vector
perturbation, is employed.
The DE method is a simple yet powerful algorithm for global optimization [40]. It is not difficult
to modify the main operators and try for improvements of the method. In the present work, we use
DE to find optimum solutions for the dimensional synthesis problem of four mechanisms. The first
three correspond to planar four-bar and the last one to a six-bar mechanism. The paper is organized
as follows: In section II we present the classical DE method, which is used throughout this work;
section III presents notation and conventions. In Section IV we employ the idea of hybrid-task for
the synthesis of mechanisms as was introduced in [38]. The problem of this section presents us with
the difficulty of mixing angles with lengths. This difficulty is addressed by introducing a factor that,
on the one hand, defines consistently the objective function and on the other hand, allows for proper
weighing of the involved errors. This is important in order to fulfill the task of function and motion
generation. Moreover, this problem is used to show an easy and effective way to handle the order
defect problem. The proposed method avoids entirely both individual penalization and space search
discretization.
Section V deals with the prescribed timing path generation for 18 points and 10 design variables.
We introduce a transformation which constructs an elitist population, in the sense of satisfying the
Grashof and crank conditions, avoiding a probabilistic or penalization approach. This problem has
been presented by other authors [26, 27, 39]. To avoid some controversies related to the values of the
objective function that each of them report, we have written a Fortran 90 program that evaluates
the objective function.

3
Start

Specify the DE parameters

Initialize the population

Generate a trial population:
Mutation, Crossover, Selection

No

Termination criteria
satisfied?
Yes
Stop

FIG. 1. Flowchart for the DE algorithm.

The ideas introduced in Sections IV and V allow us to solve in a single manner the path generation
problem without prescribed timing for 18 target points and 27 design variables, which is the problem
described in Section VI. In Section VII an Ackerman mechanism is optimized. Finally, we present our
conclusions in Section VIII.

II.

CLASSICAL DE

Below, the original version of the method is outlined [41].
1. The population:
Px,g = (xi,g ), i = 1, ...m; g = 0, ...gmax
xi;g = (xji;g ), j = 1, ...D;

(1)

where D, m and gmax represent the dimensionality of x, the number of individuals and the
number of generations respectively. In [42] it is mentioned that a good choice for m is 10D.
However, to balance the speed and reliability in [43] values from 2D to 40D are suggested.
2. Initialization of population:
xji;0 = randj (0, 1) · (bjU − bjL ) + bjL .
Vectors bU and bL are the parameter limits and randj (0, 1) is a random number in [0, 1) generated
for each parameter.

4
3. Mutation:
vi;g = xr0 ;g + F · (xr1 ;g − xr2 ;g ).

(2)

The main difference between DE and other evolutionary algorithms like GA comes from 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 is a scale factor greater than zero. Even
though upper limits for F do not exist, values greater than 1 are rarely chosen in the literature
[40–42, 44].
4. Crossover:
A dual recombination of vectors is used to generate the trial vector:
(
j
if(randj (0, 1) 6 Cr or j = jrand )
vi;g
j
ui;g = ui;g =
j
xi;g otherwise.

(3)

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

xi;g+1 =

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

(4)

The method just described is known as DE/rand/1/bin. There are variants of it. For example, when
F is chosen to be a random number, the variant is called dither. In this work we will use the exposed
method with the dither variant where F ∈ [0; 1). Fig. 1 shows the flowchart for the DE algorithm.
III.

MECHANISM SYNTHESIS PROBLEM: NOTATION AND CONVENTIONS

The simplicity of a 4-bar mechanism, (easy to manufacture and highly reliable) makes it a very
important mechanism with a large number of industrial applications. Its use ranges from simple
devices such as windshield-wiping mechanisms and door-closing mechanisms to more complicated ones
such as rock crushers, sewing machines, round balers, and suspension systems of automobiles [30].
In this section the notation and conventions used throughout this work are established. The only
exception is in Section VII where we will deal with a 6-bar mechanism.
A four-bar linkage shown in Fig. 2 consists of four rigid links and four revolute joints. The set
of variables that describes the mechanism (the design variable vector) will be put into the vector X
whose components will be enclosed within braces. Usually, in the synthesis of a mechanism there are
two sets of points (or coordinates), desired and generated points, allocated in the vectors rd and rgen ,
respectively. A vector error E = rd − rgen is proposed and the objective function is defined as the
square of its Euclidean norm.
f ob = |E|2 .

(5)

If quantities are not dimensionally homogeneous, constants with appropriate units must be introduced
so that equations have compatible units. In this work, there are quantities with different units, and
some constants are chosen so that f ob is dimensionless. For example, in the problem of motion
generation, we have to fit angles and coordinates, so the quadratic error will be
X

E2 =
c(xi;d − xi;gen )2 + c(yi;d − yi;gen )2 + (θi;d − θi;gen )2 ,
(6)
i

where xi;d[gen] , yi;d[gen] and θi;d[gen] are the coordinates and angles of the desired [generated] point i.

5

FIG. 2. Four-bar linkage notation.

The constant c has a numerical value equal to 1 and is introduced for consistency with units. The
objective function is given by
X

f ob =
fc2 (xi;d − xi;gen )2 + fc2 (yi;d − yi;gen )2 + (θi;d − θi;gen )2 .
(7)
i

The fc constant is introduced for consistency with units but is not necessarily 1. Such a constant can
be defined by the user as a weight factor. In this work it is chosen as the inverse of the longest distance
between the coordinates. As a matter of illustration, for the points P = {(1, 1), (2, 3), (−5, −1)} we
construct the set Uxy = {1, 2, 3, −5, −1} (i.e., the union of the coordinates) and take fc = 1/dr with
dr = max(Uxy ) − min(Uxy ). In this case min(Uxy ) = −5, max(Uxy ) = 3 thus fc = 1/8. Notice that
this definition is motivated by the curvature concept. For example, in the case of a circle with radius
r, we have s/r = θ or ks = θ where k is the curvature, s the arc length subtended by the angle θ.
In general, f ob 6= 0 and its minimization process is what generates values for the parameters of a
possible mechanism. In the analysis of mechanisms, two conditions are important. They are known as
the crank and Grashof conditions (CG):
min(r1 , r2 , r3 , r4 ) = crank,
2 min(r1 , r2 , r3 , r4 ) + 2 max(r1 , r2 , r3 , r4 ) < r1 + r2 + r3 + r4 .

(8)
(9)

In our case r2 is the crank, see Fig 2. Whenever we refer to a transformation acting on a vector,
|xi is used instead of x. Usually such transformations are carried out by subroutines or functions in
Fortran 90 and by functions in C++. For linear transformations, the matrix representation can be
used. In this work, all the algorithms for the synthesis of mechanisms were implemented in Fortran
90. The compiler used was ifort and the calculations were made in an intel Core 2 Duo processor with
velocity of 2.53 GHz, 4 GB of memory and a bus velocity of 1.07 GHz.
IV.

HYBRID TASK SYNTHESIS

In this section we analyze the problem presented by McGarva [45]. We address the problem from
the viewpoint of hybrid tasks as proposed by Smaili and Diab in [38]. The problem has three tasks:
function generation, motion generation and path generation. Table I (as presented in [38]) summarizes
the variables used in this study.
The design variable vector is defined as
X = {x0 , y0 , r1 , r2 , r3 , r4 , rcx , rcy , ψ1 , ψ2 , ψ3 , ψ4 , ψ5 , ψ6 , ψ7 , γ}.

(10)

6
Desired point, i 1

2

3

4

5

6

7

8

9

10

4.67
2.6
122
*
138

4.38
2.2
137
*
143

4.04
1.67
151
*
147

Function points

xid
yid
ψi + γ
θi + γ
φi + γ

7.03
5.99
21
*
108
11

6.95
5.45
36
*
110
12

6.77
5.03
50
*
113
13

6.4
4.6
65
*
117
14

5.91
4.03
79
*
121
15

5.43
3.56
93
*
126
16

4.93
2.94
108
*
132
17

Motion function

xid
yid
ψi + γ
θi + γ
φi + γ

3.76
1.22
N
−13
*
18

3.76
1.97
N
−7
*
19

3.76
2.78
N
−2
*
20

3.76
3.56
N
2
*
21

3.76
4.34
N
7
*
22

3.76
4.91
N
11
*
23

3.76
5.47
N
14
*
24 25

Path point

xid
yid
ψi + γ
θi + γ
φi + γ

3.8
5.98
266
*
*

4.07
6.4
281
*
*

4.53
6.75
295
*
*

5.07
6.85
309
*
*

5.05
6.84
324
*
*

5.89
6.83
338
*
*

6.41
6.8
353
*
*

6.92
6.58
367
*
*

TABLE I. Hybrid-tasks problem; N: Generated crank angle values, *: Non-prescribed values.

In addition to the CG restrictions, we have the following constraints for motion generation:
ψmin < ψ j < ψmax ; j = {1, 2, ..., 7}
ψ k < ψ k+1 ; k = {1, 2, .., 6}.

(11)

In this case the objective function consists of three parts:
f ob = f obf unc + fg
obmot + fg
obpath

(12)

where, in the usual approach of the least square method, the partial objective functions are defined as:
X
f obf unc =
(θi;d − θi;gen )2f unc
(13)
i

fg
obmot

X

=
fc2 (xi;d − xi;gen )2mot + fc2 (yi;d − yi;gen )2mot + (θi;d − θi;gen )2mot

(14)

i

fg
obpath =

X

fc2 (xi;d − xi;gen )2path + fc2 (yi;d − yi;gen )2path

(15)

i

The evaluation of the weight factor fc is explained in Sec. III. In this case we have max(Uxy ) = 7.03,
min(Uxy ) = 1.22 thus fc2 = 0.02 with units of length−2 .
The following values were tested for Cr: 0.05, 0.1, ..., 0.9. It turns out that 0.3 gives the best results.
The number of individuals and generations were m = 250, gmax = 15 000, respectively. The evaluation
of f ob resulted in a value of 6.99 × 10−3 with the design variables shown in Table II.
x0
y0
r1
r2
r3
r4
rcx
rcy
-8.0339 1.07673 13.2425 1.96639 7.71759 7.57298 13.4593 3.13037
ψ1
ψ2
ψ3
ψ4
ψ5
ψ6
ψ7
γ
3.5639 3.83348 4.05641 4.22857 4.48498 4.71726 4.92507 5.83047
TABLE II. Parameter values of an optimal mechanism. Hybrid-tasks synthesis.

For the searching space we have used the limits:
xmin = {−15, −15, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3, 0}
xmax = {15, 15, 15, 15, 15, 15, 15, 15, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 5.03, 2π}.

7
A.

Constraints Management

In this case the CG conditions do not play any active role. We can just verify that they are met
after the minimum of f ob is obtained. Concerning the requirement of the constraints of Eq. (11),
previous methods are based on the discretization of the search space for ψ angles [38],
h
i
j
j
ψ j ∈ ψmin
, ψmax
.
(16)
The best fitting angle is selected from this range. In our case this discretization is not applied, and
individuals xi,g are chosen so as to comply with the restriction of Eq. (11). To this end, a random
vector of angles within desired limits is generated, and its coordinates written in ascending order. This
idea has been used in [30, 39]. Thus the method here is not exactly a classic DE because the evolution
of individuals is manipulated. However, it is clear that the results will be the same. The only thing
it does is to accelerate the evolutionary process. Symbolically, if |ψi represents a vector of random
d represents a transformation that puts them in ascending order, then
numbers and sort
h
ij
d
ψ j = sort|ψi
.
(17)
There are several ways to implement Eq. (17). In particular, it can be done in the crossover part.
(
j
ṽi;g
if(randj (0, 1) 6 Cr or j = jrand )
j
ui;g = ui;g =
(18)
j
x̃i;g otherwise.
where
h
ij
d
r̃j = sort|ri
.

(19)

d will act only on those components that we choose to order.
The transformation sort

FIG. 3. Optimal mechanism and the corresponding coupler curve. Hybrid-tasks synthesis.

For the ordering of the ψ angles, we have used the heap sort method [46–48] as it is efficient enough
and easy to implement.
Penalizing angles ψ j is not very efficient because the probability of having a set of size n randomly
ordered is low if n is large. For example, the probability to throw in seven random numbers between 0
and 1 (or any other continuum interval) in an ordered way is 1/7!, which is about 2 × 10−4 . We thus
end up with a method without individuals to evolve unless the number of individuals in the initial

8
population were extremely large, which would lead to a grossly inefficient method. Proceeding as [38]
is a brilliant possibility and the results so obtained are very good. However, discretizing the searching
space could prevent us from locating the minimum. Fig. 3 shows the mechanism obtained. In B a
program in M athematica R that makes an animation of the mechanism is shown step by step.
V.

A CLASSICAL COMPARISON: PATH GENERATION FOR 18 TARGET POINTS
AND 10 DESIGN VARIABLES

Recently, in [39] a hybrid method (GA-DE) was proposed that can synthesize a four-bar mechanism
and the problem of prescribed timing path generation for 18 points, (previously introduced by [26] and
[27]) is addressed. We will optimize this problem by using a DE algorithm. The objective function
value is lower than the reported values of previous references. It is worth mentioning that the values
for the links of the mechanism generated are of the same order of magnitude of the generation path
dimensions.

FIG. 4. Optimal mechanism and the corresponding coupler curve. Prescribed timing path generation.

A.

The problem

The target points are:
xd = {0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.02, 0, 0, 0.03, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.6}
yd = {1.1, 1.1, 1.1, 1, 0.9, 0.75, 0.6, 0.5, 0.4, 0.3, 0.25, 0.2, 0.3, 0.4, 0.5, 0.7, 0.9, 1}

(20)
(21)

The design variable vector is
X = {x0 , y0 , r1 , r2 , r3 , r4 , rcx , rcy , γ, ψ 0 }

(22)

and the precribed timing is defined by
ψk = ψ0 +

π
(k − 1); k = {1, 2, ..., 18}.
9

(23)

9
Figure 2 shows the design variables.
This problem has been discussed in Refs. [26, 27, 39]. They do not show the explicit form of f ob,
and there is a controversy concerning the numerical values of the objective function. We show in Table
III the values that according to [39] the other two references should have obtained.
Wen-Yi Lin [39]
−2

f ob = 1.08613 × 10

Kunjur and Krishnamurty [26] Cabrera etal [27]
f ob = 1.09034 × 10−2

f ob = 3.48391 × 10−2

TABLE III. Values for the objective function reported by [39].

Here we get the following values of f ob for the design variable vectors that they report: [39],
f ob = 1.0306 × 10−2 ; [26], f ob = 1.0214 × 10−2 ; [27], f ob = 3.3748 × 10−2 . They are slightly
different from the values of [39], perhaps because of rounding errors. With the purpose of avoiding
any misunderstanding, in A we show a Fortran 90 program that evaluates f ob.
Table IV shows the values for the design variable vector for which the objective function is 9.088 ×
10−3 . The values 0.1, 0.2, . . . , 0.9 were tested for Cr. It turns out that 0.3 gives the best results.
Figure 4 shows the optimum mechanism and its path.
x0
y0
r1
r2
r3
r4
rcx
rcy
γ
ψ0
0.27892 0.11673 1.08913 0.42259 0.96444 0.58781 0.39137 0.42950 0.32195 0.86323

TABLE IV. Parameter values of an optimal mechanism with f ob = 9.088 × 10−3 .

In order to obtain the last result for f ob, we proceed in two steps. First, we choose parameter values
inside the interval [−1.5, 1.5] for x0 and y0 . For the remaining parameters we choose values in [0, 1.5],
and we evaluate f ob over and over until we find a design variable vector for which f ob 6 5 × 10−2 .
Second, from the obtained parameters, the searching space is reduced to
vxmin = {0.2, 0.1, 0.8, 0.3, 0.7, 0.4, 0.2, 0.3, 0.1, 0.7}
vxmax = {0.3, 0.3, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1}.

(24)
(25)

The value of f ob = 9.088 × 10−3 was obtained for 200 individuals and 11 817 generations. We could
reach smaller values of f ob if the generation number and/or individual number were increased, but
improvements are not considerable. For example, for 200 individuals and 30 000 generations we obtain
f ob = 9.06 × 10−3 . Moreover, by making a third refinement of the searching space, we obtain f ob =
9.03 × 10−3 for the design variables shown in Table V.
x0
y0
r1
r2
r3
r4
rcx
rcy
γ
ψ0
0.26439 0.16956 1.04028 0.42446 0.89397 0.60308 0.36129 0.38864 0.26873 0.90493

TABLE V. Parameter values of an optimal mechanism with f ob = 9.03 × 10−3 .

B.

On steps 1 and 2

In this work we subdivide the optimization task in two steps. In the first step we use an elitist
population in the sense of choosing only those individuals that satisfy the CG condition. To this end
we construct a transformation that takes an individual that does not satisfy the CG condition and
turns it into one that does. Then, in the second stage (with the result for the possible mechanism
obtained in this first stage) we refine the searching space, remove the CG condition and re-run the
optimization program. The process terminates when some criteria have been met and the individual
satisfies the CG condition.

10
C.

On the construction of the elitist population

In order to construct individuals satisfying the CG conditions we could proceed in a random way,
but this would be inefficient. In this work we proceed as follows:
Assuming the links of the four-bar mechanism belong to the interval (0,1), four random numbers
are generated (the links) and they are sorted in ascending order. At this point, proceeding randomly
would not be a bad choice since the probability of satisfying the CG condition for the sorted list is
0.5. However a better choice will be to construct a transformation T̂ that makes the CG condition
fulfill the ascending order list, |xi. There are many possible forms for the transformation T̂ . We define
T̂ = F̂ R̂, where R̂ is defined as the transformation that inverts the components of a vector and F̂ a
reflection plus a translation. Symbolically,
R̂|x1 , x2 , x3 , x4 i = |x4 , x3 , x2 , x1 i,
F̂ |xi = −|xi + |1i,
|1i = |1, 1, 1, 1i .

(26)
(27)
(28)

If the upper limit for the links is |Li, we replace |1i by |Li.
Notice that R̂ is a linear transformation, whereas F̂ is not. Once the vector that satisfies the Grashof
condition has been constructed, the crank is taken as the lesser of the elements; thus the conditions
CG will be satisfied.
For example, suppose that we have the four numbers xr = {0.38, 0.98, 0.25, 0.19} which do not
satisfy the CG condition. After sorting them we have |xi = |0.19, 0.25, 0.38, 0.98i and
R̂|xi = |0.98, 0.38, 0.25, 0.19i,
F̂ R̂|xi = | − 0.98, −0.38, −0.25, −0.19i + |1, 1, 1, 1i,
F̂ R̂|xi = |0.02, 0.62, 0.75, 0.81i.
By choosing the crank as 0.02, we can see that xg = {0.02, 0.62, 0.75, 0.81} satisfies the CG conditions
since min(xg) = 0.02, max(xg) = 0.81 and 0.02 + 0.81 < 0.62 + 0.75.
In general, suppose we have four positive numbers less or equal than 1 that are sorted in ascending
order, but that do not satisfy the CG conditions. Let |xi = |x1 , x2 , x3 , x4 i be the vector containing
such numbers. We have x1 + x4 > x2 + x3 since the numbers are sorted and do not comply Eq. (9).
Clearly −x1 − x4 < −x2 − x3 and (1 − x1 ) + (1 − x4 ) < (1 − x2 ) + (1 − x3 ). For these four constructed
numbers, the minimum is (1 − x4 ) and the maximum is (1 − x1 ) so the CG conditions are satisfied if
we chose the crank as (1 − x4 ).
VI.

PATH GENERATION WITHOUT PRESCRIBED TIMING FOR 18 TARGET POINTS
AND 27 DESIGN VARIABLES

It is interesting to synthesize the above mechanism without the prescribed timing Eq. (23). Finding
a minimum for the objective function is now more difficult. We have 27 design variables and the order
defect problem appears hard to solve.
Let
X ={x0 , y0 , r1 , r2 , r3 , r4 , rcx , rcy , γ, ψ 0 , ψ 1 , ψ 2 , ψ 3 , ψ 4 , ψ 5 , ψ 6 , ψ 7 , ψ 8 , ψ 9 , ψ 10 , ψ 11 , ψ 12 , ψ 13 , ψ 14 ,
ψ 15 , ψ 16 , ψ 17 }

(29)

be the design variable vector. Besides the CG conditions we also have the requirement
ψ k < ψ k+1 ; k = {0, 1, .., 16}
which prevents the order defect.

(30)

11
The use of penalization for the restriction Eq. (30) is not effective. Practically all the individuals
would be penalized as the probability of finding one that would not is 1/18! – a very small probability.
If we discretize the searching space for angles then there is no guarantee that the minimum will lie in
the generated intervals. However, if we adopt the approach stated in subsection IV A, the problem is
easily solved and in a consistent manner.
The time used for the algorithm was 110 seconds and this was the longest time for all the programs
run in this study. The running times for the other cases, were between 30 and 80 seconds. The value
of the objective function was f ob = 3.69 × 10−3 for the design variable vector whose components are
shown in Table VI.
x0
y0
r1
r2
r3
r4
rcx
rcy
γ
0.22922 -0.63525 2.27468 0.44667 2.18422 0.72409 1.02937 0.82440 0.58183
ψ0
ψ1
ψ2
ψ3
ψ4
ψ5
ψ6
ψ7
ψ8
0.78140 1.09985 1.34998 1.68045 2.00009 2.35036 2.70304 2.95102 3.22683
ψ9
ψ 10
ψ 11
ψ 12
ψ 13
ψ 14
ψ 15
ψ 16
ψ 17
3.58801 4.11376 4.35829 4.70801 5.07939 5.35914 5.76271 6.21586 6.49216
TABLE VI. Parameter values of an optimal mechanism. Path generation without prescribed timing.

The searching interval for the angles ψ was
0 < ψ j < 2π; j = {0, 1, ..., 17}.

(31)

It is well known that DE can yield individuals that do not belong to the searching interval. This is the
case for the last angle. However, since there is no order defect the values of table VI are an acceptable
solution for the problem.
Once again the result is obtained in two steps. First, we choose an elitist initial population that
satisfies the CG conditions. Then, the CG condition is removed and the searching space is restricted
according to the solution obtained in the first step.
It is worthwhile to mention that we tried to optimize the f ob function using the DE method without
the transformations of sections IV A and V C but the method was not capable of finding the minimum.
VII.

ACKERMAN STEERING LINKAGE SYNTHESIS

In this section DE is used for the synthesis of an Ackerman steering. For the deduction of the
equations used and a detailed treatment of the problem see [49].
It is known that when a vehicle is moving very slowly there is a kinematic condition between the
inner and outer wheel that allows it to turn slip-free. The condition is called the Ackerman condition
and is written as follows:
w
cot δo − cot δi = ,
(32)
l
where w and l represent the width and length of the vehicle, δo and δi are the rotation angles of the
wheels (Figure 5).
In general it is desirable for a mechanism to satisfy the Ackerman condition. Unfortunately, there
is no four-bar mechanism that can fulfill the Ackerman condition perfectly. However, it is possible
to synthesize a six-bar mechanism to work closely to the Ackerman condition and be exact at a few
points.
A six-bar Watt’s mechanism can be used to design the vehicle steering. The sizes In this case are
w = 1 m, l = 1.8 m and the minimum radius R = 2.5 m. The position of the center of mas with
respect to the rear axle is a = 0.45 m.
We have
p
R = a2 + l2 cot2 δM ,
(33)

12

Outer
whell

Inner
whell

Center of
rotation

FIG. 5. Vehicle diagram.

FIG. 6. Six-bar Watt’s mechanism.

with
δM =

cot δo + cot δi
2

therefore δM = 37.2731◦ , R1 = l cos δM and consequently R1 = 2.36514 m.
From trigonometry we have




l
l
δi = arctan
;
δ
=
arctan
o
R1 − w2
R1 + w2

(34)

(35)

so we obtain that δi and δ0 must lie in the ranges −32.1387◦ 6 δi 6 43.9818◦ and −43.9818◦ 6 δo 6
32.1387◦ in order to achieve the desired turning radius.
Unlike previous examples where the number of points is finite, in this case it is possible to use an
arbitrary number of desired angles. Therefore, it is convenient to define the objective function as
f ob =

|E|2
,
n

(36)

where E is the vector containing the n differences between δ2 and δack . The angle δack is the steering
angle δo from Eq. (32).

13
h

x

ξ

0.298192 -0.472091 0.219837
TABLE VII. Parameter values for the multi-link Ackerman steering.
0.5
0.4
0.3
0.2
0.1
0.5

-0.5
-0.1

FIG. 7. Optimal mechanism. Ackerman steering linkage synthesis.

The design variable vector is
X = {h, x, ξ} .

(37)

Fig. 6 shows the mechanism.
With a search space of 0.1 6 h 6 0.45, −0.5 6 x 6 0.2, 13◦ 6 ξ 6 30◦ , and a working range
(−35◦ , 45◦ ) with steps of 0.1◦ for δ1 we obtain the values shown in Table VII for the design variables.
The objective function is found to be f ob = 7.6 × 10−5 .
The obtained mechanism is illustrated in Fig. 7. Table VIII shows the values of the desired angles
and generated angles.

VIII.

CONCLUSIONS

Dimensional synthesis of mechanisms is a subject of great relevance in the field of mechanical design.
Among the great variety of optimization methods available, those that employ evolutionary algorithms
have seen an increase in use due to the excellent results that they allow.
In this work we have presented a methodology that uses differential evolution to solve the dimensional
synthesis problem of four mechanisms. With the use of a heuristic deduction, we have determined a
weight factor that allows us to solve the hybrid-tasks problem in an efficient manner. Two transformations were implemented in the differential evolution algorithm. The first one deals with the order
defect problem and was coded in the crossover part of the differential evolution algorithm. With this
transformation, the penalization approach and the use of big populations are avoided. In addition, the
chance of not finding the minimum of the objective function has disappeared as the need of discretization of the search space is also avoided. The second transformation constructs elitist populations in the
sense that their individuals satisfy the Grashof and crank conditions. Therefore, a random generation
and/or a penalization procedure are avoided, which makes this method more efficient.
Something that deserves mention is the amazing speed of convergence of the differential evolution
method which for generations as large as 80 000, the total CPU time was less than two minutes in a
single processor.

14
δ1

δack
◦

−30
−20◦
−10◦
0◦
10◦
20◦
30◦
40◦

δ2
◦

−40.471
−24.567◦
−11.069◦
0◦
9.117◦
16.822◦
23.571◦
29.720◦

−40.254◦
−23.700◦
−10.810◦
0◦
9.303◦
17.332◦
24.174◦
29.869◦

TABLE VIII. Desired (δack ) and generated (δ2 ) angles.

ACKNOWLEDGMENTS

We want to thank F. Larios for proofreading the manuscript. We also thank PROMEP and Conacyt
for support.

Appendix A: Fortran 90 objective function (f ob). Path generation for 18 target points and 10
design variables

What follows is the code for the f ob function in Fortran 90.
Double precision function fob(x0,y0,r1,r2,r3,r4,rcx,rcy,gamma,psi0)
Implicit None
Integer, Parameter:: Np=18
Double precision, Parameter :: Pi=3.14159265358979d0
Double precision :: x0,y0,r1,r2,r3,r4,rcx,rcy,gamma,psi0,L1,L2,L3,xd(Np), &
yd(Np),psi(Np),KA(Np),KB(Np),KC(Np),theta(Np),px(Np),py(Np),Ex(Np),Ey(Np),&
Ex2,Ey2
Integer :: k
xd=(/0.5d0, 0.4d0, 0.3d0, 0.2d0, 0.1d0, 0.05d0, 0.02d0, 0d0, 0d0,0.03d0,&
0.1d0, 0.15d0, 0.2d0, 0.3d0, 0.4d0, 0.5d0, 0.6d0, 0.6d0/)
yd = (/1.1d0, 1.1d0, 1.1d0, 1d0, 0.9d0, 0.75d0, 0.6d0,0.5d0,0.4d0,0.3d0,&
0.25d0, 0.2d0, 0.3d0, 0.4d0, 0.5d0, 0.7d0, 0.9d0, 1d0/)
L3=(r4**2 - r1**2 - r2**2 - r3**2)/(2d0*r2*r3)
L2=r1/r3
L1=r1/r2
Do, k=1,Np
psi(k) = psi0 + (k-1)*Pi/9d0
Enddo
KA = Dcos(psi) - L1 + L2*Dcos(psi) + L3
KB = -2d0*Dsin(psi)
KC = L1 + (L2 - 1)*Dcos(psi) + L3
theta = 2d0*Datan2(-KB - Dsqrt(KB**2-4d0*KA*KC),2d0*KA)
px = x0 + Dcos(gamma)*(r2*Dcos(psi) + rcx*Dcos(theta) - rcy*Dsin(theta)) - &

15
Dsin(gamma)*(r2*Dsin(psi) + rcx*Dsin(theta) + rcy*Dcos(theta))
py = y0 + Dsin(gamma)*(r2*Dcos(psi) + rcx*Dcos(theta) - rcy*Dsin(theta)) + &
Dcos(gamma)*(r2*Dsin(psi) + rcx*Dsin(theta) + rcy*Dcos(theta))
Ex = xd-px
Ey = yd-py
Ex2 = Dot_Product(Ex,Ex)
Ey2 = Dot_Product(Ey,Ey)
fob = Ex2 + Ey2
Return
End

Appendix B: Mathematica R hybrid task animation

nparam = {-8.0339,1.07673,13.2425,1.96639,7.71759,7.57298,13.4593,3.13037,
3.5639,3.83348,4.05641,4.22857,4.48498,4.71726,4.92507,5.83047};
vparam = {x0,y0,r1,r2,r3,r4,rcx,rcy,psi1,psi2,psi3,psi4,psi5,psi6,psi7,gamma};
supersolanima = Thread[Rule[vparam, nparam]]
r0 = {x0, y0} /. supersolanima;
xs = {7.03,6.95,6.77,6.4,5.91,5.43,4.93,4.67,4.38,4.04,3.76,3.76,3.76,3.76,
3.76,3.76,3.76,3.8,4.07,4.53,5.07,5.05,5.89,6.41,6.92};
ys = {5.99,5.45,5.03,4.6,4.03,3.56,2.94,2.6,2.2,1.657,1.22,1.97,2.78,3.56,
4.34,4.91,5.47,5.98,6.4,6.75,6.85,6.84,6.83,6.8,6.58};
DatT = Thread[{xs, ys}];
Dat = Take[DatT, {11, 25}];
L3 = (r4^2 - r1^2 - r2^2 - r3^2)/(2 r2 r3);
L2 = r1/r3;
L1 = r1/r2;
KA = Cos[psi] - L1 + L2 Cos[psi] + L3;
KB = -2 Sin[psi];
KC = L1 + (L2 - 1) Cos[psi] + L3;
theta[psi_] = 2 ArcTan[(-KB - Sqrt[KB^2 - 4 KA KC])/
(2 KA)] /.supersolanima;
Px[psi_] = (x0 + Cos[gamma] (r2 Cos[psi] + rcx Cos[theta[psi]] rcy Sin[theta[psi]]) - Sin[gamma] (r2 Sin[psi] + rcx Sin[theta[psi]] +
rcy Cos[theta[psi]])) /.supersolanima;
Py[psi_] = (y0 + Sin[gamma] (r2 Cos[psi] + rcx Cos[theta[psi]] rcy Sin[theta[psi]]) + Cos[gamma] (r2 Sin[psi] + rcx Sin[theta[psi]] +
rcy Cos[theta[psi]])) /.supersolanima;

16

PrPl[psi_] = Thread[{Px[psi], Py[psi]}];
B[psi_] = r0 + {r2 Cos[psi + gamma], r2 Sin[psi + gamma]} /.supersolanima;
Cc[psi_] = B[psi] + {r3 Cos[theta[psi] + gamma], r3 Sin[theta[psi] +
gamma]} /.supersolanima;
CoorD = (r0 + r1 {Cos[gamma], Sin[gamma]}) /. supersolanima;
gDat = ListPlot[Dat];
linkb[psi_]
linkc[psi_]
linkd[psi_]
linke[psi_]
linkf[psi_]

:=
:=
:=
:=
:=

Graphics[{Thick,
Graphics[{Thick,
Graphics[{Thick,
Graphics[{Thick,
Graphics[{Thick,

Line[{B[psi], r0}]}];
Line[{Cc[psi], B[psi]}]}];
Line[{Cc[psi], CoorD}]}];
Line[{PrPl[psi], B[psi]}]}];
Line[{PrPl[psi], Cc[psi]}]}];

linka = Graphics[{Thickness[.01], EdgeForm[Thick],RGBColor[0.75, 0.75, 0.75],
Polygon[{r0, {r0[[1]], CoorD[[2]]}, CoorD}]},PlotRange -> {{-10, 10},
{-4.8, 8}}];
gr = ListPlot[{Table[PrPl[psi], {psi, 0, 2 Pi, Pi/50}], Dat},Joined ->
{True, False}, PlotStyle ->{{PointSize[Medium],AbsoluteThickness[1.1]}},
PlotRange ->{{-10, 10}, {-5.2, 7.5}}, Frame -> True, Axes -> False,
FrameLabel -> {"X", "Y"}, AspectRatio -> Automatic];
Animate[Show[{linka, gr, linkb[psi], linkc[psi], linkd[psi], linke[psi],
linkf[psi], gDat}, Axes -> True, AspectRatio -> Automatic], {psi, 0, 2 Pi}]

[1] R. S. Hartenberg and D. J, Kinematic Synthesis of Linkages (McGraw Hill, New York, 1964).
[2] F. Freudenstein, Transactions of the ASME 76, 483 (1954).
[3] C. H. Suh and C. W. Radcliffe, Kinematics and Mechanism Desig (John Wiley and Sons, New York,
1978).
[4] G. N. Sandor and A. G. Erdman, Advanced Mechanism Design, Vol. 2 (Prentice Hall, 1984).
[5] A. P. Morgan and C. W. Wampler, Journal of Mechanical Design 112, 544 (1990).
[6] K. Al-Widyan, J. Angeles, and J. J. Cervantes-Sánchez, Proceedings of Design Engineering Technical
Conferences(September 2002).
[7] J. S. Bourrelle, C. Chen, S. Caro, and J. Angeles, Proc. 12th IFToMM World Congress, Besanon(June
18-21 2007).
[8] J. L. Blechschmidt and J. J. Uiker, Mechanisms, Transmissions, and Automation in Design 108, 543
(December 1986).
[9] H. H. Mabie and C. F. Reinholtz, Mechanisms and Dynamics of Machinery, 4th ed. (Wiley, New York,
1987).
[10] R. L. Norton, Kinematics and Dynamics of Machinery (McGraw Hill Inc., New York, 2009).
[11] J. J. Cervantes-Sánchez, H. I. Medellı́n-Castillo, J. M. Rico-Martı́nez, and E. J. González-Galván, Mechanism and Machine Theory 44, 103 (January 2009).
[12] G. S. Singh and S. Beohar, Journal of the Franklin Institute 334, 377 (1997).
[13] C. H. Suh and C. W. Radcliffe, ASME Paper , (1966).
[14] H. Nolle and K. H. Hunt, Journal of Mechanisms 6, 267 (1971).
[15] C. H. Suh and A. W. Mecklenburg, Mechanism and Machine Theory 8, 479 (Winter 1973).

17
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39]
[40]
[41]
[42]
[43]
[44]
[45]
[46]
[47]
[48]
[49]

R. I. Alizade, I. G. Novruzbekov, and G. N. Sandor, Mechanism and Machine Theory 10, 327 (1975).
J. Angeles, A. Alivizatos, and R. Akhras, Mechanism and Machine Theory 23, 343 (1988).
R. Akhras and J. Angeles, Mechanism and Machine Theory 25, 97 (1990).
Z. Liu and J. Angeles, Journal of Mechanical Design 114, 569 (1992).
S. Krishnamurty and D. A. Turcic, Mechanism and Machine Theory, 599(1992).
H. Martinez-Alfaro, H. Valdez, and J. Ortega, ASME Design Engineering Technical Conferences, (1998).
S. Krishnamurty, Proceedings of 3rd Applied Mechanisms and Robotics Conference, 94(1993).
A. Vasiliu and B. Yannou, Mechanism and Machine Theory 36, 299 (2001).
G. Galán-Marı́n, F. J. Alonso, and J. M. Del Castillo, Mechanism and Machine Theory 44, 1132 (2009).
W. E. Fang, Proceeding of 23rd Biennial Mechanisms Conference(1994).
A. Kunjur and S. Krishnamurthy, Journal of Applied Mechanisms and Robotics 4, 18 (1997).
J. Cabrera, A. Simon, and M. Prado, Mechanism and Machine Theory 37, 1165 (October 2002).
R. Starosta, Journal of Theoretical and Applied Mechanics 46, 395 (2008).
N. Nariman-Zadeh, M. Felezi, A. Jamali, and M. Ganji, Mechanism and Machine Theory 44, 180 (2009).
S. K. Acharyya and M. Mandal, Mechanism and Machine Theory 44, 1784 (September 2009).
P. S. Shiakolas, D. Koladiya, and J. Kebrle, Journal of Inverse Problems in Engineering 10, 485 (2002).
P. S. Shiakolas, D. Koladiya, and J. Kebrle, Mechanism and Machine Theory 40, 319 (March 2005).
J. Cabrera, F. Nadal, J. P. M. noz, and A. Simon, Mechanism and Machine Theory 42, 791 (July 2007).
R. R. Bulatović and S. R. Dordević, Mechanism and Machine Theory 44, 235 (January 2009).
I. Ullah and S. Kota, Journal of Mechanical Design 119, 504 (1997).
M. A. Laribi, A. Mlika, L. Romdhane, and S. Zeghloul, Mechanism and Machine Theory 39, 717 (July
2004).
A. Smaili, N. Diab, and N. Atallah, Mechanical Design 127, 917 (2005).
A. Smaili and N. Diab, Mechanism and Machine Theory 42, 115 (January 2007).
W. Y. Lin, Mechanism and Machine Theory 45, 1096 (August 2010).
J. Brest, S. Greiner, B. Boskovic, M. Mernik, and V. Zumer, IEEE Trans. on Evol. Comp. 10, 646
(December 2006).
K. V. Price, R. M. Storn, and J. A. Lampinen, Differential Evolution: A Practical Approach to Global
Optimization (Springer, Germany, 2005).
R. Storn, Biennial Conference of the North American Fuzzy Information Pro- cessing Society (NAFIPS),
IEEE, 519(1996).
J. Ronkkonen, S. Kukkonen, and K. V. Price, IEEE Congress on Evolutionary Computation, 506(2005).
R. Storn, Advances in Differential Evolution; Differential Evolution Research Trends and Open Questions, Vol. 143 (Springer, 2008) pp. 1–31.
J. R. McGarva, Mechanism and Machine Theory 29, 223 (1994).
W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran 77 The
Art of Scientific Computing, 2nd ed., Vol. 1 (Cambridge Univerity Press, 1995).
W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, The Art of Parallel Scientific Computing, 2nd ed., Vol. 2 (Cambridge Univerity Press, 2002).
W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, The Art of Scientific Computing, 2nd
ed., Vol. 2 (Cambridge Univerity Press, 2007).
R. N. Jazar, Vehicle Dynamics: Theory and Application (Springer, 2008).

