2018年认证杯SPSSPRO杯数学建模A题(第一阶段)海豚与沙丁鱼全过程文档及程序

2018年认证杯SPSSPRO杯数学建模

探究海豚猎捕时沙丁鱼群的躲避运动模型

A题 海豚与沙丁鱼

原题再现：

  沙丁鱼以聚成大群的方式来对抗海豚的捕食。由于水下光线很暗，所以在距离较远时，海豚只能使用回声定位方法来判断鱼群的整体位置，难以分辨每个个体。鱼群的行动是有协调性的，在没有外部威胁或障碍物时，鱼群常常会聚成接近球形的形态。而当海豚接触甚至冲进鱼群，鱼群则会进行协同的躲避，所以不易在大鱼群中追踪一个目标。沙丁鱼的这种群体行为降低了其被海豚捕食的概率。
  第一阶段问题： 请你建立合理的数学模型来描述沙丁鱼群在遇到一条海豚捕食时的运动规律。

整体求解过程概述(摘要)

  沙丁鱼为细长的银色小型鱼，游泳迅速，没有外部威胁时鱼群常聚成接近球形的形态，而当捕食者进攻时会进行协同的躲避，从而降低其被捕食的概率。本文基于沙丁鱼群遭遇一条海豚捕食时的情景，研究沙丁鱼群对应的运动规律。考虑到海豚和沙丁鱼群的同向、相向运动，以及二者之间的距离，分析海豚进击模型。当距离较远时，海豚通过回声定位来判断鱼群的整体位置，此时沙丁鱼群并未意识到威胁的存在，以接近球形的形态稳定行进。考虑个体视野范围以及速率差异，改进 Boid 模型，建立无序聚集运动模型；并划分感知区域，降维构建有序环绕运动模型。
  当海豚进攻沙丁鱼群时，鱼群进行协同躲避。考虑二者的相向、同向运动，分析四种外部空间结构（中空包围状、沙漏状、密集球状、驱赶聚集状）的动态变化。构造个体沙丁鱼安全场，划分为紧急逃离区域和适度逃离区域，分析各区域内影响个体沙丁鱼运动速度矢量的多种因素，采用主层次分析法提取出不同空间结构下不同区域内的主因子，建立四种空间结构下鱼群的躲避模型。考虑沙丁鱼空间密度等因素，引入随机森林算法对海豚捕食模型进行训练，动态优化四种躲避模型下沙丁鱼被捕食的概率结果，分析沙丁鱼群最优逃逸方式及运动规律。
  采集大自然中海豚捕食沙丁鱼群的场景，提取图片中的数据信息进行模型的验证，结果表明：沙丁鱼群遭遇一条海豚捕食时会优选呈中空包围状和密集球状进行逃逸。

问题分析：

  出于生存、避险、觅食等原因，鱼类的行为方式往往表现为群体行为，沙丁鱼往往以聚成大群的方式来对抗海豚的捕食。鱼群的行动是有协调性的，在没有外部威胁或障碍物时，鱼群常常会聚成接近球形的形态。而当海豚接触甚至冲进鱼群，一方面当近处的鱼快速逃避时，会给猎手海豚一个严重的心理错觉，以为远处的鱼还没发觉，于是扑向另外的对象。然而，鱼本身的侧线反馈机制会使远处的鱼逃避得更快，从而导致海豚找不到具体的目标，一无所获。
  此外，鱼类集群运动机理的研究涉及团体动力学、集群行为学等领域，围绕该方向的研究主要分为两大流派：一类是利用理论推导和应用已有的结论，试图去理解鱼类集群行为的内在原因，另一类是从鱼类的实际集群运动和个体鱼的特征出发，通过搭建模型去探索鱼类实际集群运动中存在的规律。由上述分析可以将海豚的捕食过程分为发现目标后的匀速前进阶段、变加速阶段、猛击捕食阶段。基于海豚的捕食过程，可以沙丁鱼群的存在状态分为稳定状态和躲避状态。由于鱼群的形成是由个体逐渐汇集而成，所以可以将沙丁鱼的稳定状态分为无序聚集状态和有序环绕状态。沙丁鱼在逃避过程中，沙丁鱼群群体的形状、空间分布、结构等是处于不断变化的过程中，鱼群在遇到外部攻击时，鱼群的外部空间结构会出现动态变化，会呈现驱赶聚集状、中空包围状、沙漏状和密集球状。
  建立数学模型来描述沙丁鱼群在遇到一条海豚捕食时的运动规律的实质就在于深入分析沙丁鱼群稳定集群模型、海豚的进击捕食模型以及沙丁鱼群多种的逃避模型，进而实现动态分析、动态寻优的过程。本文基于分析已有资料，采用从个体到局部、从局部到整体的思想，以 Boid 模型为基础，通过部分改良提出沙丁鱼群的无序聚集运动模型；采用降维的方法将三维模型降为二维模型，提出沙丁鱼稳定环绕状态模型以分析沙丁鱼的有序环绕运动。结合海豚的进击捕食模型，基于主成分分析方法简化沙丁鱼的运动方程，建立模型模拟驱赶聚集状、中空包围状、沙漏状和密集球状四种逃避轨迹，分析沙丁鱼在四种逃避状态下的被捕食概率，探索有利于沙丁鱼群体进化的最优逃避方式。

