动物集群运动行为模型系列之五-毕业论文.doc

动物集群运动行为模型

动物集群运动模型
摘要

本文主要模拟了鱼群的集群运动、鱼群躲避捕食者追捕的运动情况以及鸟群觅食运动的模拟，以此研究动物个体间的信息传递机制，同时也是对群体智能的初步探索。
针对问题一，需要我们给出对鱼群集群运动的模型，并编写程序将运动模拟出来，对此我们建立了Boid模型，根据模型给出的准则以及算法，我们通过matlab编程，在忽略阻力等因素下分别模拟出在平面以及空间鱼群的运动，并得出密度必须大于一定值时，鱼群才能最终达到同步。鱼群的整个集群运动从刚开始的随机产生的各个个体的不均匀无规则分布到逐渐的聚拢成群再到最后的一致方向的前进。
针对问题二，我们在问题一的模型的基础上给出了鱼群躲避捕食者的模型，制定了鱼个体的适度逃离区域和加速逃离区域，分析捕食者与鱼个体的关系，给出进一步的模型，通过编写程序得到模拟的结果，得到了对鱼群躲避捕食者的运动的合理的动态模拟，并且给出了模型的改进方向。
针对问题三，我们更加倾向于研究鸟群的觅食行为，因此我们将问题改成鸟群的觅食模拟，将鸟群的觅食行为转化为求最优解的问题，这正好与问题中提到了有一部分个体掌握食物源位置信息相对应。针对问题，我们建立了粒子群优化模型，通过PSO算法，通过鸟群寻找食物的最短路径的最优解的问题的分析，我们利用优化算法来模拟了鸟群在山间的觅食行为，得到了鸟群可以绕过我们设定的障碍物（山峰）到达食物点。

关键字：动物集群运动 Boid模型 PSO算法 鸟群觅食

一、问题重述
在动物界，大量集结成群进行移动或者觅食的例子并不少见，这种现象在食草动物、鸟、鱼和昆虫中都存在。这些动物群在运动过程中具有很明显的特征：群中的个体聚集性很强，运动方向、速度具有一致性。通过数学模型来模拟动物群的集群运动行为以及探索动物群中的信息传递机制一直是仿生学领域的一项重要内容。
通过观察附件中给出的图片和视频资料，或者在网上搜索相关资料观察，思考动物集群运动的机理，建立数学模型刻画动物集群运动、躲避威胁等行为，例如，可以考虑以下问题的分析建模：
1. 建立数学模型模拟动物的集群运动。
2. 建立数学模型刻画鱼群躲避黑鳍礁鲨鱼的运动行为。
3. 假定动物群中有一部分个体是信息丰富者(如掌握食物源位置信息，掌握迁徙路线信息)，请建模分析它们对于群运动行为的影响，解释群运动方向决策如何达成。
建议与说明：
1.在上述问题的讨论中，如果能适时分析动物群中的信息传递机制无疑是更好的。
2.如果对问题2和问题3之外的其他集群运动行为更感兴趣，也可将这两个问题替换为你所感兴趣的问题来讨论。
3.建模过程中的数据资料可以在网上查询或者自行合理设定。若果感到在三维空间讨论问题太复杂，可以先在二维空间讨论，再推广至三维空间。
4.最好能对你所做的机理分析模型给出计算机仿真方法以便于实际情况对比评价。

二、模型假设
1.忽略障碍、阻力以及其它无关次要因素对于集群运动的影响
2.问题一鱼群中每个个体运动的速度都是恒定一样的
3.鱼群集群运动的模拟中不考虑障碍物的存在
4.忽略其它种群对本文所研究种群的影响
5.不考虑集群中个体的体积，都按粒子处理

三、符号说明
鱼群的总数
集群中每个个体的位置矢量
集群中每个个体的速度矢量

集群中每个个体运动的速度
排斥区域的半径
一致区域的半径
吸引区域的半径
惯性权重
粒子数
空间维数
最小速度
最大速度
粒子的位置矢量
捕食者的位置矢量
个体在时刻的预期方向
、 学习因子（加速因子）
、 均匀分布在（0,1）之间的随机数
在第次迭代时粒子的位置表示
在第次迭代时粒子的速度表示
个体极值
全局极值

