AN INTRODUCTION OF KRILL HERD ALGORITHM FOR ENGINEERING OPTIMIZATION

. A new metaheuristic optimization algorithm, called Krill Herd (KH), has been recently proposed by Gandomi and Alavi (2012). In this study, KH is introduced for solving engineering optimization problems. For more verification, KH is applied to six design problems reported in the literature. Further, the performance of the KH algorithm is com- pared with that of various algorithms representative of the state-of-the-art in the area. The comparisons show that the results obtained by KH are better than the best solutions obtained by the existing methods.


Introduction
The engineering optimization problems are commonly non-linear. They have different design variables under complex constraints. These constraints can be considered as simple bounds or as non-linear relationships. The non-linearity of such optimization problems often leads to a multimodal response landscape (Yang 2010;Gandomi et al. 2013a). Consequently, only global optimization algorithms can be used to obtain optimal solutions. Metaheuristic algorithms can be defined as upper level general methodologies. They can be used as guiding strategies in designing underlying heuristics to handle engineering optimization problems (Gandomi et al. 2013a;Talbi 2009). The main characteristics of metaheuristics are: (1) intensification and (2) diversification (Yang 2009). Intensification searches around the current best solutions and selects the best candidates or solutions. Diversification guarantees that the algorithm can explore the search space more efficiently.
The main goals of developing modern metaheuristic methods are to solve problems faster, to solve large problems, and to obtain robust algorithms (Talbi 2009). The most typical types of metaheuristics are genetic algorithms (GA) and particle swarm optimization (PSO). The efficiency of metaheuristic algorithms is related to the fact that they imitate the best features in nature.
Krill Herd (KH) algorithm is a new metaheuristic search algorithm. This algorithm is based on simulating the herding behaviour of krill individuals using a Lagrangian model and crossover. This algorithm is developed by Gandomi and Alavi (2012) and the preliminary studies show that it is very promising and could outperform existing algorithms (Gandomi et al. 2013c(Gandomi et al. , 2013d. In this paper, the KH algorithm is further validated against various engineering optimization problems. The introduced search strategy is compared with other popular optimization algorithms. Finally, the unique features of KH are discussed and topics for further studies are proposed.

Lagrangian model of the krill herding
Predators remove individuals, reduce of the average krill density, and distance the krill swarm from the food location. Therefore, predation can be considered as the initialization of the optimization algorithm. The fitness of each individual in the natural system, is supposed to be the distances from the food centre and the highest density of the krill swarm (Gandomi, Alavi 2012). Hofmann et al. (2004) proposed three effective factors of individual krill position as: i. movement induced by other krill individuals; ii. foraging activity; and iii. random diffusion, which can be formulated during the time and for n dimensional space, using the following Lagrangian model: where N i , F i and D i are respectively the motions i, ii and iii.

Motion induced by other krill individuals
For the i th krill individual, the induced motion is formulated as: (2) and , ( where N max is the maximum induced speed, α i is the direction of motion, ω n is the inertia weight, N i old is the last induced motion, α i local is the local effect provided by the neighbors and α i target is the target direction effect provided by the best krill individual. According to Hofmann et al. (2004), N max considered to be equal to 0.01 (ms -1 ).
The effect of the neighboursin a krill movement individual can be formulated as (Gandomi, Alavi 2012): where K worst and K best are, respectively, the worst and the best fitness values of the krill individuals; K i represents the fitness value of the i th krill individual; K j is the fitness of j th neighbour and j ∈ {1, 2, …, NN} and X is the positions in the search domain. ε is suggested to be a small number (Gandomi, Alavi 2012). In this study, this parameter is equal to 10 -6 .
For choosing a neighbour krill for the i th krill individual, a sensing distance (d s,i ) is defined using: , where N is the number of the krill individuals and i ∈ {1, 2, …, N}. Based on this equation, if the distance of two krill individuals is less than sensing distance, they are neighbors (Gandomi, Alavi 2012). The effect of the best fitness krill into the i th individual krill is formulated as: , where, C best is an empirical effective coefficient as: where rand is a uniform random value between 0 and 1, I and I max are, respectively, the actual and maximum number of iterations (Gandomi, Alavi 2012).

Foraging motion
There are two main terms in the foraging motion, the food attraction and the previous krill experience, which can be formulated as follows (Gandomi, Alavi 2012): and (11) where V f is the foraging speed, ω f is the inertia weight, F i old is the last foraging motion, β i food is the food attraction and β i best is the effect of the best fitness of the i th krill during its history. Foraging speed is taken 0.02 (ms -1 ) based on (Price 1989).
In each iteration, the centre of food of can be defined like centre of mass as follows: (12) and, the food attraction for the can be formulated as: , where C food is empirically defined as: The effect of the best fitness of the i th krill individual during the history is defined as: , where K ibest is the best previously visited position by the i th krill individual.

Physical diffusion
The physical diffusion is a random process which can be formulated as follows (Gandomi, Alavi 2012): where D max is the maximum diffusion speed, and δ is the uniform random directional vector between -1 and 1. Based on the suggested values in Morin et al. (1988), D max ∈ [0.002, 0.010](ms -1 ).

Motion process of the KH algorithm
Using the three explained motions, the position vector of a krill individual during the time interval from t to t+Δt is formulated as (Gandomi, Alavi 2012): , where Δt can be obtained from: where NV is the number of variables, LB j and UB j are lower and upper bounds of the j th variables, respectively, and C t is a constant number which is considered as 0.5 in this study.

Crossover
As it is evaluated in the original paper, crossover is an effective process in the KH algorithm. By generating a uniformly distributed random vector values between 0 and 1, the m th component of X i , X i,m is manipulated as (Gandomi, Alavi 2012): where Cr i is crossover probability of the i th kril individual.

Constraint handling
In order to solve the problem simpler, nonlinear constraints in the penalty function approach might be collapsed with the cost function into a response functional. This results in transformation of the constrained optimization problem into an unconstrained optimization one.
The following example clarifies the issue. Assuming that there are some nonlinear equality constraints ф i and some inequality constraints ψ j , the response functional ∏ can be defined as follows ): where: 1 ≤ µ i and 0 ≤ v i . The coefficients of penalty terms should be large enough; their values may depend on the specific optimization problem. The contribution of any equality constraints function to the response functional ∏ is null but increases notably as soon as the constraint is violated. The same applies to inequality constraints when they become critical ). If integer/discrete design variables are involved in an optimization problem, the variable is rounded to the nearest integer/discrete value.

Implementation and numerical experiments
Engineering optimization problems are complex, sometimes even the optimal solutions of interest do not exist. In order to see how the KH algorithm performs, four standard engineering test problems are solved. It should be noted that because of the random nature of the KH algorithm, 50 trials with independent population initializations have been made to obtain a better conclusion of the performance. Figure 1 presents an example for designing a uniform column of tubular section to carry a compressive load P = 2500 kgf at minimum cost (Rao 1996). The column is made of a material with a yield stress (σ y ) of 500 kgf/cm 2 , a modulus of elasticity (E) of 0.85 × 10 6 kgf/cm 2 , and a density (ρ) equal to 0.0025 kgf/cm 3 . The length (L) of the column is 250 cm. The stress included in the column should be less than the buckling stress (constraint g 1 ) and the yield stress (constraint g 2 ). The mean diameter of the column is restricted between 2 and 14 cm (constraint g 3 and g 4 ), and columns with thickness outside the range 0.2-0.8 cm are not commercially available (constraint g 5 and g 6 ). The cost of the column includes material and construction costs (Hsu, Liu 2007). It is taken as the objective function. The optimization model of this problem is given as follows:

Case I. Tubular column design
Minimize: .
Subject to: where: 2 ≤ d ≤ 14 and 0.2 ≤ t ≤ 0.8.  Table 1 illustrates the statistical results for the best objective value by KH when 10,000 searches havebeen done in each run. The statistical values presented in this table clearly show the proposed algorithm is successful in this case. Table 2 compares the results obtained by KH with those reported in the literature (Rao 1996;Hsu, Liu 2007;Gandomi et al. 2013b;Rocha, Fernandes 2009). It can be observed from Table 2 that the best objective values by Rao (1996) and Hsu and Liu (2007) and are not feasible because the second constraint (g 2 ) is violated. The Result of KH algorithm is also better than results obtained in Gandomi et al. (2013a), Rocha and Fernandes (2009). Therefore, KH algorithm provides better results than other algorithms.

Case II. Three-bar truss design
This case considers a 3-bar planar truss structure shown in Figure 2. This problem was first presented by Nowcki (1974) and it is one of the benchmark structural engineering problems . The volume of a statically loaded 3-bar truss is to be minimized subject to stress (σ) constraints on each of the truss members. The objective is to evaluate the optimal cross sectional areas (A 1 , A 2 ). The mathematical formulation is given as below: Minimize: where H is shown in Figure 2 and it is equal to 100 cm.
This design problem is a nonlinear fractional programming problem. The statistical values of the best solution obtained by the KH algorithm are given in Table 3. From this table, the optimized costs corresponding to worst and best designs are very close to each other. The solution by KH is (A 1 , A 2 ) = (0.78867, 0.40902) with the objective value equal to 263.97156 after 5,000 function evaluation in each run. Table 4 presents the solutions obtained by KH and those reported by other methods reported in the literature. As it is seen, the best objective value reported by Tsai (2005) is not feasible because the first constraint (g 1 ) is violated. Hence, it can be concluded that the results obtained by KH are better than those of previous studies for this benchmark problem.

Case III. Speed reducer design
KH is applied to the design of a speed reducer which is a benchmark structural optimization problem   (Fig. 3), with the face width (b), module of teeth (m), number of teeth on pinion (z), length of shaft 1 between bearings (l 1 ), length of shaft 2 between bearings (l 2 ), diameter of shaft 1 (d 1 ), and diameter of shaft 2 (d 2 ). The objective is to minimize the total weight of the speed reducer. The constraints involve limitations on the bending stress of the gear teeth, surface stress, transverse deflections of shafts 1 and 2 due to transmitted force, and stresses in shafts 1 and 2.
The mathematical formulation can be summarized as minimizing the following function: where: 2.6 ≤ b ≤ 3.6, 0.7 ≤ m ≤ 0.8, 17 ≤ z ≤ 28, 7.3 ≤ l 1 ≤ 8.3, 7.8 ≤ l 2 ≤ 8.3, 2.9 ≤ d 1 ≤ 3.9, and 5.0 ≤ d 1 ≤ 5.5. The corresponding statistical values of the Best KH model are presented in Table 5. From this table, the ratio between the optimized costs corresponding to worst and best designs is 1.00 and it shows KH algorithm has successfully found the optimum design in all runs. Table 6 presents a comparison of the results obtained by KH and other methods. As it is seen, the KH results are better than those reported by Hsu and Liu (2007) Montes and Ocana (2008). Although the best objective values derived by Kuang et al. (1998), Li and Papalambros (1985), and Azarm and Li (1989) are better than those of KH, the reported values are not feasible. The result provided by Yang and Gandomi (2012) is best one and the results obtained in this study are the third best fitness value. It should be noted that some previous studies consider the simple bound of l 2 like l 1 so they are not considered in the comparison study. Figure 4 shows a schematic representation of the helical compression spring design problem ). The spring is subject to an axially guided constant compression load. It must be designed for mini-  Notes: a DE is differential evolution; b BA is bat algorithm; c SC is society and civilization; d EA is evolutionary algorithm; e NP is non-linear programming; f Hybrid electromagnetism like algorithm; g Violated set.   Table 7  .

