**3. Heuristic algorithms**

Optimization studies drawn attention in 1960's. It was developed a lot of algorithm based on the more sophisticated mathematical background. These algorithms were called exact methods. But, exact methods produced certain solutions only for the limited scope of application. As a result, attention turned to heuristic algorithms. (Fisher et al., 1998). Heuristics algorithms try to find acceptable solutions to the problems using some heuristic knowledge and most of them simulate real life. They use not only pure mathematics, but also algorithms using with basic formulations. While the most important property of heuristic algorithms is that they are designed for the unconstrained optimization problems, they can also be adapted to the constrained optimization problems. The algorithms are designed to find the minimum value of the objective function within the bounds of the constraints . If a solution doesn't satisfy the constraints, this solution is not acceptable, even if the value of the objective function is minimum. So, if there are constraints in a minimization problem, penalty function is added to objective function. The total function is called as fitness function. Namely, if constraints are in feasible region, then there is no penalty and penalty function is equal to zero. If the constraints are not in feasible region, the fitness function is penalized by penalty function.

Heuristic algorithms don't guarantee finding the optimal solutions. They try to find acceptable solutions near to optimum in a reasonable time. These algorithms were studied for both discrete and continuous optimization problems since the 1970's. Researchers have tried to develop adaptive and hybrid heuristics algorithms. By this aim, tens of algorithms were developed in last decade. Heuristic algorithms can be classified as single solution

Heuristic Optimization Algorithms in Robotics 315

(Kalra et al., 2006), trajectory generation for non-holonomic mobile manipulators (Chen & Zalzala, 1997), generation of walking periodic motions for a biped robot (Selene et al, 2011), the solution of the forward kinematics of 3RPR planar parallel manipulator(Chandra & Rolland, 2011), kinematic design optimization of a parallel ankle rehabilitation robot (Kumar

GA starts to run with a lot of possible solutions according to the initial population which are randomly prepared. Each individual is encoded as fixed-length binary string, characterbased or real-valued encodings. Encoding is an analogy with an actual chromosome (Lee & El-Sharkawi, 2008). A binary chromosome represents the design variables with a string containing 0s and 1s. Design variables are converted to binary arrays with a precision p. The

Where, **p** is a precision selected by the program designer. **xku** and **xkl** are respectively upper and lower bounds of the design variables **B.N** is the number of bits representing design variables. If the number of bits is 4, the binary values of the design variables are represented

Initial population is generated randomly after coding the variables. Initial population is represented as a matrix. Typical population size varies between 20 and 300. Each row of the

Strings are evaluated through iterations, called generations. During each generation, the strings are evaluated using fitness function. Fitness function measures and evaluates the coded variables in order to select the best fitness strings (Saruhan, 2006). Then, it tries to find optimum solutions by using genetic operators. By this way, the best solutions are selected and worsts are eliminated. GA runs according to unconstrained optimization procedure. So, the constrained continuous optimization problems are transformed into unconstrained continuous optimization problem by penalizing the objective function value with the

Fitness function (FF) values are calculated for each individual as it is seen in equation 3, 4 and 5. In a minimization problem, penalty function is added to objective function. The fitness function (FF), is the sum of objective function (OF) and the penalty function (PF) consisting constraint functions (CF). After elitism, selection, crossover, and mutation are operated, a new population is generated according to the fitness function values. The future

OF= f (x ,x ,...x ) obj 1 2 n (3)

l m 2 2 PF R (max[0,g (x ,x ,...x ,a)]) Rj(max[0, h (x ,x ,...x , b) tol]) i i12 n j12 n i 1 j 1 

 FF=OF+PF (5) In the above equations, **a** and **b** are constant values**, l** is the number of inequality constraints, **m** is the number of equality constraints, **n** is the number of variables and **tol** is the tolerance value, for the equality constraints. If **hj**'s value is smaller than the tolerance value, it is accepted that equality constraints have been satisfied. If the problem variables satisfy the

quadratic penalty function. The total function is called fitness function.

)/p)+1) (2)

(4)

