cuckooSearch.jl 6.18 KB
 Fredrik Bagge Carlson committed Sep 07, 2015 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ``````using Devectorize """ `cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e5, timeout = Inf)`\n `n` = Number of nests (or different solutions) `pa=0.25` Discovery rate of alien eggs/solutions Change this if you want to get better results Based on implementation by @inproceedings{yang2009cuckoo, title={Cuckoo search via L{\'e}vy flights}, author={Yang, Xin-She and Deb, Suash}, booktitle={Nature \& Biologically Inspired Computing, 2009. NaBIC 2009. World Congress on}, pages={210--214}, year={2009}, organization={IEEE} } http://www.mathworks.com/matlabcentral/fileexchange/29809-cuckoo-search--cs--algorithm """ `````` 18 ``````function cuckoo_search(f,X0, Lb,Ub;n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e3, timeout = Inf) `````` Fredrik Bagge Carlson committed Sep 09, 2015 19 `````` X00 = deepcopy(X0) `````` Fredrik Bagge Carlson committed Sep 07, 2015 20 21 22 23 24 25 26 27 28 29 30 31 32 `````` nd=size(X0,1); X0t = X0' Lb = Lb' Ub = Ub' if !all(isfinite(Lb)) Lb=X0t-0.99999*abs(X0t); end if !all(isfinite(Ub)) Ub=X0t+0.99999*abs(X0t); end # Random initial solutions nest = zeros(n,nd) `````` Fredrik Bagge Carlson committed Sep 09, 2015 33 `````` nest[1,:] = X0t `````` Fredrik Bagge Carlson committed Sep 07, 2015 34 35 `````` for i=2:n nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb)); `````` Fredrik Bagge Carlson committed Sep 09, 2015 36 37 `````` DEBUG && @assert !any(nest[i,:] .> Ub) DEBUG && @assert !any(nest[i,:] .< Lb) `````` Fredrik Bagge Carlson committed Sep 07, 2015 38 39 `````` end # Get the current best `````` Fredrik Bagge Carlson committed Sep 09, 2015 40 `````` `````` Fredrik Bagge Carlson committed Sep 07, 2015 41 42 `````` fitness=10^20*ones(n,1); fmin,bestnest,nest,fitness=get_best_nest(f,nest,nest,fitness); `````` Fredrik Bagge Carlson committed Sep 09, 2015 43 44 `````` DEBUG && println("f(X0) = \$(f(X00)), f(bestnest) = \$(fmin)") DEBUG && @assert X00 == X0 `````` Fredrik Bagge Carlson committed Sep 07, 2015 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 `````` N_iter=0; t0 = time() ## Starting iterations while fmin>Tol && N_iter < max_iter # Generate new solutions (but keep the current best) new_nest=get_cuckoos(nest,bestnest,Lb,Ub); fnew,best,nest,fitness=get_best_nest(f,nest,new_nest,fitness); # Update the counter N_iter += n; if fnew timeout display("Cuckoo search: timeout \$(timeout)s reached (\$(time()-t0)s)") break end # Discovery and randomization new_nest=empty_nests(nest,Lb,Ub,pa) ; # Evaluate this set of solutions fnew,best,nest,fitness=get_best_nest(f,nest,new_nest,fitness); # Update the counter again N_iter += n; # Find the best objective so far if fnew timeout display("Cuckoo search: timeout \$(timeout)s reached (\$(time()-t0)s)") break end end ## End of iterations ## Post-optimization processing ## Display all the nests println("Total number of iterations=",N_iter); `````` Fredrik Bagge Carlson committed Sep 09, 2015 81 `````` println("f(bestnest) = \$(fmin)") `````` Fredrik Bagge Carlson committed Sep 07, 2015 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 `````` squeeze(bestnest',2),fmin end ## --------------- All subfunctions are list below ------------------ ## Get cuckoos by ramdom walk function get_cuckoos(nest,best,Lb,Ub) # Levy flights n=size(nest,1); # Levy exponent and coefficient # For details, see equation (2.21), Page 16 (chapter 2) of the book # X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010). beta=3/2; sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta); for j=1:n s=nest[j,:]; # This is a simple way of implementing Levy flights # For standard random walks, use step=1; ## Levy flights by Mantegna’s algorithm u=randn(size(s))*sigma; v=randn(size(s)); betai = 1/beta @devec step=u./abs(v).^betai; # In the next equation, the difference factor (s-best) means that # when the solution is the best solution, it remains unchanged. stepsize=0.01*step.*(s-best); # Here the factor 0.01 comes from the fact that L/100 should the typical # step size of walks/flights where L is the typical lenghtscale; # otherwise, Levy flights may become too aggresive/efficient, # which makes new solutions (even) jump out side of the design domain # (and thus wasting evaluations). # Now the actual random walks or flights s=s+stepsize.*randn(size(s)); # Apply simple bounds/limits nest[j,:]=simplebounds(s,Lb,Ub); `````` Fredrik Bagge Carlson committed Sep 09, 2015 116 117 `````` DEBUG && @assert !any(nest[j,:] .> Ub) DEBUG && @assert !any(nest[j,:] .< Lb) `````` Fredrik Bagge Carlson committed Sep 07, 2015 118 `````` end `````` Fredrik Bagge Carlson committed Sep 09, 2015 119 `````` return nest `````` Fredrik Bagge Carlson committed Sep 07, 2015 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 ``````end ## Find the current best nest function get_best_nest(f,nest,newnest,fitness) # Evaluating all new solutions for j=1:size(nest,1) fnew=f(squeeze(newnest[j,:]',2)); if fnew<=fitness[j] fitness[j]=fnew; nest[j,:]=newnest[j,:]; end end # Find the current best (fmin,K) = findmin(fitness) ; best=nest[K,:]; `````` Fredrik Bagge Carlson committed Sep 09, 2015 135 `````` return fmin,best,nest,fitness `````` Fredrik Bagge Carlson committed Sep 07, 2015 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 `````` end ## Replace some nests by constructing new solutions/nests function empty_nests(nest,Lb,Ub,pa) # A fraction of worse nests are discovered with a probability pa n=size(nest,1); # Discovered or not -- a status vector K=rand(size(nest)).>pa; # In the real world, if a cuckoo’s egg is very similar to a host’s eggs, then # this cuckoo’s egg is less likely to be discovered, thus the fitness should # be related to the difference in solutions. Therefore, it is a good idea # to do a random walk in a biased way with some random step sizes. ## New solution by biased/selective random walks stepsize=rand()*(nest[randperm(n),:]-nest[randperm(n),:]); new_nest=nest+stepsize.*K; `````` Fredrik Bagge Carlson committed Sep 09, 2015 152 153 154 155 `````` for j = 1:size(nest,1) new_nest[j,:]=simplebounds(new_nest[j,:],Lb,Ub); end return new_nest `````` Fredrik Bagge Carlson committed Sep 07, 2015 156 157 158 159 160 ``````end # Application of simple constraints function simplebounds(s,Lb,Ub) # Apply the lower bound `````` Fredrik Bagge Carlson committed Sep 09, 2015 161 162 `````` I = s.Ub; s[J] = Ub[J]; return s `````` Fredrik Bagge Carlson committed Sep 07, 2015 167 168 169 170 171 172 173 174 175 176 177 178 ``````end # # ## You can replace the following by your own functions # # # A d-dimensional objective function # # function fobj(u) # # ## d-dimensional sphere function sum_j=1^d (u_j-1)^2. # # # with a minimum at (1,1, ...., 1); # # z=sum((u-1).^2); # # end # # dims = 10 # # cuckoo_search(fobj,zeros(dims),-10*ones(dims),10*ones(dims))``````