四、问题分析
本问题是一个动物集群运动的模型问题，动物的集群运动包括很多，其中有觅食、追尾、躲避捕食者等等运动，问题一需要我们考虑动物集群的运动模型，也就是鱼群的游弋、鸟群的飞翔等行为，是不需要考虑觅食、追尾等行为活动的，我们通过建立Boid模型进行鱼群集群运动的模拟。问题二需要我们给出鱼群躲避捕食者的运动模拟，要解决问题，那就需要我们在问题一模型的基础上给出鱼的逃逸模式，然后对逃逸运动进行模拟。问题三是需要我们模拟集群运动中存在领导者时的集群运动的模拟，可以运用和问题一一样的思路，但是我们对于鸟群的觅食行为更感兴趣，所以我们转而对鸟群觅食进行建模，我们选择PSO算法，通过模型求最优解的过程对鸟群觅食行为进行模拟，从而建立起了比较合适的模型。针对这些问题，我们主要的工作是首先建立合适的模型，通过我们建立的模型根据模型的算法，我们可以编写程序得到对集群运动的模拟。

五、模型的建立与求解

5.1鱼群集群运动的模拟

5.1.1模型的建立

我们根据问题一的要求，通过查阅资料得出了Boid模型可以解决类似的问题，这里我们就选用Boid模型作为此文的模型。
Boid群模型包括三个简单的指导规则，它们描述了一个单一的“Boid”如何基于位置和邻近个体的速度进行分离、内聚和排序活动的。模型的三个规则（Reynolds聚结规则）如下：
（1） 群中心定位：试图与邻近的群个体保持接近；
（2） 避免障碍：避免与邻近的群个体发生冲突；
（3） 速度匹配：试图与邻近的群个体速度匹配；
设系统有个个体组成，它们的位置和速度矢量分别为、，每个个体在三维空间中按照恒定的速度运动，为个体在时刻的预期方向。在每一步，每个个体可以感知到三个不重叠的区域中其他个体的位置和角度，这些信息用于计算，这三个区域分别为：排斥区域，一致区域，吸引区域，也称为Three-circle模型，其模型的三个区域如图（图1）所示：

图1 Boid模型的三个区域示意图，zor为排斥区域，zoo为一致区域，zoa为吸引区域；为视野盲区，取值范围为
个体的运动规则为：首先，每个个体尽量与排斥区域（以该个体为中心，以为半径的球）中的其他个体保持最小距离，并记其中的个体数为，则个体的预期方向按照下面的方式调整

其中。其次，如果，则个体的预期方向受“一致区域”（以个体位中心，处于和之间的球形区域，除去该个体后面的角度为的盲区）及“吸引区域”（以个体位中心，处于和之间的球形区域，除去该个体后面的角度为的盲区）中的个体的影响，记相应区域中的邻居个数分别为、，则可定义、如下：

如果，则；同样，如果，则；如果两者都不为0，则定义如果经过上面的运算后所得到的，或者在三个区域中都没有个体，则。设旋转速率为，即每一时步个体所能转过的最大角度为，如果与之间的角度差小于，则，否则，个体向期望的方向旋转角度，这样就得到了个体下一步的运动方向。

5.1.2模型的求解

根据模型，我们利用matlab编程，得到二维空间（平面）和三维空间（立体）的鱼群运动的模拟。

5.1.2.1平面鱼群运动的模拟

程序见附录1。我们取了个体为100，个体的旋转速率为0.5，得到了运动的模拟动画。
下面我们截取了整个动态模拟过程中的3幅图。

图2 第63时步时动态截图

图3 第95时步时动态截图

图4第104时步时动态截图

5.1.2.2空间鱼群运动的模拟

程序见附录2。我们取了个体为100，个体的旋转速率为0.5，得到了运动的模拟动画。
下面我们截取了整个动态模拟过程中的3幅图。

图5个体初始时刻的状态截图

图6个体开始聚集收拢的截图