B. N≥log2(((xku -Xkl

bit number is calculated by equation 2 (Lin & Hajela, 1992).

et al., 2011).

as 0000, 0001, 0010, …, 1111.

population matrix is called as an individual.

of population depends on the evolutionary rules.

algorithms and population based algorithms. The first category gathers Local Search (LS), Greedy Heuristic (GH) (Feo et al., 1995), Simulated Annealing (SA) (Kirkpatrick et al., 1983), Tabu Search (TS) (Laguna, 1994) and Iterated Local Search (ILS) (Glover & Garry, 2002). The second category, which is more and more studied, groups evolutionary algorithms such as Genetic Algorithms (GA)(Goldberg, 1989), Ant Colony Optimization (ACO) (Dorigo et al, 1996), Artificial Bee Colony (ABC)( Basturk & Karaboga, 2006), Scatter Search (SS) (Marti et al., 2006), Immune Systems (IS) (Farmer et al., 1986), Differential Evolution Algorithms (DEA) (Storn, 1996), Particle Swarm Optimization (PSO) (Keneddy & Eberhart, 1995), Harmony Search (HS) (Geem & Kim, 2001) and Gravitational Search Algorithm (GSA) (Rashedi et al., 2009). Evolutionary algorithms simulate natural phenomena or social behaviors of animals. So, the algorithms can be subcategorized to these criteria. While SA, GSA, HS simulates natural phenomena, ACO, ABC, PSO simulates the social behaviors of animals. The taxonomy of global optimization algorithms can be found in the related literature (Weise, 2008).

During the last years, hybrid optimization approaches have been popular. Firstly, cooperations have been realized between several heuristic algorithms. Later, cooperations between heuristic algorithms and exact methods have been realized. Since they gather the advantages of both types of algorithm, hybrid algorithms give more satisfied results. (Jourdan et al., 2009).

In single solution heuristic algorithms, the algorithm starts with an initial solution called current solution. This solution can be created randomly and consists of the values of design variables of the optimization problems. The objective function value is calculated with these initial values. If the initial solution satisfies the constraints, fitness function will be the same with the objective function. The second solution, called candidate solution, is created with random values close to the initial solution. Between the two fitness functions, the one which has a minimum value is selected as candidate solution. As long as the iteration proceeds, it is expected that algorithm converges to the global minimum value. But algorithms can be trapped to the local minimums. Heuristic algorithms develop different strategies to avoid to get trapped the local minimums. So, exploration and exploitation are of great importance for finding the global minimum value with the heuristics algorithms. In the first iterations, the algorithm searches the global solution space with big steps in order not to get trapped the local minimums. This is called exploration. In the next iterations, the algorithm searches the solution space with small steps in order to find best minimum. This is called exploitation. It is desired to have a balance between exploration and exploitation in the heuristic algorithms.

## **3.1 Genetic algorithms**

GA was developed in Michigan University by Holland and later it got popularity with the efforts of Goldberg. GA is an adaptive global search method that mimics the metaphor of natural biological evolution. It operates on a population of potential solutions by applying the principle of survival of the fittest to achieve an optimal solution.

In the literature, GA is one of the most applied heuristic optimization algorithms. GA has been applied successfully to a huge variety of optimization problems, nearly in all disciplines, such as materials science, aircraft applications, chemistry, construction, seismology, medicine and web applications. GA has also been applied to the problems in Robotic such as the solution of multimodal inverse kinematics problem of industrial robots

algorithms and population based algorithms. The first category gathers Local Search (LS), Greedy Heuristic (GH) (Feo et al., 1995), Simulated Annealing (SA) (Kirkpatrick et al., 1983), Tabu Search (TS) (Laguna, 1994) and Iterated Local Search (ILS) (Glover & Garry, 2002). The second category, which is more and more studied, groups evolutionary algorithms such as Genetic Algorithms (GA)(Goldberg, 1989), Ant Colony Optimization (ACO) (Dorigo et al, 1996), Artificial Bee Colony (ABC)( Basturk & Karaboga, 2006), Scatter Search (SS) (Marti et al., 2006), Immune Systems (IS) (Farmer et al., 1986), Differential Evolution Algorithms (DEA) (Storn, 1996), Particle Swarm Optimization (PSO) (Keneddy & Eberhart, 1995), Harmony Search (HS) (Geem & Kim, 2001) and Gravitational Search Algorithm (GSA) (Rashedi et al., 2009). Evolutionary algorithms simulate natural phenomena or social behaviors of animals. So, the algorithms can be subcategorized to these criteria. While SA, GSA, HS simulates natural phenomena, ACO, ABC, PSO simulates the social behaviors of animals. The taxonomy of global optimization algorithms can be found in the related

During the last years, hybrid optimization approaches have been popular. Firstly, cooperations have been realized between several heuristic algorithms. Later, cooperations between heuristic algorithms and exact methods have been realized. Since they gather the advantages of both types of algorithm, hybrid algorithms give more satisfied results.

In single solution heuristic algorithms, the algorithm starts with an initial solution called current solution. This solution can be created randomly and consists of the values of design variables of the optimization problems. The objective function value is calculated with these initial values. If the initial solution satisfies the constraints, fitness function will be the same with the objective function. The second solution, called candidate solution, is created with random values close to the initial solution. Between the two fitness functions, the one which has a minimum value is selected as candidate solution. As long as the iteration proceeds, it is expected that algorithm converges to the global minimum value. But algorithms can be trapped to the local minimums. Heuristic algorithms develop different strategies to avoid to get trapped the local minimums. So, exploration and exploitation are of great importance for finding the global minimum value with the heuristics algorithms. In the first iterations, the algorithm searches the global solution space with big steps in order not to get trapped the local minimums. This is called exploration. In the next iterations, the algorithm searches the solution space with small steps in order to find best minimum. This is called exploitation. It is desired to have a balance between exploration and exploitation in the heuristic

GA was developed in Michigan University by Holland and later it got popularity with the efforts of Goldberg. GA is an adaptive global search method that mimics the metaphor of natural biological evolution. It operates on a population of potential solutions by applying

In the literature, GA is one of the most applied heuristic optimization algorithms. GA has been applied successfully to a huge variety of optimization problems, nearly in all disciplines, such as materials science, aircraft applications, chemistry, construction, seismology, medicine and web applications. GA has also been applied to the problems in Robotic such as the solution of multimodal inverse kinematics problem of industrial robots

the principle of survival of the fittest to achieve an optimal solution.

literature (Weise, 2008).

(Jourdan et al., 2009).

algorithms.

**3.1 Genetic algorithms** 

(Kalra et al., 2006), trajectory generation for non-holonomic mobile manipulators (Chen & Zalzala, 1997), generation of walking periodic motions for a biped robot (Selene et al, 2011), the solution of the forward kinematics of 3RPR planar parallel manipulator(Chandra & Rolland, 2011), kinematic design optimization of a parallel ankle rehabilitation robot (Kumar et al., 2011).

GA starts to run with a lot of possible solutions according to the initial population which are randomly prepared. Each individual is encoded as fixed-length binary string, characterbased or real-valued encodings. Encoding is an analogy with an actual chromosome (Lee & El-Sharkawi, 2008). A binary chromosome represents the design variables with a string containing 0s and 1s. Design variables are converted to binary arrays with a precision p. The bit number is calculated by equation 2 (Lin & Hajela, 1992).

$$\text{B. N2} \text{log}\_2(((\text{x}\_k \text{ -} \text{X}\_k \text{ })/\text{p}) + 1) \tag{2}$$

Where, **p** is a precision selected by the program designer. **xku** and **xkl** are respectively upper and lower bounds of the design variables **B.N** is the number of bits representing design variables. If the number of bits is 4, the binary values of the design variables are represented as 0000, 0001, 0010, …, 1111.

Initial population is generated randomly after coding the variables. Initial population is represented as a matrix. Typical population size varies between 20 and 300. Each row of the population matrix is called as an individual.

Strings are evaluated through iterations, called generations. During each generation, the strings are evaluated using fitness function. Fitness function measures and evaluates the coded variables in order to select the best fitness strings (Saruhan, 2006). Then, it tries to find optimum solutions by using genetic operators. By this way, the best solutions are selected and worsts are eliminated. GA runs according to unconstrained optimization procedure. So, the constrained continuous optimization problems are transformed into unconstrained continuous optimization problem by penalizing the objective function value with the quadratic penalty function. The total function is called fitness function.

Fitness function (FF) values are calculated for each individual as it is seen in equation 3, 4 and 5. In a minimization problem, penalty function is added to objective function. The fitness function (FF), is the sum of objective function (OF) and the penalty function (PF) consisting constraint functions (CF). After elitism, selection, crossover, and mutation are operated, a new population is generated according to the fitness function values. The future of population depends on the evolutionary rules.

$$\mathbf{OF} = \mathbf{f}\_{\mathbf{obj}}(\mathbf{x}\_1, \mathbf{x}\_2, \dots \mathbf{x}\_{\mathbf{n}}) \tag{3}$$

$$\text{PF} = \sum\_{\mathbf{i}=1}^{\text{L}} \text{R}\_{\mathbf{i}}(\max[0, \mathbf{g}\_{\mathbf{i}}(\mathbf{x}\_{1}, \mathbf{x}\_{2}, \dots \mathbf{x}\_{\mathbf{n}}, \mathbf{a})])^{2} + \sum\_{\mathbf{j}=1}^{\text{m}} \text{R}(\max[0, \left|\mathbf{h}\_{\mathbf{j}}(\mathbf{x}\_{1}, \mathbf{x}\_{2}, \dots \mathbf{x}\_{\mathbf{n}}, \mathbf{b})\right| - \text{tol}))^{2} \tag{4}$$

$$\text{FF} = \text{OC} + \text{PF} \tag{5}$$

In the above equations, **a** and **b** are constant values**, l** is the number of inequality constraints, **m** is the number of equality constraints, **n** is the number of variables and **tol** is the tolerance value, for the equality constraints. If **hj**'s value is smaller than the tolerance value, it is accepted that equality constraints have been satisfied. If the problem variables satisfy the

Heuristic Optimization Algorithms in Robotics 317

et al, 1983). Because of the fact that SA is an effective algorithm and it is easy to implement for the difficult engineering optimization problems, it is also popular for last studies in several different engineering areas. SA is applied to various optimization problem in Robotics such as path planning problem solution(Gao & Tian,2007), the generation of the assembly sequences in robotic assembly(Hong & Cho, 1999), trajectory optimization of advanced launch system (Karsli & Tekinalp, 2005), optimal robot arm PID control, the placement of serial robot manipulator and collision avoidance planning in multi-robot. SA combines local search and Metropolis Algorithm. Algorithm starts with a current solution at initial temperature. At each temperature, algorithm iterates n times, defined by the program designer. In the iterations for each temperature, a candidate solution that is close to the current solution is produced randomly. If candidate solution is better than the current solution, the candidate solution is replaced with the current solution. Otherwise, the candidate solution is not rejected at once. The random number is produced between 0-1 and compared with the acceptance probability P. P is an exponential function as given by equation 6. If produced random number is smaller than P, the worse solution is accepted as current solution for exploration. Search process progresses until stopping criteria is satisfied. Maximum run time, maximum iteration number, last temperature e.g. may be the stopping

f(x\_cand) f(x\_curr) ( )

(6)

e T

For an object function with one design variable, **x\_cand** is candidate design variable, **x\_curr** is current design variable and **f** is object function, **T** is temperature. At high temperatures, P is close to one. Since the most random number is smaller than one, the possibility of bad solution's acceptance is high for exploration in the first steps of the algorithm. T is decreased along the search process. At low temperatures, P is close to zero. Since the random numbers are bigger than zero, the possibility of bad solution's acceptance approach to zero. When P is equal to zero, the process converges to local search method for exploitation. The pseudo

P =

criteria.

code of SA is given in Figure 4.

Fig. 4. The pseudo codes of SA Algorithm

constraints, gi and (hj –tol) values will be negative and PF will be zero. **Ri** and **Rj** are the penalty coefficients and these values affect the algorithm performance directly.

Selection is to choose the individuals for the generation based on the principle of survival of the best according to the Darwin's theory. The main idea in selection is to create parents from the chromosomes having better fitness function values. Namely the parents are selected from the chromosomes of the first population as directly proportional to their fitness function values. The better fitness function value of a chromosome, the more chance to be a parent. With the selection method, the individuals having worse fitness values disappear in the next generations. The most popular selection methods are roulette wheel and tournament rank. The successful individuals are called parents.

After selection, crossover, which represents mating of two individuals, is applied to the parents. The parents in GA create two children with crossover. There are various kinds of crossover process such as one point crossover, multiple point crossovers and uniform crossover. The crossover operator is applied with a certain probability, usually in the range 0.5-1. This operator allows the algorithm to search the solution space. Namely, exploration is performed in GA with crossover operator. In a basic crossover, the children called offspring are created from mother and father genes with a given crossover probability. Mutation is an another operation in GA and some gens are changed with this operation. But this is not common. Generally mutation rate of binary encoding is specified smaller than 0.1. In contrast to the crossover, mutation is used for the exploitation of the solution space. The last operation is the elitism. It transfers the best individual to the next generation. The evaluation of the previous population with the operators stated above, the new population is generated till the number of generation. Fitness function values are calculated in each new population and the best resulted ones are paid attention among these values. Until the stopping criteria are obtained, this process is repeated iteratively. The stopping criteria may be the running time of the algorithm, the number of generation and for fitness functions giving the same best possible values in a specified time. The pseudo codes of GA are given in Figure 3.


Fig. 3. The pseudo codes of GA

#### **3.2 Simulated annealing**

SA uses an analogy of physical annealing process of finding low energy states of solids and uses global optimization method. SA is inspired by Metropolis Algorithm and this approach was firstly submitted to an optimization by Kirkpatrick and his friends in 1983 (Kirkpatrick

constraints, gi and (hj –tol) values will be negative and PF will be zero. **Ri** and **Rj** are the

Selection is to choose the individuals for the generation based on the principle of survival of the best according to the Darwin's theory. The main idea in selection is to create parents from the chromosomes having better fitness function values. Namely the parents are selected from the chromosomes of the first population as directly proportional to their fitness function values. The better fitness function value of a chromosome, the more chance to be a parent. With the selection method, the individuals having worse fitness values disappear in the next generations. The most popular selection methods are roulette wheel

After selection, crossover, which represents mating of two individuals, is applied to the parents. The parents in GA create two children with crossover. There are various kinds of crossover process such as one point crossover, multiple point crossovers and uniform crossover. The crossover operator is applied with a certain probability, usually in the range 0.5-1. This operator allows the algorithm to search the solution space. Namely, exploration is performed in GA with crossover operator. In a basic crossover, the children called offspring are created from mother and father genes with a given crossover probability. Mutation is an another operation in GA and some gens are changed with this operation. But this is not common. Generally mutation rate of binary encoding is specified smaller than 0.1. In contrast to the crossover, mutation is used for the exploitation of the solution space. The last operation is the elitism. It transfers the best individual to the next generation. The evaluation of the previous population with the operators stated above, the new population is generated till the number of generation. Fitness function values are calculated in each new population and the best resulted ones are paid attention among these values. Until the stopping criteria are obtained, this process is repeated iteratively. The stopping criteria may be the running time of the algorithm, the number of generation and for fitness functions giving the same best possible values in a specified time. The pseudo codes of GA are given

SA uses an analogy of physical annealing process of finding low energy states of solids and uses global optimization method. SA is inspired by Metropolis Algorithm and this approach was firstly submitted to an optimization by Kirkpatrick and his friends in 1983 (Kirkpatrick

penalty coefficients and these values affect the algorithm performance directly.

and tournament rank. The successful individuals are called parents.

in Figure 3.

Fig. 3. The pseudo codes of GA

**3.2 Simulated annealing** 

et al, 1983). Because of the fact that SA is an effective algorithm and it is easy to implement for the difficult engineering optimization problems, it is also popular for last studies in several different engineering areas. SA is applied to various optimization problem in Robotics such as path planning problem solution(Gao & Tian,2007), the generation of the assembly sequences in robotic assembly(Hong & Cho, 1999), trajectory optimization of advanced launch system (Karsli & Tekinalp, 2005), optimal robot arm PID control, the placement of serial robot manipulator and collision avoidance planning in multi-robot.

SA combines local search and Metropolis Algorithm. Algorithm starts with a current solution at initial temperature. At each temperature, algorithm iterates n times, defined by the program designer. In the iterations for each temperature, a candidate solution that is close to the current solution is produced randomly. If candidate solution is better than the current solution, the candidate solution is replaced with the current solution. Otherwise, the candidate solution is not rejected at once. The random number is produced between 0-1 and compared with the acceptance probability P. P is an exponential function as given by equation 6. If produced random number is smaller than P, the worse solution is accepted as current solution for exploration. Search process progresses until stopping criteria is satisfied. Maximum run time, maximum iteration number, last temperature e.g. may be the stopping criteria.

$$\mathbf{P} = \frac{-(\frac{\mathbf{f(x\_c)} \mathbf{cand}}{\mathbf{T}} - \mathbf{f(x\_c)} \mathbf{r})}{\mathbf{T}} \tag{6}$$

For an object function with one design variable, **x\_cand** is candidate design variable, **x\_curr** is current design variable and **f** is object function, **T** is temperature. At high temperatures, P is close to one. Since the most random number is smaller than one, the possibility of bad solution's acceptance is high for exploration in the first steps of the algorithm. T is decreased along the search process. At low temperatures, P is close to zero. Since the random numbers are bigger than zero, the possibility of bad solution's acceptance approach to zero. When P is equal to zero, the process converges to local search method for exploitation. The pseudo code of SA is given in Figure 4.

Fig. 4. The pseudo codes of SA Algorithm

Annealing process simulates optimization. Annealing is realized slowly in order to keep the system of the melt in a thermodynamic equilibrium. Convergence to the optimal solution is controlled by the annealing process.

Proper cooling scheme is important for the performance of SA. The proper annealing process is related with the initial temperature, iteration for each temperature, temperature decrement coefficient and stopping criteria. All these criteria can be found in related article (Blum & Roli, 2001). In general, the temperature is updated according to equation 7.

$$\mathbf{T}\_{k+1} = \alpha \mathbf{T}\_k \tag{7}$$

Heuristic Optimization Algorithms in Robotics 319

current velocity and the next velocity of the *i*th particle. The next position of a particle is

PSO simulates social swarm's behaviour. The velocity formula consists of three important components. First part presents the past velocity. The particle generally continues in the

thinking. This part is also called self-knowledge and the last part **2rand()(gbest-xi**

presents the social part of the particle swarm(Kennedy, 1997). When implementing the particle swarm algorithm, some considerations must be paid attention in order to facilitate the convergence and prevent going to infinity of the swarm. These considerations are limiting the maximum velocity, selecting acceleration constants, the constriction factor, or

The velocity and position formula with constriction factor (K) has been introduced by Clerc(Clerc, 1999) and Eberhart and Shi (Eberhart & Shi, 2001) as given by equation 10,11

k1 k k k <sup>k</sup> v K(v rand()(p x ) rand()(g x )) i i1 ii 2 best best <sup>i</sup>

k1 k k1 x xv i ii

<sup>2</sup> K , , 4 1 2 <sup>2</sup> 2 4

Where, **K** is the constriction factor, **1** and **2** represent the cognitive and social parameters, respectively. **rand** is uniformly distributed random number. is set to 4.1. The constriction factor produces a damping effect on the amplitude of an individual particle's oscillations, and as a result, the particle will converge over time. The pseudo code of PSO is given in Figure 5.

GSA was introduced by Esmat Rashedi and his friends in 2009 as an evolutionary algorithm. It is one of the recent optimization techniques inspired by the law of gravity. GSA was

   (10)

> 

**<sup>k</sup>** is the best value of *i*th

**k+1** are respectively the

**k+1** are respectively the

**k)** presents private

(12)

**k)**

**<sup>k</sup>** and **xi**

**<sup>k</sup>** and **vi**

(11)

**k-xi**

Where, **k** is the iteration number, **rand** is a random number, **pbesti**

same direction. This is called habit. The second part **1rand()(pbesti**

particle and **gbest** is the global best value of all particles. **xi**

current position and the next position of the *i*th particle. **vi**

specified by its current position and its velocity.

the inertia constant (Valle et al., 2008).

Fig. 5. The pseudo codes of PSO

**3.4 Gravitational search algorithm** 

and 12.

Where, **Tk+1** is the next temperature, **Tk** is the current temperature and is the temperature decrement coefficient, generally selected between 0.80 and 0.99.

Even if SA tries to find global minimum, it can be trapped by local minima. So, in order to overcome this drawback and develop the performance of SA, researchers have developed adaptive or hybridized SA with other heuristics, such as fuzzy logic, genetic algorithms, support vector machines, and distributed/parallel algorithms.

#### **3.3 Particle swarm optimization**

PSO was introduced by James Kennedy and Russell Eberhart in 1995 as an evolutionary computation technique (Kennedy & Eberhart, 1995) inspired by the metaphor of social interaction observed among insects or animals. It is simpler than GA because PSO has no evolution operators such as crossover and mutation (Lazinca, 2009).

PSO have many advantages such as simplicity, rapid convergence, a few parameters to be adjusted and no operators (Li & Xiao; 2008). So, PSO is used for solving discrete and continuous problems. It was applied to a wide range of applications such as function optimization, Electrical power system applications, neural network training, task assignment and scheduling problems in operations research, fuzzy system control and pattern identification. PSO was also applied to a lot of different Robotic applications. PSO was applied mobile robot navigation (Gueaieb & Miah, 2008), Robot Path Planning (Zargar & Javadi, 2009) and Swarm Robotic Applications and Robot Manipulators applications.

PSO algorithm is initialized with a group of particles. Each particle is a position vector in the search space and its dimension represents the number of design variables. PSO is started with random initial positions and velocities belong to each particle. For each particle, fitness function's value is computed. Global and local best values are updated. Local best value is the best value of the current particle found so far and the global best value is the best value of all local bests found so far. Velocity and positions are updated according to the global best and local bests. Each particle flies through the search space, according to the local and global best values. Along the iterations, particles converge to the global best. The convergence of the particles to the best solution shows the PSO performance. The rate of position change of the *i*th particle is given by its velocity. In the literature there are different velocity formulas. The first one proposed by Keneddy is given by equation 9. Particles' velocity and positions are updated according to equation 8 and 9.

$$\mathbf{x}\_1^{\mathbf{k}+1} = \mathbf{x}\_1^{\mathbf{k}} + \mathbf{v}\_1^{\mathbf{k}+1} \tag{8}$$

$$\mathbf{v\_{i}^{k+1}} = \mathbf{v\_{i}^{k}} + \boldsymbol{\varphi\_{1}} \text{rand()} (\mathbf{p\_{besti}^{k}} - \mathbf{x\_{i}^{k}}) + \boldsymbol{\varphi\_{2}} \text{rand()} (\mathbf{g\_{best}} - \mathbf{x\_{i}^{k}}) \tag{9}$$

Annealing process simulates optimization. Annealing is realized slowly in order to keep the system of the melt in a thermodynamic equilibrium. Convergence to the optimal solution is

Proper cooling scheme is important for the performance of SA. The proper annealing process is related with the initial temperature, iteration for each temperature, temperature decrement coefficient and stopping criteria. All these criteria can be found in related article

 Tk+1 = T k (7) Where, **Tk+1** is the next temperature, **Tk** is the current temperature and is the temperature

Even if SA tries to find global minimum, it can be trapped by local minima. So, in order to overcome this drawback and develop the performance of SA, researchers have developed adaptive or hybridized SA with other heuristics, such as fuzzy logic, genetic algorithms,

PSO was introduced by James Kennedy and Russell Eberhart in 1995 as an evolutionary computation technique (Kennedy & Eberhart, 1995) inspired by the metaphor of social interaction observed among insects or animals. It is simpler than GA because PSO has no

PSO have many advantages such as simplicity, rapid convergence, a few parameters to be adjusted and no operators (Li & Xiao; 2008). So, PSO is used for solving discrete and continuous problems. It was applied to a wide range of applications such as function optimization, Electrical power system applications, neural network training, task assignment and scheduling problems in operations research, fuzzy system control and pattern identification. PSO was also applied to a lot of different Robotic applications. PSO was applied mobile robot navigation (Gueaieb & Miah, 2008), Robot Path Planning (Zargar & Javadi, 2009) and Swarm Robotic Applications and Robot Manipulators applications. PSO algorithm is initialized with a group of particles. Each particle is a position vector in the search space and its dimension represents the number of design variables. PSO is started with random initial positions and velocities belong to each particle. For each particle, fitness function's value is computed. Global and local best values are updated. Local best value is the best value of the current particle found so far and the global best value is the best value of all local bests found so far. Velocity and positions are updated according to the global best and local bests. Each particle flies through the search space, according to the local and global best values. Along the iterations, particles converge to the global best. The convergence of the particles to the best solution shows the PSO performance. The rate of position change of the *i*th particle is given by its velocity. In the literature there are different velocity formulas. The first one proposed by Keneddy is given by equation 9. Particles'

k1 k k1 x xv i ii

k1 k k k <sup>k</sup> v v rand()(p x ) rand()(g x )) i i1 ii 2 best best <sup>i</sup>

 (9)

(8)

(Blum & Roli, 2001). In general, the temperature is updated according to equation 7.

decrement coefficient, generally selected between 0.80 and 0.99.

support vector machines, and distributed/parallel algorithms.

evolution operators such as crossover and mutation (Lazinca, 2009).

velocity and positions are updated according to equation 8 and 9.

controlled by the annealing process.

**3.3 Particle swarm optimization** 

Where, **k** is the iteration number, **rand** is a random number, **pbesti <sup>k</sup>** is the best value of *i*th particle and **gbest** is the global best value of all particles. **xi <sup>k</sup>** and **xi k+1** are respectively the current position and the next position of the *i*th particle. **vi <sup>k</sup>** and **vi k+1** are respectively the current velocity and the next velocity of the *i*th particle. The next position of a particle is specified by its current position and its velocity.

PSO simulates social swarm's behaviour. The velocity formula consists of three important components. First part presents the past velocity. The particle generally continues in the same direction. This is called habit. The second part **1rand()(pbesti k-xi k)** presents private thinking. This part is also called self-knowledge and the last part **2rand()(gbest-xi k)** presents the social part of the particle swarm(Kennedy, 1997). When implementing the particle swarm algorithm, some considerations must be paid attention in order to facilitate the convergence and prevent going to infinity of the swarm. These considerations are limiting the maximum velocity, selecting acceleration constants, the constriction factor, or the inertia constant (Valle et al., 2008).

The velocity and position formula with constriction factor (K) has been introduced by Clerc(Clerc, 1999) and Eberhart and Shi (Eberhart & Shi, 2001) as given by equation 10,11 and 12.

$$\mathbf{v\_{i}^{k+1}} = \mathbf{K(v\_{i}^{k} + \rho\_{1}\mathbf{rand(l)}(\mathbf{p\_{besti}^{k} - x\_{i}^{k}) + \rho\_{2}\mathbf{rand(l)}(\mathbf{g\_{best} - x\_{i}^{k}))} \tag{10}$$

$$\mathbf{x\_{i}^{k+1}} = \mathbf{x\_{i}^{k}} + \mathbf{v\_{i}^{k+1}} \tag{11}$$

$$\mathbf{K} = \frac{2}{\left| 2 - \wp - \sqrt{\wp^2 - 4\wp} \right|} \quad \text{{ $\wp = \wp\_1 + \wp\_2$ } \text{  $\wp > 4$ }} \tag{12}$$

Where, **K** is the constriction factor, **1** and **2** represent the cognitive and social parameters, respectively. **rand** is uniformly distributed random number. is set to 4.1. The constriction factor produces a damping effect on the amplitude of an individual particle's oscillations, and as a result, the particle will converge over time. The pseudo code of PSO is given in Figure 5.


Fig. 5. The pseudo codes of PSO

#### **3.4 Gravitational search algorithm**

GSA was introduced by Esmat Rashedi and his friends in 2009 as an evolutionary algorithm. It is one of the recent optimization techniques inspired by the law of gravity. GSA was originally designed for solving continuous problems. Since it is quite newly developed algorithm, it hasn't been applied to a lot of area. GSA was applied to filter modelling (Rashedi et al., 2011), Post-Outage Bus Voltage Magnitude Calculations (Ceylan et al., 2010), clustering (Yin et al., 2011), Scheduling in Grid Computing Systems (Barzegar et al., 2009).

GSA algorithm is based on the Newtonian gravity law: ''Every particle in the universe attracts every other particle with a force that is directly proportional to the product of their masses and inversely proportional to the square of the distance between them" (Rashedi et al., 2009).

In GSA, particles are considered as objects and their performances are measured by their masses. GSA is based on the two important formulas about Newton Gravity Laws given by the equation 13 and 14. The first one is as given by the equation 13. This equation is the gravitational force equation between the two particles, which is directly proportional to their masses and inversely proportional to the square of distance between them. But in GSA instead of the square of the distance, only the distance is used. The second one is the equation of accelaration of a particle when a force is applied to it. It is given by equation 14.

$$\mathbf{F} = \mathbf{G} \frac{\mathbf{M}\_1 \mathbf{M}\_2}{\mathbf{R}^2} \tag{13}$$

$$\mathbf{a} = \frac{\mathbf{F}}{\mathbf{M}} \tag{14}$$

Heuristic Optimization Algorithms in Robotics 321

Mai = Mpi = Mii = Mi, i= 1, 2, · · · ,N (16)

(17)

(18)

(19)

(20)

(21)

(22)

d dd v (t 1) rand v (t) a (t) i ii i (23)

fitness (t) worst(t) m (t) <sup>i</sup> <sup>i</sup> best(t) worst(t)

> m (t) M (t) <sup>i</sup> i N

In the equation 17, *fitnessi(t)* represents the fitness value of the *i*th mass at time t, *worst(t)* is the maximum value of all the fitness values and *best(t)* is the minimum value of the all fitness values for a minimization problem. For maximization problems, the opposite of this expression is applied. In GSA, active gravitational mass **Ma**, passive gravitational mass **Mp**

Gravitiational constant which is propotional with the time is also updated according to the

G(t)=G e <sup>T</sup> <sup>0</sup>

Where, **α** is a user specified constant, **T** is the total number of iterations, and **t** is the current iteration. This equation is similar to the SA acceptance probability equation given by the equation 6. According to the formula the gravity force is decreasing by the time. The

> M (t)M (t) <sup>d</sup> pi aj d d F (t) G(t) (x (t) x (t)) ij j i R (t) <sup>ε</sup> ij

> > d d N F (t) rand F (t) <sup>i</sup> <sup>j</sup> ij j Kbest,j i

**Fijd(t)** is the gravity force acting on mass *i* from mass *j* at time t, **G(t)** is the gravitational constant at time t, **Mpi** is the passive gravitational mass of object i, **Maj** is the active gravitational mass of object j, is a small constant, **Rij** is the Euclidean distance between two

number between 0 and 1, **Kbest** is the set of first K objects with the best fitness value and

The new positions and velocities of the particles are calculated by the equation 23 and 24.

d <sup>d</sup> F (t) <sup>i</sup> a (t) <sup>i</sup> M (t) <sup>i</sup>

**d(t)** is the force that acts on particle *i* in dimension *d*, **randj** is a random

 

biggest mass.The acceleration of an object *i* in *d*th dimension is as follows.

t α

 

and inertial mass **Mi** are accepted as they are equal to each other.

graviational force and total forces are updated as follows.

iteration number given by the equation 19.

objects i and j . **Fi**

m (t) <sup>j</sup> j 1

**G** is gravitational constant, **M1** and **M2** are masses and **R** is distance, **F** is gravitational force, **a** is acceleration. Based on these formulas, the heavier object with more gravity force attracts the other objects as it is seen in Figure 6.

Fig. 6. The Newton Gravitational Force Representation

GSA algorithm initializes the position of the masses(X) given by equation 15.

$$\mathbf{X}\_{\mathbf{i}} \equiv (\mathbf{x}\_{\mathbf{i}}^1, \mathbf{x}\_{\mathbf{i}}^2, \dots, \mathbf{x}\_{\mathbf{i}}^d, \mathbf{x}\_{\mathbf{i}}^n) \quad \mathbf{i} = \mathbf{1}, \mathbf{2}, \mathbf{3}, \dots, \mathbf{N} \tag{15}$$

Where, **Xi** is the position of *i*th mass. **n** gives the dimension of the masses. N gives the number of particles. **xi <sup>d</sup>** represents *i*th particle position in the *d*th dimension.

The first velocities of the masses are accepted as zero. Until the stopping criterion (maximum iteration), algorithm evaluates the fitness functions of the objects. After the best and worst values are found for the current iteration, the masses of the particles are updated according to the equation 16, 17 and 18.

originally designed for solving continuous problems. Since it is quite newly developed algorithm, it hasn't been applied to a lot of area. GSA was applied to filter modelling (Rashedi et al., 2011), Post-Outage Bus Voltage Magnitude Calculations (Ceylan et al., 2010), clustering (Yin et al., 2011), Scheduling in Grid Computing Systems (Barzegar et al., 2009). GSA algorithm is based on the Newtonian gravity law: ''Every particle in the universe attracts every other particle with a force that is directly proportional to the product of their masses and inversely proportional to the square of the distance between them" (Rashedi et

In GSA, particles are considered as objects and their performances are measured by their masses. GSA is based on the two important formulas about Newton Gravity Laws given by the equation 13 and 14. The first one is as given by the equation 13. This equation is the gravitational force equation between the two particles, which is directly proportional to their masses and inversely proportional to the square of distance between them. But in GSA instead of the square of the distance, only the distance is used. The second one is the equation of accelaration of a particle when a force is applied to it. It is given by equation 14.

M M1 2 F G <sup>2</sup> <sup>R</sup>

F a

**G** is gravitational constant, **M1** and **M2** are masses and **R** is distance, **F** is gravitational force, **a** is acceleration. Based on these formulas, the heavier object with more gravity force attracts

(13)

M (14)

n,) i=1,2,3,…,N (15)

al., 2009).

the other objects as it is seen in Figure 6.

Xi=(xi

according to the equation 16, 17 and 18.

number of particles. **xi**

Fig. 6. The Newton Gravitational Force Representation

GSA algorithm initializes the position of the masses(X) given by equation 15.

2, … , xi

d, , xi

Where, **Xi** is the position of *i*th mass. **n** gives the dimension of the masses. N gives the

**<sup>d</sup>** represents *i*th particle position in the *d*th dimension. The first velocities of the masses are accepted as zero. Until the stopping criterion (maximum iteration), algorithm evaluates the fitness functions of the objects. After the best and worst values are found for the current iteration, the masses of the particles are updated

1, xi

$$\mathbf{M}\_{\rm ui} = \mathbf{M}\_{\rm pi} = \mathbf{M}\_{\rm ii} = \mathbf{M}\_{\rm i} \text{ i} = \mathbf{1} \text{ } \mathbf{2}, \text{ } \cdot \text{ } \cdot \text{ } \cdot \text{ } \cdot \text{ N} \tag{16}$$

$$\mathbf{m\_{i}(t)} = \frac{\text{fitness}\_{i}(t) - \text{worst(t)}}{\text{best(t)} - \text{worst(t)}} \tag{17}$$

$$\mathbf{M}\_{\mathbf{i}}(\mathbf{t}) = \frac{\mathbf{m}\_{\mathbf{i}}(\mathbf{t})}{\mathbf{N}} \tag{18}$$
 
$$\sum\_{\mathbf{j}=1}^{\sum} \mathbf{m}\_{\mathbf{j}}(\mathbf{t})$$

In the equation 17, *fitnessi(t)* represents the fitness value of the *i*th mass at time t, *worst(t)* is the maximum value of all the fitness values and *best(t)* is the minimum value of the all fitness values for a minimization problem. For maximization problems, the opposite of this expression is applied. In GSA, active gravitational mass **Ma**, passive gravitational mass **Mp** and inertial mass **Mi** are accepted as they are equal to each other.

Gravitiational constant which is propotional with the time is also updated according to the iteration number given by the equation 19.

$$\mathbf{G(t) = G\_0 e^{-\mathbf{a} \frac{\mathbf{t}}{T}}} \tag{19}$$

Where, **α** is a user specified constant, **T** is the total number of iterations, and **t** is the current iteration. This equation is similar to the SA acceptance probability equation given by the equation 6. According to the formula the gravity force is decreasing by the time. The graviational force and total forces are updated as follows.

$$\mathbf{F\_{i\dot{j}}^{\mathbf{d}}(t) = G(t)} \frac{\mathbf{M\_{pi}(t)M\_{a\dot{j}}(t)}}{R\_{i\dot{j}}(t) + \varepsilon} (\mathbf{x\_{\dot{j}}^{\mathbf{d}}(t) - x\_{\dot{i}}^{\mathbf{d}}(t)}) \tag{20}$$

$$\mathbf{F\_{i}^{d}}(\mathbf{t}) = \sum\_{\mathbf{j} \in \mathbf{Kbest}, \mathbf{j} \neq \mathbf{i}}^{\mathbf{N}} \mathbf{rand}\_{\mathbf{j}} \mathbf{F\_{i}^{d}}(\mathbf{t}) \tag{21}$$

**Fijd(t)** is the gravity force acting on mass *i* from mass *j* at time t, **G(t)** is the gravitational constant at time t, **Mpi** is the passive gravitational mass of object i, **Maj** is the active gravitational mass of object j, is a small constant, **Rij** is the Euclidean distance between two objects i and j . **Fi d(t)** is the force that acts on particle *i* in dimension *d*, **randj** is a random number between 0 and 1, **Kbest** is the set of first K objects with the best fitness value and biggest mass.The acceleration of an object *i* in *d*th dimension is as follows.

$$\mathbf{a\_{i}}^{\mathbf{d}}(\mathbf{t}) = \frac{F\_{\mathbf{i}}^{\mathbf{d}}(\mathbf{t})}{\mathbf{M\_{i}}(\mathbf{t})} \tag{22}$$

The new positions and velocities of the particles are calculated by the equation 23 and 24.

$$\mathbf{v\_1^d(t+1) = rand\_1 \times v\_1^d(t) + a\_1^d(t)}\tag{23}$$

$$\mathbf{x}\_{\mathbf{i}}^{\mathbf{d}}(\mathbf{t}+1) = \mathbf{x}\_{\mathbf{i}}^{\mathbf{d}}(\mathbf{t}) + \mathbf{v}\_{\mathbf{i}}^{\mathbf{d}}(\mathbf{t})\tag{24}$$

Heuristic Optimization Algorithms in Robotics 323

0001 *R P <sup>T</sup>*

Fig. 9. The nearly Puma robot's link parameters and coordinate systems replacement

**I αi-1 ai-1 di θ<sup>i</sup>** 1 0 0 h θ<sup>1</sup> 2 -π/2 0 d1 θ<sup>2</sup> 3 0 l1 0 θ<sup>3</sup> 4 0 l2 0 θ<sup>4</sup> 5 -π/2 l3 0 θ<sup>5</sup> 6 π/2 0 l4 θ<sup>6</sup>

(25)

Fig. 8. DH parameters between two coordinate systems

Table 1. The DH parameters of the robot

According to the law of motion, particles try to move towards the heavier objects. The heavier masses represent good solutions and they move slowly for the exploitation. The pseudo code of GSA is given in Figure 7.

Fig. 7. The pseudo codes of GSA
