%% 蚁群算法求解水下无线网络节点的最大覆盖面积问题
close all;clc;
%% 仿真数据, shengdu代表首层水面高度；numebr代表节点数量；r_tongxin表示通信距离；
number=10; shengdu=400; fanwei=70; r_tongxin=70;
C=gallery('uniformdata',[number 3],10).*fanwei;
% C=rand(number,3).*fanwei;
C(:,3)=shengdu;
%% 仿真数据预处理，保证节点间连通率为1--------------------这个地方有待完善，对于不连通的点可以进行处理！！！！！！！
result=count_liantonglv(C,C(1,:),r_tongxin);
%% 仿真数据按照水下深度加密并展示 ------暂时猜测由于cengshu的原因，导致不同深度的节点不具备通信条件，即单层的高度应该小于节点通信半径，
cengshu=20; C=repmat(C,cengshu,1); % 将首层水面节点的坐标按照规定的层数进行复制；
for i=number:number:size(C,1) % 将每一层的节点高度进行调节，从最高按层降低至0；
C(i-number+1:i,3)=shengdu-(i/number-1)*shengdu/(cengshu-1);
end
figure(1);
plot3(C(:,1),C(:,2),C(:,3),'k.','markersize',26);xlabel('x 轴');ylabel('y 轴');zlabel('z 轴');
title('城市分布图');grid on;
% M为问题的规模 M个城市
[M,N] = size(C);
% 用来记录任意两个城市之间的距离
distance = zeros(M,M);
for m=1:M
for n=1:M
distance(m,n) = sqrt(sum((C(m,:)-C(n,:)).^2));
end
end
%% 参数预定义
m = M*1.2; % 蚂蚁的个数 一般取10-50
alpha = 1; % 计算转移概率中信息素的重要程度 一般取【1,4】
beta = 1; % 计算转移概率中启发式英子的重要程度 一般取【3,5】
rho = 0.3; % 相邻时间的信息素蒸发系数
G = 100; % 迭代次数上限
Q = 10; % 信息素增加系数
Eta = distance; % 启发式因子，这里直接取两点之间的距离
Tau = ones(M,M); % 信息素矩阵 存储着每两个城市之间的信息素的数值
Tabu = zeros(m,M); % 禁忌表，记录每只蚂蚁走过的路程
gen = 1;
R_best = zeros(G,M); % 各代的最佳路线,节点总数应该只有水面一层的数量；
L_best = 0.*ones(G,1); % 每一代的无效体积覆盖率-初始假设为无穷大
%% 开始迭代计算
while gen<G
% 将m只蚂蚁随机放到M个城市上，每只蚂蚁都是从水面开始的，保证和sink链接
% random_pos = [];
% for i=1:(ceil(m*cengshu/M))
% random_pos = [random_pos,randperm(number)];
% end
% Tabu(:,1) = (random_pos(1,1:m))';
Tabu(:,1) = 1;
%% 访问第几个城市，相同的XY只能访问一次；
for i=2:(number)
for j=1:m % m只蚂蚁
visited = Tabu(j,1:(i-1)); % 访问第i个城市前，第j个蚂蚁已经访问过的城市；
visit_P = zeros(1,(M-(i-1)*cengshu)); % 蚂蚁j访问剩下的城市的概率
unvisited= (1:1:M);
we1=ismember(C(:,1),C(visited,1)); %将a中第一列等于b的为1，不等于的为0；
we2=ismember(C(:,2),C(visited,2)) ;
we11=find(we1==1&we2==1);
unvisited(we11)=[];
%% 计算待选择城市的概率;退档机制：由visited的end开始逐步向前返回，直至有点可以前进；
for num_visited=size(visited,2):-1:1
if ~isempty(find(sqrt(sum((C(unvisited,:)-repmat(C(visited(num_visited),:),size(unvisited,2),1)).^2,2))<r_tongxin, 1))% 判断当前经过点，与未经过点存在路径则计算概率；
for k=1:length(unvisited) % Tau(visited(end),unvisited(k))访问过的城市的最后一个与所有未访问的城市之间的信息素
visit_P(k) = ((Tau(visited(num_visited),unvisited(k)))^alpha)*(Eta(visited(num_visited),unvisited(k))^beta);
% 假如这两个点无法通过多跳实现通信，表示访问概率为0
if(norm(C(visited(num_visited),:)-C(unvisited(k),:))>r_tongxin)
visit_P(k)=0;
end
end
break;
else
if num_visited==1
fprintf('what is wrong');
end
end
end
temm=visit_P;
visit_P = visit_P/sum(visit_P); % 访问每条路径的概率的大小,如果概率为0，需要解决selected为空集的问题。
%% 按照概率选择下一个要访问的城市
% 这里运用轮盘赌选择方法 这里也可以选择选择概率最大的路径去走， 这里采用轮盘赌选择法。
Pcum = cumsum(visit_P); % 求得visit_P的累加值
selected = find(Pcum>=rand); % 求得大于随机数的概率在数组中的位置
to_visited = unvisited(selected(1)); % 确定第j个蚂蚁将要访问的下一城市
Tabu(j,i) = to_visited; % 添加到禁忌表
end
end
% if gen>=2
% Tabu(1,:) = R_best(gen-1,:);
% end
%% 记录m只蚂蚁迭代的最佳路线的体积有效覆盖率，
L = zeros(1,m);sum_ratio=0;
for i=1:m
L(i) = series_overlapping(C(Tabu(i,1:number),:),r_tongxin,[0 fanwei;0 fanwei;0 shengdu]);
end
L_best(gen) = max(L); % 记录每一代中路径的无效覆盖率最小值；
pos = find(L==L_best(gen)); % pos代表了最短路径是第几只蚂蚁；
R_best(gen,:) = Tabu(pos(1),:); % 最优的路径
%% 更新信息素的值
Delta_Tau = zeros(M,M);
for i=1:m % m只蚂蚁
for j=1:(number-1) % M座城市,实际蚂蚁单次只走完number个；
Delta_Tau(Tabu(i,j),Tabu(i,j+1)) = Delta_Tau(Tabu(i,j),Tabu(i,j+1)) + Q/(1/L(i)); % m只蚂蚁的信息素累加 这里采用的是论文中ant-cycle模型
end
% Delta_Tau(Tabu(i,M),Tabu(i,1)) = Delta_Tau(Tabu(i,M),Tabu(i,1)) + Q/L(i);%如果注释则表示加上终点返回起点的那一段路径的信息素
end
Tau = (1-rho).*Tau+Delta_Tau; % 更新路径上的信息素含量
%% 禁忌表清零
Tabu = zeros(m,M);
%% 绘制最佳路线
for i=1:(number-1)
plot3([C(R_best(gen,i),1),C(R_best(gen,i+1),1)],[C(R_best(gen,i),2),C(R_best(gen,i+1),2)],[C(R_best(gen,i),3),C(R_best(gen,i+1),3)],'bo');
hold on;
end
hold on;
grid on;
% 绘制最佳路线的起点、终点
% plot3(C(R_best(gen,1),1),C(R_best(gen,1),2),C(R_best(gen,1),3),'ro');
% plot3(C(R_best(gen,M),1),C(R_best(gen,M),2),C(R_best(gen,M),3),'ro');
title(['最高覆盖率：',num2str(max(L_best(gen)))]);
hold off;
pause(0.05);
gen = gen+1;
end
%% 绘图
L_plot=[];
for i=1:length(L_best)
L_plot=[L_plot;max(L_best(1:i))];
end
figure(2);
plot(L_plot,'ko-');
grid on;
title('路径长度变化曲线');
xlabel('迭代次数');
ylabel('路径长度数值');
% 计算全部球体的有效体积比例，即【（全部球体的实际覆盖体积/全部球体的理论总体积】；
% sphere ： [n x 3]，代表n个球体的球心坐标 x，y，z；需要注意的是这里的xyz必须都大于0；
% r_tongxin ： 球体的半径，这里等于通信半径；
% quyu ： [x0 x1;y0 y1;z0 z1]代表节点所处的区域；
% 由于实验区域是以[0 0 0]作为起点的长方体区域，因此区域内坐标全部是大于0的，这样便于设定边界；
function sum_ratio=series_overlapping(sphere,r_tongxin,quyu)
% 测试数据，可以注释
%
% sphere= [41.76061184 68.54127777 400
% 15.29626893 183.1294293 400
% 52.18824665 128.58125 400
% 23.28042305 28.98108094 400
% 115.5835083 168.5242946 400
% 98.9211924 94.68295859 400
% 65.6103901 83.84985899 400
% 177.7892674 179.3786168 400
% 98.6246862 13.36017467 400
% 190.1825888 66.2612519 400
% 131.311037 52.59624236 400
% 104.4405789 124.7195736 400
% 154.1250962 125.8762652 400
% 56.89105971 174.6295173 400
% 138.1844883 73.15866107 400
% 174.5788959 34.34024643 400
% 64.10875284 146.0318037 400
% 7.237784561 91.72803281 400
% 164.6614527 133.5909421 400
% 3.18172403 43.11732495 330
% 0.301033791 78.88285038 330
% 182.0978432 155.2369154 330
% 14.72034335 160.4334129 330
% 185.693325 174.4790427 330
% 88.71079591 171.6445411 330
% 189.5333513 146.4151184 330
% 135.9935366 39.00809681 330
% 4.861887192 57.21357568 330
% 158.3052095 37.78816506 330
% 124.264238 11.38644003 330
% 89.85630654 192.5063974 260
% 84.61629932 182.7453295 260
% 137.1361678 19.7331577 260
% 0.932791444 143.7630805 260
% 125.6079739 8.719560836 260
% 192.8158001 6.22558221 190
% 9.054379029 157.1439911 190
% 71.26686176 189.490607 190
% 69.38345592 193.9794298 120
% 196.8496581 3.019990816 120]
% r_tongxin=70;
% quyu=[0 200;0 200;0 400];
%% 初始化参数
buchang=5; % 初始化小正方形的步长；
row=size(sphere,1);% 看一下球体的数量；
posi=[];
for i=1:row
posi=[posi;lapping(sphere(i,:),r_tongxin,quyu)]; % 统计每个球体覆盖的正方体中心坐标，汇总存入posi；
end
result=unique(posi,'rows');%选取独一的覆盖正方体；
total=quyu(1,2)*quyu(2,2)*quyu(3,2)/(buchang^3);
sum_ratio=size(result,1)/total;% 这里ratio一直等于1，原因是每个球体对应的正方体坐标在小数点后都不相同；不对，不止小数点，个位数都应该是0开始才能匹配上；
end
% 输出单个球体所占的正方体坐标；这里的buchang最好为5或者10，不然会有bug；
% sphere ： [1 x 3]，代表n个球体的中心坐标 x，y，z；
% r_tongxin ： 球体的半径，这里等于通信半径
% buchang ： 小正方形的步长（由于是将空间分为较多个小正方体来统计）；
function zuobiao=lapping(sphere,r_tongxin,quyu)
buchang=5; % 小正方体边长默认为5m；
jia=2;
% 将球体所在的正方形网格区域划分为n个边长为buchang的小正方体,total则表示正方形的总数；
nn=ceil(2*r_tongxin/buchang)+jia;% 这里加jia是由于网格分布需要从球心坐标向下取整的坐标来划分，会损失右部边缘部分；
total=nn^3;
XYZ_INTI=floor(sphere(1,:))-floor(rem(sphere(1,:),10));
% 在此基础上获取全部正方体的中心坐标存入 y_matrix ，从而判断首个球体的覆盖正方体；
y_matrix=zeros(total,3);
x1 = XYZ_INTI(1)-r_tongxin+buchang/2:buchang:XYZ_INTI(1)-r_tongxin+buchang/2+buchang*(nn-1);
y1 = repmat(x1,nn^2,1);
z1 = y1(:);
x2 = XYZ_INTI(2)-r_tongxin+buchang/2:buchang:XYZ_INTI(2)-r_tongxin+buchang/2+buchang*(nn-1);
y2 = repmat(x2,nn,1);
z2 = y2(:);
z2 = repmat(z2,nn,1);
x3 = XYZ_INTI(3)-r_tongxin+buchang/2:buchang:XYZ_INTI(3)-r_tongxin+buchang/2+buchang*(nn-1);
z3 = x3(:);
z3 = repmat(z3,nn^2,1);
y_matrix(:,1:3)=[z1,z2,z3];
% 依次判其余球体是否覆盖首个球体的正方形，覆盖结果保存至sph_matrix,其中每列表示一个球体的情况，覆盖为1，不覆盖为0。
cc=buchang/2;
list=y_matrix(:,1)-cc>=quyu(1,1)&y_matrix(:,1)+cc<=quyu(1,2)&y_matrix(:,2)-cc>=quyu(2,1)&y_matrix(:,2)+cc<=quyu(2,2)...
&y_matrix(:,3)-cc>=quyu(3,1)&y_matrix(:,3)+cc<=quyu(3,2);% list代表在限定区域范围内的坐标位置
last_matrix=y_matrix(list,:);% 区域范围内的坐标
list2=sqrt(sum((last_matrix-repmat(sphere,size(last_matrix,1),1)).^2,2))<=r_tongxin;
zuobiao=last_matrix(list2,:);
end
% function result=count_liantonglv(x_jiedian,x_sink,cr)
% result ： 表示连通率
% x_jiedian[n 3] ： 表示水下网络区域的全部节点坐标
% x_sink[3 1] ： 表示水面sink节点的坐标
% cr ： 表示节点的通信半径
function result=count_liantonglv(x_jiedian,x_sink,cr)
% % 数据以下为例：
% x_sink=[250 250 500];
% x_jiedian=gallery('uniformdata',[400 3],4).*500;
% cr=77;
% 统计 (一跳)通信的节点，新建数组liantong，对应每个节点的连同结果，1代表连同，2代表失联；
number=size(x_jiedian);
liantong=zeros(number(1),1)
sum=0;
for i=1:number(1)
% liantong(i)=ganzhi_moxing(x_jiedian(i,:),x_sink,cr,1)
liantong(i)=norm(x_jiedian(i,:)-x_sink)<50;
if(liantong(i)==1)
sum=sum+1;
end
end
result=sum/number(1);
%如果不存在(一跳)节点，则停止统计多跳节点；
true=1;
if(result==0)
true=0;
end
%统计(多跳)通信的节点
while(true)
success_epoch=0; % success_epoch 表示单次循环中成功连通的节点数量；
for i=1:number(1)
if(liantong(i)==0)% i 表示还没有连通的节点
for j=1:number(1)
if(liantong(j)==1)% j 表示已经连通的节点
liantong(i)=ganzhi_moxing(x_jiedian(i,:),x_jiedian(j,:),cr,1);% 判断未连通节点中，哪些节点可以连通已知节点
if(liantong(i)==1)
success_epoch=success_epoch+1;
break;
end
end
end
end
end
if(success_epoch==0)
break;
end
end
% 更新连通率
sum=0;
for i=1:number(1)
if(liantong(i)==1)
sum=sum+1;
end
end
result=sum/number(1);
end 0
点赞
收藏
分享
基于蚁群算法的水下传感网节点部署
相关推荐
0 条评论