图7个体一致向一个方向运动截图
从上面的各个图中，我们可以看出，鱼群的整个集群运动从刚开始的随机产生的各个个体的不均匀无规则分布到逐渐的聚拢成群再到最后的一致方向的前进。整个模拟过程是比较合理的，比较符合实际，由于论文上无法显示动态的图，所以我们采取了截图的方式。
从我们的模拟中，我们还发现，取不同的密度，得到的模拟结果是不一样的，在高密度的情景下，我们得出鱼群在经过有限的运动时间后，会最终达到同步，即运动方向达到一致。而且只有在密度大于一定范围时系统才能最终达到一致同步，随着个体密度的减小，我们可以推测存在一个临界密度，只有在大于临界密度的时候系统才能最终达到同步。密度越大，越容易达到同步。密度约小，鱼群只能形成多个小规模的鱼群，无法形成整个鱼群运动一致。
取不同的旋转速率时，所得到的模拟动画也不一样。通过输入不同的旋转速率，我们得出结论：旋转速率越大，鱼群总是在原位置附近徘徊，无法达到一致同步的趋势；旋转速率较小时，能够达到最终的一致同步，即使不能达到一致同步，也能以小规模的集群同步运动。

5.2鱼群躲避捕食者的运动的模拟

5.2.1模型的建立

在问题一的模型的基础上，为考虑有外来捕食者的情况下群集的应急机制，设计了如下的粒子个体应急措施。如下图所示

图8 捕食者模型示意图
粒子个体的应急区域分为两个层次：适度逃离区和加速逃离区。图阴影部分即为加速逃离区，内圆与外圆间的环形部分是适度逃离区。图下半部分显示了捕食者分别位于两个区域时，粒子个体的受力大小和方向。为简化模型，规定加速逃离区同原模型的排斥区重合。当有捕食者进入到粒子个体的加速逃离区时，“逃跑”就是个体的最高策略，个体会沿捕食者与个体所在的直线告诉逃离捕食者，速度的大小同捕食者和个体间的距离存在一定的函数关系，此时个体将不考虑吸引、排斥异己方向同步等的作用。当有捕食者位于粒子个体的适度逃离区时，个体会受到捕食者与个体所在直线远离捕食者方向的作用，同时还将受到吸引、排斥以及方向同步等的作用。因此，个体受到的实际作用是两者作用的叠加。还应指出，当捕食者位于粒子个体的适度逃离区时，个体会加大方向同步区域的大小，期望尽快与其他粒子的方向保持同步。粒子的逃离速度公式如下：

(1) (2)
图9捕食者与粒子距离同逃离速度（1）与同步区域（2）的函数关系
其中是加速逃离区半径的大小，与排斥区域半径大小一样，即；是适度逃离区的半径大小，与吸引区域半径一样，即；是正常情况下同步区域半径大小，是此半径的上界（不超过）。
当捕食者位于粒子的加速逃离区时，粒子的实际速度就是其逃离速度，当捕食者位于粒子的适度逃离区时，粒子的实际速度是其逃离速度与在吸引、排斥以及方向同步等的作用下的速度的合成。

5.2.2模型的求解

通过matlab编程，我们得出了对鱼群逃逸的运动的模拟（程序见附录3）。
截图如下：

图10鱼群躲避捕食者的运动模拟截图

图11鱼群躲避捕食者的运动模拟截图

图12鱼群躲避捕食者的运动模拟截图（绿色圆圈为捕食者，红色叉为鱼）

结论：我们通过对鱼群躲避捕食者的模型建立，然后通过matlab编程得到了对鱼群运动的模拟，总的可以看出鱼群的逃逸运动。

5.3对鸟群的觅食运动的模拟

5.3.1模型的分析与建立

通过对问题的分析，我们建立了PSO模型对鸟群的觅食行为进行模拟。

5.3.1.1 PSO模型的简介

粒子群优化算法是基于群体的演化算法，其思想来源于人工生命和演化计算理论。Reynolds对鸟群飞行的研究发现，鸟仅仅是追踪它有限数量的邻居，但最终的整体结果是整个鸟群好像在一个中心的控制之下，即复杂的全局行为是由简单规则的相互作用引起的。PSO即源于对鸟群捕食行为的研究，一群鸟在随机的搜寻食物，如果区域只有一块食物，那么找到食物的最简单有效的策略就是搜寻目前离食物最近的鸟的周围区域。人们通常是以他们自己及他人的经验来作为决策的依据，这就构成了PSO的一个基本概念。