Case IV. Helical compression spring design
The design shear stress caused by the compression load should be lower than the allowable maximum shear stress (S) of the material (g 1 ). The free length of the spring should be shorter than the maximum specified value L free (g 2 ). The wire diameter must not be less than the specified minimum diameter d min (g 3 ). The outerdiameter of the coil should be smaller than the specified maximum diameter D max (g 4 ). The inner coil diameter must be at least three times less than the wire diameter to avoid a lightly wound spring (g 5 ). The deflection under the given load δ must be less than the specified maximum deflection under preload δ pm (g 6 ). The combined deflection must be consistent with the coil free length L free (g 7 ).
The deflection from preload to maximum load must be greater than the specified working deflection δ w (g 8 ). The cost function and constrained of the problem to be minimized are the spring volume, expressed as: Minimize: .
Subjected to: where: The values assigned to constant terms involved in the spring design problem statement are listed in Table 8. The optimization results obtained by KH are presented in Table 9. The optimization process was completed within 30,000 function evaluations. The ratio between the optimized costs corresponding to worst and best designs is 1.14.

Conclusions
The new KH algorithm is utilized to solve engineering optimization problems. The results indicate that the KH algorithm is very efficient for solving engineering problems. More, it performs superior to different existing algorithms. It can be because of this fact that there are fewer parameters to be fine-tuned in KH than in other algorithms. As agents generally contribute to the moving of each other based on their fitnesses, therefore the violated agents cannotaffect a lot on the others. In addition, a neighbour agent has an attractive/repulsive effect on the movement of the agent, therefore the neighbours with the better fitness, violated or not violated, attract the agent and other neighboursrepulse it. As an instance, if the agent isslightly violating the constraints, the agents with more violations repulse it and the feasible agents or agents with less violation attract the agent. These effects can act as an effective local search for each krill individual and the results generally prove it. KH can be viewed as a system with multiple interacting Markov chains selected and biased towards global optimality. This powerful optimization method can be extended to study multi-objective optimization applications with various constraints, including NP-hard problems (Gandomi, Alavi 2012).