模型假设：

  1. 假设本文中的个体沙丁鱼的一般游动是随机性的；
  2. 假设每个沙丁鱼个体能够感知其在群体中的位置；
  3. 该沙丁鱼群不受环境变化影响（地震、海啸等自然灾害）；
  4. 假设沙丁鱼群中的个体是存在部分差异的，主要表现为尺寸上的大、小两种；
  5. 本文中模型所涉及的鱼群数量在短期内是不发生变化的（被海豚捕食情况除外），即不考虑个体沙丁鱼的突发死亡等非一般状况。

论文缩略图：

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：(代码和文档not free)

clear;
clc;
n=400;
x(1,:)=10*rand(1,n);
y(1,:)=10*rand(1,n);
drc(1,:)=(rand(1,n)-0.5)*2*pi;
d=[];
for k=1:80
 for i=1:n
 for j=1:n
 d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
 end
 end
 for i=1:n
 d(i,i)=inf;
 end
 for i=1:n
 a=0;b=1;c=0;e=0;h=0;l=0;
 for j=1:n
 if d(i,j)<5
 a=a+1;
 c=c+x(j);e=e+y(j);
 m=atan(y(j)/x(j));
 l=l+m;
 end
 if d(i,j)<1
 b=b+1;
 g=atan((y(j)-y(i))/(x(j)-x(i)));
 h=h+g;
 end
 end
 D2=atan(e/a-y(i))/(c/a-x(i));
 D3=1/a;
 D4=h/b;
 drc(k+1,i)=0.4*drc(k,i)+0.3*D2+0.2*D3+0.1*D4;
 v=0.5;
 end
 for i=1:n
 x(k+1,i)=x(k,i)+v*cos(drc(k+1,i))*1;
y(k+1,i)=y(k,i)+v*sin(drc(k+1,i))*1;
 if x(k+1,i)>30
 x(k+1,i)=x(k+1,i)-30;
 elseif x(k+1,i)<0
 x(k+1,i)=x(k+1,i)+30;
 end
 if y(k+1,i)>30
 y(k+1,i)=y(k+1,i)-30;
 elseif y(k+1,i)<0
 y(k+1,i)=y(k+1,i)+30;
 end
 end
end
for i=1:k
 %pause(0,1)
 plot(x(i,:),y(i,:),'.')
 axis([0 30 0 30])
 getframe
end

clear;
clc;
t=1;n=500;
x=[];y=[];xs=[];
x(1,:)=30*rand(1,n);y(1,:)=30*rand(1,n);
draction(1,:)=(rand(1,n)-0.5)*2*pi;
d=[];a=0;
xs(1,:)=30*rand(1,1);ys(1,:)=30*rand(1,1);
dractions(1,1)=(rand(1,1)-0.5)*2*pi;
for k=1:200
 dractions(k+1,1)=dractions(k,1);
 vs(k+1,1)=0.4;
 xs(k+1,1)=xs(k,1)+vs(k+1,1)*cos(dractions(k+1,1))*1;
 ys(k+1,1)=ys(k,1)+vs(k+1,1)*sin(dractions(k+1,1))*1;
 if xs(k+1,1)>30
 xs(k+1,1)=xs(k+1,1)-30;
 elseif xs(k+1,1)<0
 xs(k+1,1)=xs(k+1,1)+30;
 end
 if ys(k+1,1)>30
 ys(k+1,1)=ys(k+1,1)-30;
 elseif ys(k+1,1)<0
 ys(k+1,1)=ys(k+1,1)+30;
 end