5.3.1.2算法的原理

PSO算法将群体中的每个个体视为多维搜索空间中的一个没有质量和体积的粒子，这些粒子在搜索空间中以一定的速度飞行并根据粒子本身的飞行经验以及同伴的飞行经验对自己的飞行速度进行动态调整，即每个粒子通过统计迭代过程中自身的最优解和群体的最优值来不断修正自己的前进方向和速度大小，从而形成群体寻优的正反馈机制。PSO算法就是这样依据每个粒子对环境的适应度将个体逐步移到较优的区域，并最终搜索、寻找到问题的最优解。PSO算法具有鲜明的生物社会背景：认知过程和社会行为，即在寻求一致的认知过程中，个体往往记住它们的信念，同时考虑其它同伴的信念，当个体察觉同伴的信念较好时，将进行适应性调整。

5.3.1.3模型的建立

在PSO算法中，用粒子的位置表示待优化问题的解，每个粒子性能的优劣程度取决于待优化问题目标函数确定的适应值，每个粒子由一个速度矢量决定其飞行方向和速率大小。设在一个维的目标搜索空间中，有个粒子组成一个群体，其中，在第次迭代时粒子组成一个群体，其中，在第次迭代时粒子的位置表示为，相应的飞行速度表示为。开始执行PSO算法时，首先随机初始化个粒子的位置和速度，然后通过迭代寻找最优解，在每一次迭代中，粒子通过跟踪两个极值来更新自己的速度和位置：一个极值是粒子本身迄今为止搜索到的最优解，称为个体极值，表示为；另一个极值是整个粒子群到目前为止找到的最优解，称为全局极值，表示为。在第次迭代计算时，粒子根据下列规则来更新自己的速度和位置：

取大值可使算法具有较强的全局搜索能力，取小值则算法倾向于局部搜索。一般的做法是将初始取0.9，并使其迭代次数的增加而线性递减至0.4，这样下去可以先侧重于全局搜索，使搜索空间快速收敛于某一区域，然后采用局部精细搜索以获得高精度的解；一般取2。另外粒子的每一维的速度取值范围为 。如果当前粒子的加速度导致它在某一维的速度超过该维上的最大速度，则该维的速度被限制值为最大速度。位置的取值范围为。
5.3.1.4 PSO算法流程

(1) 随机初始化粒子群体的位置和速度；通常是在允许的范围内随机产生的，每个粒子的 pbest 坐标设置为其当前位置，且计算出其相应的个体极值(即个体的适应度值)，而全局极值即全局的适应度值)就是个体极值中最好的，记录该最好值的粒子序号，并将gbest 设置为该最好粒子的当前位置；
(2) 计算每个粒子的适应值；
(3) 对每个粒子，将其适应值与个体极值进行比较，如果较优，则更新当前的个体极值；
(4) 对每个粒子，将其适应值与全局极值进行比较，如果较优，则更新当前的全局极值；
(5) 根据更新位置速度的公式，更新每个粒子的位置和飞行速度；
(6) 如未达到预先设定的停止准则(通常设置为最大迭代次数)，则返回步骤(2)，达到则停止计算；

5.3.2模型的求解

种群大小取30，粒子大小2，，均取2。
根据上述算法，我们通过matlab编程，得出了对鸟群觅食的模拟（程序见附录4）。截图如下：

图13鸟群绕过山峰觅食

图14鸟群飞向食物点

