RelHead2(i,1) = RelHead2(i,1) + pi; end
while RelHead2(i,1) > pi
RelHead2(i,1) = RelHead2(i,1) - pi; end end
noise = (rand(m,1)*eta) -eta/2;
newHeading2(:,1) = heading2(:,1) + RelHead2(:,1)./sum(neighbours2,2) + noise(:,1);
newPos2(:,1) = pos2(:,1) + cos(newHeading2(:,1))*speed; newPos2(:,2) = pos2(:,2) + sin(newHeading2(:,1))*speed; for k=1:m
while newPos2(k,1) < 0
newPos2(k,1) = newPos2(k,1) + fieldsize; end
while newPos2(k,1) >= fieldsize
newPos2(k,1) = newPos2(k,1) - fieldsize; end
while newPos2(k,2) < 0
newPos2(k,2) = newPos2(k,2) + fieldsize; end
while newPos2(k,2) >= fieldsize
newPos2(k,2) = newPos2(k,2) - fieldsize; end end
heading2 = newHeading2; pos2 = newPos2;
meanHeading2(time,1) = mean(heading2,1); deviationMean2(time,1) =
mean(abs(meanHeading2(time,1)-heading2(:,1)));
fieldsize = sqrt(n/density); pos1= fieldsize*rand(n,2); heading1 = 2*pi*rand(n,1); relposX1 = zeros(n); relposY1 = zeros(n); neighbours1 = zeros(n); RelHead1 = zeros(n,1); B = zeros(n,1);
newHeading1= zeros(n,1); newPos1 = zeros(n,2); meanHeading1= zeros(TIME,1);
deviationMean1 = zeros(n,TIME); %?-à′3ìDòóDfor run=1:RUNS for
21
time=1:TIME￡???????é???μ?2?·?
for i=1:n for j=1:n for s=1:m
relposX1(i,j) = abs(pos1(i,1) - pos1(j,1)); relposY1(i,j) = abs(pos1(i,2) - pos1(j,2)); relposX3(j,s) = abs(pos2(s,1) - pos1(j,1)); relposY3(j,s) = abs(pos2(s,2) - pos1(j,2));
if(sqrt(relposX1(i,j)^2 + relposY1(i,j)^2) <= range | sqrt(relposX3(j,s)^2 + relposY3(j,s)^2)<=2*range) neighbours1(i,j)=1; end end end end for i=1:n for j=1:n
if(neighbours1(i,j)==1)
relheading1(i,j) = heading1(j,1)-heading1(i,1); end end end
RelHead1 = sum(relheading1,2); for i=1:n
while RelHead1(i,1) < -pi
RelHead1(i,1) = RelHead1(i,1) + pi; end
while RelHead1(i,1) > pi
RelHead1(i,1) = RelHead1(i,1) - pi; end end
noise = (rand(n,1)*eta) -eta/2;
newHeading1(:,1) = heading1(:,1) + RelHead1(:,1)./sum(neighbours1,2) + noise(:,1);
newPos1(:,1) = pos1(:,1) + cos(newHeading1(:,1))*speed; newPos1(:,2) = pos1(:,2) + sin(newHeading1(:,1))*speed; for k=1:n
while newPos1(k,1) < 0
newPos1(k,1) = newPos1(k,1) + fieldsize; end
while newPos1(k,1) >= fieldsize
22
newPos1(k,1) = newPos1(k,1) - fieldsize; end
while newPos1(k,2) < 0
newPos1(k,2) = newPos1(k,2) + fieldsize; end
while newPos1(k,2) >= fieldsize
newPos1(k,2) = newPos1(k,2) - fieldsize; end end
heading1 = newHeading1; pos1 = newPos1;
meanHeading1(time,1) = mean(heading1,1); deviationMean1(time,1) =
mean(abs(meanHeading1(time,1)-heading1(:,1)));
scatter(pos2(:,1),pos2(:,2), 'og');hold on scatter(pos1(:,1),pos1(:,2), 'xr');hold off
axis([0 fieldsize 0 fieldsize]); xlabel(time ); M(time) = getframe; end
runDev(:,run) = deviationMean1(:,1); runDev(:,run) = deviationMean2(:,1); end
%plot(runDev)
movie2avi(M,'boidtest.avi', 'quality',100);
附录4
function [pso F] = pso_2D()
pop_size = 30; % pop_size 种群大小
gbest = zeros(1,part_size+1); % gbest 当前搜索到的最小的值 max_gen = 60; % max_gen 最大迭代次数 region=zeros(part_size,2); % 设定搜索空间范围 region=[-3,3;-3,3];
rand('state',sum(100*clock)); % 重置随机数发生器状态 arr_present = ini_pos(pop_size,part_size); % present 当前位置,随机初始化,rand()的范围为0~1
v=ini_v(pop_size,part_size); % 初始化当前速度 pbest = zeros(pop_size,part_size+1);
% pbest 粒子以前搜索到的最优值，最后一列包括这些值的适应度
23
part_size = 2; % part_size 粒子大小,
w_max = 0.7298; % w_max 权系数最大值 w_min = 0.4;
v_max = 2;
% **最大速度,为粒子的范围宽度 c1 =2; % 学习因子1 c2 = 2; % 学习因子2
best_record = zeros(1,max_gen); % best_record记录最好的粒子的适应度。 arr_present(:,end)=ini_fit(arr_present,pop_size,part_size);
pbest = arr_present; %初始化各个粒子最优值
[best_value best_index] = min(arr_present(:,end)); %初始化全局最优，即适应度为全局最小的值，根据需要也可以选取为最大值 gbest = arr_present(best_index,:); x=[-3:0.01:3]; y=[-3:0.01:3];
z=@(x,y) 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... - 1/3*exp(-(x+1).^2 - y.^2); for i=1:max_gen grid on; ezmesh(z)
,hold
on,
grid
on,plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'), text(-0.8,-1.2,6.5,'食物');title('鸟群觅食');hold off; drawnow
F(i)=getframe; pause(0.01);
w = w_max-(w_max-w_min)*i/max_gen;
reset = 0; % reset = 1时设置为粒子群过分收敛时将其打散，如果＝1则不打散
if reset==1 bit = 1;
for k=1:part_size
bit = bit&(range(arr_present(:,k))<0.1); end
if bit==1 % bit=1时对粒子位置及速度进行随机重置
arr_present = ini_pos(pop_size,part_size); % present 当前位置,随机初始化
v = ini_v(pop_size,part_size); % 速度初始化 for k=1:pop_size
arr_present(k,end) = fitness(arr_present(k,1:part_size)); end
warning('粒子过分集中！重新初始化……'); % 给出信息 display(i); end
24
end
for j=1:pop_size v(j,:) ...
=
w.*v(j,:)+c1.*rand.*(pbest(j,1:part_size)-arr_present(j,1:part_size))
+c2.*rand.*(gbest(1:part_size)-arr_present(j,1:part_size));
% 粒子速度更新 (a)
c = find(abs(v)>1.5); %**最大速度设置，粒子的范围宽度
v(c) = sign(v(c))*1.5; %如果速度大于3.14则，速度为3.14 arr_present(j,1:part_size)
arr_present(j,1:part_size)+v(j,1:part_size);
% 粒子位置更新 (b)
arr_present(j,end) = fitness(arr_present(j,1:part_size)); if
(arr_present(j,end)>pbest(j,end))&(Region_in(arr_present(j,:),region))
% 根据条件更新pbest,如果是最小的值为小于号,相反则为大于号 pbest(j,:) = arr_present(j,:); end end
[best best_index] = max(arr_present(:,end)); % 如果是最小的值为min,相反则为max
if best>gbest(end)&(Region_in(arr_present(best_index,:),region)) % 如果当前最好的结果比以前的好，则更新最优值gbest,如果是最小的值为小于号,相反则为大于号
gbest = arr_present(best_index,:); end
best_record(i) = gbest(end); display(i); end
pso = gbest; display(gbest);
function fit = fitness(present) %**需要求极值的函数,本例即peaks函数 - 10*(present(1)/5
=
fit=3*(1-present(1)).^2.*exp(-(present(1).^2) - (present(2)+1).^2) ...
-
present(1).^3
-
present(2).^5).*exp(-present(1).^2-present(2).^2) ... - 1/3*exp(-(present(1)+1).^2 - present(2).^2);
function ini_present=ini_pos(pop_size,part_size)
25
ini_present = 3*rand(pop_size,part_size+1);
%初始化当前粒子位置,使其随机的分布在工作空间 %** 6即为自变量范围
function ini_velocity=ini_v(pop_size,part_size) ini_velocity =3/2*(rand(pop_size,part_size)); %初始化当前粒子速度,使其随机的分布在速度范围内
function flag=Region_in(pos_present,region) [m n]=size(pos_present); flag=1; for j=1:n-1
flag=flag&(pos_present(1:j)>=region(j,1))&(pos_present(1:j)<=region(j,2)); end
function arr_fitness=ini_fit(pos_present,pop_size,part_size) for k=1:pop_size
arr_fitness(k,1) = fitness(pos_present(k,1:part_size)); %计算原始种群的适应度 end
26