for i=1:n
 ds(i)=sqrt((x(k,i)-xs(k,1))^2+(y(k,i)-ys(k,1))^2);
 for j=1:n
 d(i,j)=sqrt((x(k,i)-x(k,j))^2+(y(k,i)-y(k,j))^2);
 end
 end
 for i=1:n
 d(i,i)=inf;
 end
 for i=1:n
 A=0;B=0.1;
 for j=1:n
 if d(i,j)<0.3
 A=A-draction(k,j)/d(i,j);
 B=B-1/d(i,j);
 if A~=0
 draction(k+1,i)=(1-t)*draction(k,i)-t*j/B;
 else draction(k+1,i)=draction(k,i);
 end
 elseif d(i,j)<5&&d(i,j)>=0.2
 A=A-draction(k,j)/d(i,j);
 B=B-1/d(i,j);
 if A~=0
 draction(k+1,i)=(1-t)*draction(k,i)+t*j/B;
 else draction(k+1,i)=draction(k,i);
 end
 end
 end
 if ds(i)<5
 a=a+1;%prevx=x;%prevy=y;
 if x(k,i)>xs(k,1)
 draction(k+1,i)=atan((y(k,i)-ys(k,1))./(x(k,i)-
xs(k,1)));
 end
 if x(k,i)<xs(k,1)
 draction(k+1,i)=atan((y(k,i)-ys(k,1))./(x(k,i)-
xs(k,1)))+pi;
 end
 end
 if ds(i)>=5&&ds(i)<10
 draction(k+1,i)=pi/2;
 end
 if ds(i)>=10&&ds(i)<20
 if x(k,i)>xs(k,1)
draction(k+1,i)=-atan((y(k,i)-ys(k,1))./(x(k,i)-
xs(k,1)));
 end
 if x(k,i)<xs(k,1)
 draction(k+1,i)=-atan((y(k,i)-ys(k,1))./(x(k,i)-
xs(k,1)))+pi;
 end
 end
 v(k+1,i)=0.2;
 end
 %×ø±ê¸üÐÂ
 for i=1:n
 if ds(i)<10
 x(k+1,i)=x(k,i)+v(k+1,i)*cos(draction(k+1,i))*1.5;
 y(k+1,i)=y(k,i)+v(k+1,i)*sin(draction(k+1,i))*1.5;
 elseif ds(i)>=10
 x(k+1,i)=x(k,i)+v(k+1,i)*cos(draction(k+1,i))*1;
 y(k+1,i)=y(k,i)+v(k+1,i)*sin(draction(k+1,i))*1;
 end
 if x(k+1,i)>30
 x(k+1,i)=x(k+1,i)-30;
 elseif x(k+1,i)<0
 x(k+1,i)=x(k+1,i)+30;
 end
 if y(k+1,i)>30
 y(k+1,i)=y(k+1,i)-30;
 elseif y(k+1,i)<0
 y(k+1,i)=y(k+1,i)+30;
 end
 end
end
for i=1:k
 pause(0.2)
 plot(x(i,:),y(i,:),'.',xs(i,:),ys(i,:),'ro','markersize',12)
 title('沙丁鱼群躲避海豚捕食的二维图')
 xlabel('x');ylabel('y');%legend('沙丁鱼','海豚')
 axis([0 30 0 30])
 getframe;
end

clear all
clc
warning off
load data.mat %储存数据信息
a = randperm(30);
Train = data(a(1:25),:);
Test = data(a(26:end),:);
P_train = Train(:,3:end);
T_train = Train(:,2);
P_test = Test(:,3:end);
T_test = Test(:,2);
model = classRF_train(P_train,T_train);

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可