图15鸟群找到食物
结论：通过上述的截图我们可以看出，我们作出的模拟十分合理，鸟群能够绕过山峰找到食物，我们通过PSO模型，通过鸟群寻找食物的最短路径的最优解算法得到对鸟群觅食过程的模拟。
六、模型评价
对于问题一，我们采用的Boid模型是有缺陷的，它给出的规则都是局部规则，每个个体仅仅根据它周围附近球域内的个体的行为调整自己的行为，只能做到局部一致性，局限性较强，并且模型是在虚拟的没有障碍物存在的空间内进行模拟的，而实际中不可能没有障碍物，基于以上的缺点因素，我们提出了Vicsek模型，这个模型能克服上诉缺点。
对于问题二，我们在模型一的基础上建立了模型，模型得出的结果也是很合适的，不过模型对于鱼群的加速逃离区和适度逃离区的选择不是很合理，模型的改进方向可以是：划定一个必须逃离区，在必须逃离区外危险系数与捕食者和个体的距离成正比，在必须逃离区之内就只能全力逃跑。鉴于时间关系，我们没有再深入下去。
对于问题三，我们采用PSO模型，通过对鸟群觅食的最短路径最优化问题的分析，利用PSO模型作出对鸟群觅食的模拟。
参考文献
[1] 姜启源，数学模型[M].北京:高等教育出版社.1987年4月第一版；
[2] 胡守信，李柏年.基于MATLAB的数学实验[M].北京:科学出版社.2004年6月；
[3] 江铭炎，人工鱼群算法及其应用. 科学出版社. :2012年1月1日第一版；

附录
附录1
function boid1(n,eta)
%generates boids in a field
n=input('input n ');eta=input('input eta');
range =0.2;
speed =0.03;
TIME = 300;
density =10;
RUNS = 1;

fieldsize = sqrt(n/density);
pos = fieldsize*rand(n,2);
heading = 2*pi*rand(n,1);
relposX = zeros(n);
relposY = zeros(n);
neighbours = zeros(n);
RelHead = zeros(n,1);
B = zeros(n,1);
newHeading = zeros(n,1);
newPos = zeros(n,2);
meanHeading = zeros(TIME,1);
deviationMean = zeros(n,TIME);
for run=1:RUNS
for time=1:TIME
for i=1:n
for j=1:n
relposX(i,j) = abs(pos(i,1) - pos(j,1));
relposY(i,j) = abs(pos(i,2) - pos(j,2));
if(sqrt(relposX(i,j)^2 + relposY(i,j)^2) <= range)
neighbours(i,j)=1;
end
end
end
for i=1:n
for j=1:n
if(neighbours(i,j)==1)
relheading(i,j) = heading(j,1)-heading(i,1);
end
end
end
RelHead = sum(relheading,2);
for i=1:n
whileRelHead(i,1) < -pi
RelHead(i,1) = RelHead(i,1) + pi;
end
whileRelHead(i,1) > pi
RelHead(i,1) = RelHead(i,1) - pi;
end
end

noise = (rand(n,1)*eta) -eta/2;
newHeading(:,1) = heading(:,1) + RelHead(:,1)./sum(neighbours,2) + noise(:,1);
newPos(:,1) = pos(:,1) + cos(newHeading(:,1))*speed;
newPos(:,2) = pos(:,2) + sin(newHeading(:,1))*speed;
for k=1:n
whilenewPos(k,1) < 0
newPos(k,1) = newPos(k,1) + fieldsize;
end
whilenewPos(k,1) >= fieldsize
newPos(k,1) = newPos(k,1) - fieldsize;
end
whilenewPos(k,2) < 0
newPos(k,2) = newPos(k,2) + fieldsize;
end
whilenewPos(k,2) >= fieldsize
newPos(k,2) = newPos(k,2) - fieldsize;
end
end
heading = newHeading;
pos = newPos;
meanHeading(time,1) = mean(heading,1);
deviationMean(time,1) = mean(abs(meanHeading(time,1)-heading(:,1)));
scatter(pos(:,1),pos(:,2), 'xr');
axis([0 fieldsize 0 fieldsize]); xlabel(time );
M(time) = getframe;
end
runDev(:,run) = deviationMean(:,1);
end
%plot(runDev)
movie2avi(M,'boidtest.avi', 'quality',100);

附录2
function boid1(n,eta)
%generates boids in a field
n=input('input n');
eta=input('input eta');
range =0.2;
speed =0.05;
TIME = 300;
density =10;
RUNS 1;
fieldsize = sqrt(n/density);
pos = fieldsize*rand(n,3);
heading = 2*pi*rand(n,1);
relposX = zeros(n);
relposY = zeros(n);
relposZ = zeros(n);
neighbours = zeros(n);
RelHead = zeros(n,1);
B = zeros(n,1);
newHeading = zeros(n,1);
newPos = zeros(n,3);
meanHeading = zeros(TIME,1);
deviationMean = zeros(n,TIME);
for run=1:RUNS
for time=1:TIME
for i=1:n
for j=1:n
relposX(i,j) = abs(pos(i,1) - pos(j,1));
relposY(i,j) = abs(pos(i,2) - pos(j,2));
relposZ(i,j) = abs(pos(i,3) - pos(j,3));
if(sqrt(relposX(i,j)^2 + relposY(i,j)^2+relposZ(i,j)^2) <= range)
neighbours(i,j)=1;
end

