Applied Mathematics and Computation 204 (2008) 694–701
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A derivation of the number of minima of the Griewank function Huidae Cho a, Francisco Olivera a,*, Seth D. Guikema b a b
Department of Civil Engineering, Texas A&M University, 3136 TAMU, College Station, TX 77843-3136, United States Department of Geography and Environmental Engineering, Johns Hopkins University, Baltimore, MD, United States
a r t i c l e
i n f o
Keywords: Griewank function Local minima Optimization Multi-modal optimization
a b s t r a c t The Griewank function is commonly used to test the ability of different solution procedures to find local optima. It is important to know the exact number of minima of the function to support its use as a test function. However, to the best of our knowledge, no attempts have been made to analytically derive the number of minima. Because of the complex nature of the function surface, a numerical method is developed to restrict domain spaces to hyperrectangles satisfying certain conditions. Within these domain spaces, an analytical method to count the number of minima is derived and proposed as a recursive functional form. The numbers of minima for two search spaces are provided as a reference. Ó 2008 Elsevier Inc. All rights reserved.
1. Introduction The Griewank function [1] has been widely used to test the convergence of optimization algorithms [2–15] because its number of minima grows exponentially as its number of dimensions increases [7,14]. The function is defined as follows:
fn ð~ xÞ ¼
n n Y 1 X xi x2i cos pffi þ 1; 4000 i¼1 i i¼1
within ½600; 600n where n is the number of dimensions of the function. The global minimum is located at ~ 0 with a value of 0. The actual number of minima may not be important when global optimization is performed, but it needs to be known to test techniques that search for local optima. Most studies vaguely mention the number of minima of the Griewank function [7–9], and, to the best of our knowledge, no analytical derivation to determine it has been given in the literature. Knowing the number of minima is critical if the Griewank function serves as the basis for evaluating algorithms designed to find local minima as well as global ones (i.e., multi-modal optimization). In some cases [14], the number of solutions given is inconsistent with analytical results. For example, [14] compared the ability of NichePSO, nbest PSO, lbest PSO, sequential niching, and deterministic crowding based on the number of minima found through numerical searches. However, further work with another algorithm has found a different number of solutions than found by [14]. In order to address this issue and provide a consistent basis for comparing algorithms, this paper analytically derives the number of minima of the Griewank function. We develop an approach in three basic steps. First, we restrict the search space to a hyperrectangle. Second, we show that the hyperrectangle is the maximum possible hyperrectangle of the Griewank function within which local minima on the Griewank function correspond to tangent points on a simpler surface. Third, we develop an analytical approach for counting the number of the tangent points on the simpler surface. This approach yields an accurate count of the number of local minima of the Griewank function within the defined hyperrectangle.
* Corresponding author. E-mail addresses:
[email protected] (H. Cho),
[email protected] (F. Olivera),
[email protected] (S.D. Guikema). 0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.07.009
695
H. Cho et al. / Applied Mathematics and Computation 204 (2008) 694–701
Section 2 elaborates on the characteristics of the function surface and redefines the problem of counting the number of minima to make it analytically tractable. Because of the complex nature of the function surface, the domain space needs to be restricted to hyperrectangles found by the numerical method introduced in Section 2. Although the analytical method to determine the number of minima derived in Section 3 cannot be applied to arbitrary domain spaces, it should be noted that the method does not miss any minima within hyperrectangles satisfying certain conditions. As most optimization algorithms are tested within fixed hyperrectangles, it remains practical to use hyperrectangles as domain spaces for testing many optimization algorithms. 2. Redefinition of the problem The partial derivative of the Griewank function with respect to xi is
! x piffi n sin Y ~ ofn ðxÞ xi xj i ¼ cos pffi : þ pffi oxi 2000 j i j¼1;j–i It is difficult, if not impossible, to analytically solve this non-linear system involving n variables. Global and local minima have to satisfy the following conditions:
! n sin pxiffii Y x xj i 0 ~ fn;i ðxÞ ¼ cos pffi ¼ 0 for i ¼ 1; . . . ; n; þ pffi 2000 i j j¼1;j–i ! n 1 1 Y xj 00 ~ ðxÞ ¼ cos pffi > 0 for i ¼ 1; . . . ; n; ; þ fn;i 2000 i j¼1 j
ð1Þ ð2Þ
0 ~ 00 ~ where fn;i ðxÞ and fn;i ðxÞ are the first and second derivatives of fn ð~ xÞ, respectively. Note that i is an index for dimensions. Inequality (2) is required to ensure that maxima are not taken into account. By rearranging (2), we obtain Qn Qn xj xj i p ffi p ffi > 2000. Because the region of non-positive values of satisfying (1) and (2) (i.e., fn ð~ xÞ P j¼1 cos j¼1 cos j j P Pn 2 n 1 1 2 ~ j¼1 xj þ 1 at local minima) is outside of the region of its positive values (i.e., fn ðxÞ < 4000 j¼1 xj þ 1 at local minima), 4000 problem domains in this paper are restricted such that
n Y j¼1
xj cos pffi j
Since a value of follows:
!
i 2000
> 0:
ð3Þ
is small for low dimensions, not much portion of the function space is lost. Eq. (1) can be rewritten as
!#1 pffi " n xi xi i Y xj sin pffi ¼ cos pffi ; 2000 j¼1;j–i j i
ð4Þ
Q Q x x where nj¼1;j–i cos pj ffi –0 because nj¼1 cos pj ffi > 0. j j Pn 2 1 Because fn ð~ xÞ meets the surface 4000 x at the global minimum and near local minima, we will find the minima of fn ð~ xÞ j¼1 j Pn 2 1 xÞ on the simpler surface 4000 two sets by finding the tangent points of fn ð~ j¼1 xj and deriving the relationship between these Pn 2 1 of points. In the following, tangent points refer to the tangent points of the Griewank function on the surface 4000 j¼1 xj unless otherwise noted. Since we only want to know the number of minima, their exact coordinates are not of direct interest. In this paper, the number of minima is indirectly derived by counting the number of tangent points associated with them. Because the tangent point associated with the global minimum is the global minimum itself, this method also takes into account the global minimum. Therefore, problem domains have to be carefully defined so that there exists one minimum xi 0 ~ ðxÞ also tends to increase along the line 2000 and, eventually, no points satisfying for each tangent point. As i or xi increases, fn;i (1) are found, which makes global optimization easier [7]. Because there are high correlations between dimensions in highQn x dimensional problems, it is hard to determine whether or not there are local minima satisfying 0 < j¼1 cos pj ffi < 1 by j
0 ~ inspecting fn;i ðxÞ surfaces separately. It is necessary to know the maximum extent of each xi beyond which there are no local minima associated with tangent points as shown in Fig. 1. For n ¼ 1, it is trivial to check the maximum extent of x1 because 0 ðx1 Þ. For n P 2, a numerical analysis is required to estimate the corners of the hyperrectangle beyond all the points lie on f1;1 which there exist tangent points not associated with any local minima. While tangent points are evenly distributed at every pffi 2p i, local minima are not. If the boundary of a domain space is located between a tangent point and its corresponding local minimum, the number of pffitangent points is not the same as the number of local minima. For this reason, a problem domain U maximum value p offfi ki , ki;max , is defined such that the largest local minimum is defined as U ¼ ð0; 2p iki Þ where ki 2 N. The p ffi th associated with a tangent point is located in ð2p iðki;max 1Þ; 2p iki;max Þ. Using the periodicity of the sine curve, the ki local k k minimum, ~ xki ¼ ðx1i ; . . . ; xni Þ, is obtained by solving the following shifted version of (4):
696
H. Cho et al. / Applied Mathematics and Computation 204 (2008) 694–701
Fig. 1. Out-most region of one dimension of the Griewank function beyond which there exist no more minima. Note that only the tangent point on the lefthand side of the gray region is associated with a local minimum. Problem domains should be smaller than the hyperrectangle defined by the gray region for the method presented in this paper to be valid.
x0i ki pffi i
sin
! ¼
" !#1 pffi k n xj i x0i ki i þ 2piðki 1Þ Y cos pffi ; 2000 j j¼1;j–i
ð5Þ
pffi k where x0i 1 ¼ x1i and x0i ki ¼ xi i 2p iðki 1Þ. Q x For one-dimensional problems, (5) is further simplified by setting n ¼ 1 and nj¼1;j–i cos pj ffi ¼ 1. If the both sides of (5) j pffi Q x meet at x0i ¼ 32 p i, there are no local minima at this point because the value of nj¼1 cos pj ffi does not satisfy (3). By solving j
3 2
pi þ 2piða 1Þ 2000
"
n Y j¼1;j–i
!#1
xj cos pffi j
¼ 1;
where a 2 R, we obtain
ki;max
$ % ! n Y 1000 xj 1 ¼ bac ¼ cos pffi þ ; 4 pi j j¼1;j–i
where bc is the maximum integer less than or equal pffito a given number (i.e., the flooring function). However, since the Griewank function is defined within ½600; 600n , 2p iki;max must be less than or equal to xmax ¼ 600. Therefore, ki;max is
%! ! $ n Y xmax 1000 xj 1 pffi ; cos pffi þ 4 pi j 2p i j¼1;j–i
ki;max ¼ min
ð6Þ
pffi and, given a one-dimensional domain space ð0; 2p iki Þ where 1 6 ki 6 ki;max , ki is the number of local minima. In problems of more than one dimension, because the position of a local minimum in one axis is highly correlated with those in the other axes, it is not trivial to analytically solve (5) for all dimensions. The values of cosðpxiffiiÞ and ki;max for i ¼ 1; . . . ; n can be numerically estimated with the pseudo code presented in Fig. 2. The subroutine defined in Fig. 3 is used to solve (5) for each dimension at a time. x0i k found in this way may not be the correct one because the correlation between Qn x 0k dimensions is not taken into account when solving (5). An estimated value of xi is used to evaluate j¼1;j–i cos pj ffi , which is j 0k iteratively plugged into (5) to estimate the next value of xi . Once ki;max is estimated, a problem domain needs to be defined. Define a problem domain by U¼ ð0; xi;max Þ, where pffi pffi Qn xj 0 < xi;max 6 2p iki;max , such that xi;max does not have to be 2p iki where 1 6 ki 6 ki;max . When j¼1;j–i cos pffi is greater than j pffi pffi pffi 0, a local minimum is found in ð2p iki 12 p i; 2 p iki Þ because cosðpxiffiiÞ is greater than 0 satisfying (3), and (4) can hold true pffi pffi pffi pffi Q x only in this range. Likewise, when nj¼1;j–i cos pj ffi is less than 0, a local minimum is found in ð2p iki 32 p i; 2p iki p iÞ. j
Thus, xi;max needspto ffi avoid pffithese ranges because, otherwise, it is possible to find local minima not associated with tangent points at xi ¼ 2p iki p i, which means that the analytical method introduced in this paper cannot be applied. Therefore, the allowable range of xi;max is either
H. Cho et al. / Applied Mathematics and Computation 204 (2008) 694–701
Fig. 2. Pseudo code to estimate cos
697
x piffi for i ¼ 1; . . . ; n. i
Fig. 3. Pseudo code for the getxi subroutine.
pffi
pffi pffi 3 pffi 2p iki 2p i; 2p iki p i 2 or
pffi
pffi pffi 1 pffi 2p iki p i; 2p iki p i : 2 The above conditions for xi;max can be interpreted as
pffi 0 < xi;max 6 2p iki;max and
( xi;max 2 X ¼
xi jxi is a multiple of
ð7Þ $ % ) xi i _ p pffi is an even integer : i 2
p pffi 2
ð8Þ
698
H. Cho et al. / Applied Mathematics and Computation 204 (2008) 694–701
In case xi;max does not satisfy (8) because integer values for xi;max are preferred, we need to make sure that there are no local minima in
$
! % xi;max p pffi pffi i ; x i;max ; p i 2 2
ð9Þ
pffi where 0 < xi;max 6 2p iki;max and xi;max R X. This test can be done indirectly by checking whether or not the distance in the ith axis between xi;max and the closest tangent point whose coordinate is greater than xi;max is greater than the possibly largest distance between them. The closest tangent point whose coordinate is greater than xi is
&
’ xi p pffi i: ti ðxi Þ ¼ p pffi 2 i 2 Likewise, the largest distance between a local minimum and its corresponding tangent point is obtained by calculating k k ti ðxi i;max Þ xi i;max because t i ðxki Þ is the tangent point associated with xki , and the distance between them also increases as xi ink k creases. If ti ðxi;max Þ xi;max is greater than t i ðxi i;max Þ xi i;max , there must be one local minimum in ðxi;max ; ti ðxi;max ÞÞ along the ith axis, which means that there are no local minima in the range defined by (9). When xi;max satisfies all the requirements described above, a domain space can be extended to U ¼ ½xi;max ; xi;max 8i 2 f1; . . . ; ng because the negative domain space ðxi;max ; 0Þ is symmetrical to ð0; xi;max Þ, and the analytical method derived in the following section takes into account both regions implicitly. 3. Derivation of the number of minima In the previous section, the problem was redefined so that the number of tangent points is the as the number of same Qn xj minima. The cosine function is defined function j¼1 cos pffi is also restricted to in ½1; 1 and, thus, the range of the j i h Qn P Pn 2 xj 1 1 ~ p ffi has a value in ½0; 2 and fn ðxÞ in 4000 nj¼1 x2j ; 4000 ½1; 1. Consequently, 1 j¼1 cos xj þ 2 . Therefore, tangent j¼1 j Qn Pn 2 xj 1 p ffi is 1. xÞ lie on the surface 4000 points of fn ð~ j¼1 xj when j¼1 cos j pffi xi The absolute value of cos pffii is 1 when xi is a multiple of p i. The times j cosðpxiffiiÞj equals to 1 depends on the range of xi or pffi ½xi;min ; xi;max . The number of p ik, where k 2 Z, within this range is calculated as
Ni ¼
xi;max x pffi i;min pffi þ 1: p i p i
The number of xi ’s satisfying cos
Nþi ¼
x piffi ¼ 1 is i
xi;max x pffi i;min pffi þ 1 2p i 2p i
and the number of xi ’s satisfying cos
Ni ¼ Ni Nþi ¼
x piffi ¼ 1 is i
xi;max 1 xi;min 1 pffi þ pffi : 2p i 2 2p i 2
Now, the number of maxima and minima can be expressed as M n ¼ Counting the number of n-tuples in the set
Qn
j¼1 N j .
( ! ) n Y x1 xn xj n An ¼ cos pffiffiffi ; . . . ; cos pffiffiffi cos pffi ¼ 1 2 ½1; 1 j n 1 j j¼1
is a combinatorial problem where combinations take place without repetitions. Any element, cos pxiffii , of n-tuples belonging Q x to the set An must have a value of 1 or 1 because, otherwise, the absolute value of nj¼1 cos pj ffi cannot be 1. Because j Qn xj p ffi should be 1, an even number of elements in an n-tuple have a value of 1, and the other elements have a value cos j¼1 j
of 1. Therefore, the number of n-tuples in the set An can be expressed as
b2nc X n j¼0
bn2c X j¼0
n! ; ðn 2jÞ!ð2jÞ!
n is the binomial coefficient. Encode n-tuples in An as ða1 ; a2 ; . . . ; an Þ where ai is either 1 or 1. If 1 and 1 are 2j substituted with þ and symbols, respectively, n-tuples in An can be represented as (+, +, , +), (, , +, , +), (, +, , n n +, , +) (i.e., one n-tuple of and two examples of , respectively), and so on. Note that there are an even number 0 2
where
2j
¼
699
H. Cho et al. / Applied Mathematics and Computation 204 (2008) 694–701
of symbols, and the others are all þ’s. The numbers of xi values satisfying ai ¼ þ and ai ¼ are N þ i and N i , respectively. Counting all the possible ~ x vectors generating n-tuples in the set An can be done recursively in terms of n. Let Sn be the numfor n ¼ 1. For n ¼ 2, there are S1 minima when a2 is ber of minima for n-dimensional problems. The simplest form is S1 ¼ N þ 1 Qn1 Qn1 Qn Qn1 aj ¼ þ also satisfy fixed to þ because S1 number of x1 ’s satisfying j¼1 j¼1 aj an ¼ j¼1 aj ¼ þ. If a2 is fixed to , j¼1 aj must be , and the number of a1 satisfying this condition is M 1 S1 (i.e., the number of maxima for n ¼ 1). Therefore, S2 ¼ S1 N þ 2 þ ðM 1 S1 Þ N 2 . Generalizing this recursive form, the following equations are obtained:
S1 ¼ Nþ1
if n ¼ 1;
ð10Þ
Sn ¼ Sn1 Nþn þ ðMn1 Sn1 Þ Nn
if n > 1
ð11Þ
for ½xi;min ; xi;max 8i 2 f1; . . . ; ng. Now, (10) and (11) can be expanded as follows:
k lx m 1;min þ 1 if n ¼ 1; 2p 2p xn;max x min pffiffiffi n; p ffiffiffi þ 1 Sn ¼ Sn1 2p n 2p n " # $ % & ’ ! n1 Y xj;max xj;min pffi pffi þ 1 Sn1 þ p j p j j¼1 xn;max 1 xn;min 1 pffiffiffi þ pffiffiffi if n > 1 2p n 2 2p n 2 S1 ¼
jx
1;max
ð12Þ
ð13Þ
for ½xi;min ; xi;max 8i 2 f1; . . . ; ng. 4. Results and discussion k
100
Fig. 4 and Table 1 present the maximum estimated number of local minima, ki;max , and the largest local minimum, xi i;max , th on the i axis. They define hyperrectangles within which (12) and (13) can be applied. Outside these regions, the analytical method presented in this paper cannot be used to count the number of minima. Fig. 4 shows ki; max for different dimensions. For n P 43, the numerical algorithm in Fig. 2 experienced difficulties in finding ki;max , and no plots were drawn. This result might be caused by reducing the search space by 2p in all directions. However, because the number of minima within only a k small fractioni of hyperrectangles defined by xi i;max is so high even for n ¼ 3 (e.g., 1215 minima in ½28; 283 , a subspace of h k k xi i;max ; xi i;max 8i 2 f1; 2; 3g), it would be practically enough to define domain spaces for up to n ¼ 40. Table 1 shows ki;max k and xi i;max estimated for up to three-dimensional problems. Note that ki;max for the same i varies with n because of the correlation between dimensions. When defining a domain space by U¼ ½xi;max pffi ; xi;max 8i 2 f1; . . . ; ng, we need to make sure k ;max k ; max . This condition satisfies (7) because ti xi i ¼ 2p iki;max for all the cases in Table 1. Also, xi;max 0 < xi;max 6 ti xi i has to satisfy (8) or (9). As a set of examples, domain spaces U ¼ ½xmax ; xmax n were evaluated for 1 6 n 6 3 where xmax 2 f14; 28g. Note that, for the sake of simplicity, domain spaces were chosen such that all xi;max ¼ xmax . For xmax ¼ 14, (8) holds true when i = 1 or 2. The
n=1
50
n=2 n=3 n=4
20
10
n=5 n=10
n=30
2
5
n=40
1
k i.max
n=20
0
10
20
30
i Fig. 4. ki;max versus i for different problem dimensions.
40
700
H. Cho et al. / Applied Mathematics and Computation 204 (2008) 694–701
Table 1 Maximum estimated number of local minima and the largest local minimum in each dimension n
i
ki;max
xi i;max
k ti xi i;max
k k ti xi i;max xi i;max
1 2
1 1 2 1 2 3
95 94 66 88 62 51
596.60 590.28 585.82 552.45 550.04 553.78
596.90 590.62 586.46 552.92 550.92 555.02
0.30 0.34 0.64 0.47 0.88 1.24
3
k
ki;max k is the largest local minimum on the ith axis; t i xi i;max is ki;max is the maximum estimated number of local minima on the ith axis within xi 2 ð0; 600Þ; xi ki;max ki;max xi its corresponding tangent point; and t i xi is the largest distance in the ith axis between them.
Table 2 Numbers of minima for ½14; 14n and ½28; 28n n
½14; 14n
½28; 28n
1 2 3
5 31 157
9 111 1215
closest tangent point whose coordinate is greater than x3;max ¼ 14 is t 3ð14Þ, and the distance in the 3rd axis between x3;max k k and t3 ð14Þ is t 3 ð14Þ 14 ¼ 2:32. This distance is greater than t3 x33;max x33; max ¼ 1:24 for n ¼ 3 as shown in Table 1. This means that the local minimum associated with t3 ðx3;max Þ exists in ðx3;max ; t3 ðx3;max ÞÞ, not in the range defined by (9) for x3;max ¼ 14. For xmax ¼ 28, (8) holds true when i ¼ 2 or 3. A visual inspection of the x1 axis and a numerical analysis show that there are no local minima in the range defined by (9) for x1;max ¼ 28. Because xmax 2 f14; 28g satisfies the boundary conditions specified by (8) and (9), we can safely use (12) and (13) to calculate the number of minima of the Griewank function. Table 2 shows the numbers of minima for the two search spaces for up to three dimensions. 5. Conclusions It is difficult to analytically solve the derivative of the Griewank function and directly count the number of minima because of the complex nature of the function surface. The problem of counting the number of minima was redefined as counting the number of tangent points lying on a parabolic surface. A numerical method was developed to find hyperrectangles within which this approach can be applied, and the number of minima of the function was analytically derived within these domain spaces based on a recursive functional form. The maximum extents of hyperrectangles for up to three dimensions were estimated, and the numbers of minima for two search spaces were provided as a reference. The numerical and analytical methods introduced in this paper can be used to determine the exact number of minima within the domain space defined by a hyperrectangle satisfying certain conditions. The number of minima derived in this paper can serve as a sound basis for evaluating multi-modal optimization algorithms. Acknowledgement This work was supported in part by the Korean government (the Ministry of Science and Technology) through the Korea Science and Engineering Foundation Grant (No. M06-2003-000-10064-0). References [1] A.O. Griewank, Generalized descent for global optimization, Journal of Optimization Theory and Applications 34 (1) (1981) 11–39. [2] J. Kennedy, Stereotyping: Improving particle swarm performance with cluster analysis, in: Proceedings of the Congress on Evolutionary Computation, IEEE Service Center, Piscataway, New Jersey, 2000, pp. 1507–1512. [3] T. Krink, J.S. Vesterstrøm, J. Riget, Particle swarm optimisation with spatial particle extension, in: Proceedings of the Congress on Evolutionary Computation, vol. 2, Honolulu, Hawaii, 2002, pp. 1474–1479. [4] J. Riget, J.S. Vesterstrøm, A diversity-guided particle swarm optimizer–the ARPSO, Technical Report 2002-02, Department of Computer Science, Aarhus Universitet, Bgn. 540, Ny Munkegade DK-8000 Aarhus C, Denmark, 2002. [5] X.-F. Xie, W.-J. Zhang, Z.-L. Yang, A dissipative particle swarm optimization, in: Proceedings of the Congress on Evolutionary Computation, Honolulu, Hawaii, 2002, pp. 1456–1461. [6] R. Brits, A.P. Engelbrecht, F. van den Bergh, Scalability of Niche PSO, in: Proceedings of the 2003 IEEE Swarm Intelligence Symposium, 2003, pp. 228– 234. . [7] M. Locatelli, A note on the Griewank test function, Journal of Global Optimization 25 (2003) 169–174. [8] N. Khemka, C. Jacob, Exploratory toolkit for evolutionary and swarm-based optimization, in: Proceedings of the 6th International Mathematica Symposium, Banff, Alberta, Canada, 2004. [9] S. He, Q.H. Wu, J.Y. Wen, J.R. Saunders, R.C. Paton, A particle swarm optimizer with passive congregation, BioSystems 78 (2004) 135–147. [10] A. Acan, A. Gunay, Enhanced particle swarm optimization through external memory support, in: Proceedings of the Congress on Evolutionary Computation 2 (2005) 1875–1882.
H. Cho et al. / Applied Mathematics and Computation 204 (2008) 694–701
701
[11] C.K. Monson, K.D. Seppi, Bayesian optimization models for particle swarms, in: Proceedings of the Genetic and Evolutionary Computation Conference, ACM Press, NY, New York, 2005, pp. 193–200. [12] M. Meissner, M. Schmuker, G. Schneider, Optimized particle swarm optimization (OPSO) and its application to artificial neural network training, BMC Bioinformatics 7 (2006) 125–136. [13] C.K. Monson, K.D. Seppi, Adaptive diversity in PSO, in: Proceedings of the Genetic and Evolutionary Computation Conference, ACM Press, NY, New York, 2006, pp. 59–66. [14] R. Brits, A.P. Engelbrecht, F. van den Bergh, Locating multiple optima using particle swarm optimization, Applied Mathematics and Computation 189 (2007) 1859–1883. [15] X. Wang, X.Z. Gao, S.J. Ovaska, A hybrid optimization algorithm based on ant colony and immune principles, International Journal of Computer Science & Applications 4 (3) (2007) 30–44.