end
end
for i=1:n
for j=1:n
if(neighbours(i,j)==1)
relheading(i,j) = heading(j,1)-heading(i,1);
end
end
end
RelHead = sum(relheading,2);
for i=1:n
whileRelHead(i,1) < -pi
RelHead(i,1) = RelHead(i,1) + pi;
end
whileRelHead(i,1) > pi
RelHead(i,1) = RelHead(i,1) - pi;
end
end
noise = (rand(n,1)*eta) -eta/2;
newHeading(:,1) = heading(:,1) + RelHead(:,1)./sum(neighbours,2) + noise(:,1);
newPos(:,1) = pos(:,1) + cos(newHeading(:,1))*speed;
newPos(:,2) = pos(:,2) + sin(newHeading(:,1))*speed;
newPos(:,3) = pos(:,3) + tan(newHeading(:,1))*speed;

for k=1:n
whilenewPos(k,1) < 0
newPos(k,1) = newPos(k,1) + fieldsize;
end
whilenewPos(k,1) >= fieldsize
newPos(k,1) = newPos(k,1) - fieldsize;
end
whilenewPos(k,2) < 0
newPos(k,2) = newPos(k,2) + fieldsize;
end
whilenewPos(k,2) >= fieldsize
newPos(k,2) = newPos(k,2) - fieldsize;
end
whilenewPos(k,3) < 0
newPos(k,3) = newPos(k,3) + fieldsize;
end
whilenewPos(k,3) >= fieldsize
newPos(k,3) = newPos(k,3) - fieldsize;
end
end
heading = newHeading;
pos = newPos;
meanHeading(time,1) = mean(heading,1);
deviationMean(time,1) = mean(abs(meanHeading(time,1)-heading(:,1)));
scatter3(pos(:,1),pos(:,2),pos(:,3) ,'*r');
axis([0 fieldsize 0 fieldsize 0 fieldsize]); xlabel(time );ylabel('Y');zlabel('Z');
M(time) = getframe;
end
runDev(:,run) = deviationMean(:,1);
end
%plot(runDev)
movie3avi(M,'boidtest.avi', 'quality',100);

附录3
function boid1(n,eta)
%generates boids in a field
n=input('input n ');eta=input('input eta');
range =0.2;
speed =0.03;
TIME = 300;
density =10;
RUNS = 1;
m=3;

fieldsize = sqrt(n/density);
pos2=fieldsize*rand(m,2);
heading2=2*pi*rand(m,2);
elposX2 = zeros(m);
relposY2 = zeros(m);
neighbours2 = zeros(m);
RelHead2 = zeros(m,1);
B = zeros(m,1);
newHeading2= zeros(m,1);
newPos2 = zeros(m,2);
meanHeading2= zeros(TIME,1);
deviationMean2 = zeros(m,TIME);
for run=1:RUNS

for time=1:TIME
for i=1:m
for j=1:m
relposX2(i,j) = abs(pos2(i,1) - pos2(j,1));
relposY2(i,j) = abs(pos2(i,2) - pos2(j,2));
if(sqrt(relposX2(i,j)^2 + relposY2(i,j)^2) <= range)
neighbours2(i,j)=1;
end
end
end
for i=1:m
for j=1:m
if(neighbours2(i,j)==1)
relheading2(i,j) = heading2(j,1)-heading2(i,1);
end
end
end
RelHead2 = sum(relheading2,2);
for i=1:m
while RelHead2(i,1) < -pi
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); %Ô­À´³ÌÐòÓÐfor run=1:RUNS for time=1:TIME£¬ÏÖ½«ÆäÉÏÃæµÄ²¿·Ö

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
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) = deviationMean