【优化求解】基于matlab虚拟力算法求解无线网络传感覆盖问题【含Matlab源码 1187期】

199 阅读4分钟

一、简介

虚拟力覆盖算法应用于传感器网络近几年引起了越来越多的人的关注,这首先由Zou等人将虚拟力引入到传感器网络,很明显传感器节点和机器人存在着很多共同之处。所以他们在虚拟力应用于机器人领域的基础上,把传感器也看作是势场中的粒子,由于地球是个巨大的磁场和地球周围存在着电磁场,可以把节点看作是存在势场中的粒子,这样传感器节点之间就存在着相互的联系。

二、源代码

clear all
clc
N=40;%传感器节点个数
XMAX=900;%区域总长度
XMIN=100;
YMAX=800;%区域总宽度
YMIN=100;%区域总宽度
x=zeros(N,2);
x(:,1)=(XMAX-XMIN)*rand(N,1)+XMIN;%使节点坐标随机分布
x(:,2)=(YMAX-YMIN)*rand(N,1)+YMIN;
figure,
xm=[XMIN YMIN;XMIN YMAX;XMAX YMAX;XMAX YMIN];%区域的四个顶点
fill(xm(:,1),xm(:,2),[0.8,0.8,0.8]);
hold on
plot(x(:,1),x(:,2),'r.','linewidth',5);
r=90; %传感器节点的感知半径
w=0:pi/50:2*pi;
for i=1:N
    x1=x(i,1)+r*cos(w);
    y1=x(i,2)+r*sin(w);%表示圆每一个角度对应一个点,然后把这些点连接起来构成圆
    hold on
    plot(x1,y1,'g');   %画出各个传感器节点的感知范围
    text(x(i,1)+3,x(i,2),['\fontsize{8}\rm',num2str(i)]);
    hold on
    fill(x1,y1,'b')
end
axis([0 1000 0 900]);%设定坐标范围
xlabel('X/m');ylabel('Y/m');
%legend('要监测的区域',['传感器节点(','\fontsize{12}\bf',num2str(N),'\fontsize{10}\rm个)的位置'],' 传感器节点的感知范围');
hold on
plot([XMIN XMAX],[YMIN YMIN],'k','linewidth',1.5);
hold on
plot([XMIN XMIN],[YMIN YMAX],'k','linewidth',1.5);
hold on
plot([XMAX XMAX],[YMIN YMAX],'k','linewidth',1.5);
hold on
plot([XMIN XMAX],[YMAX YMAX],'k','linewidth',1.5);%把四个顶点连接起来组成一个监测区域
%--------------------------------虚拟力算法---------------------------------
%------------先对区域进行离散化---------
deta=2;%网络大小
x1=XMIN:deta:XMAX;
y1=YMIN:deta:YMAX;
[xx,yy]=meshgrid(x1,y1);
[m,n]=size(xx);
K=m*n; %总的网格点数目
xx1=reshape(xx,K,1);%网格点的横坐标
yy1=reshape(yy,K,1);%网格点的纵坐标
%hold on
%plot(xx1,yy1,'g*')
%----------------计算起初的网络覆盖率-----------------------------------------------------
[no_cover,summ,k1]=compute_cover(xx1,yy1,x,r); 
%no_cover存储没有被覆盖的格点位置,k1为求被覆盖的格点,summ为被覆盖的格点数;
q(1,1)=summ/K;
%-------------------------------------------------------
figure,
fill(xm(:,1),xm(:,2),[0.8,0.8,0.8]); %填充监测区域
hold on
plot(x(:,1),x(:,2),'ro','markerfacecolor','r','linewidth',3); %传感器节点位置
for i=1:N
    text(x(i,1)+3,x(i,2),['\fontsize{8}\rm',num2str(i)]); %标出传感器节点位置
end
axis([0 1000 0 900]);
xlabel('X/m');ylabel('Y/m');
%----------------------
R=2*r;    %传感器节点的通信半径
maxiter=100; %最大迭代次数
max_step=2.5; %传感器节点移动的最大步长(在格点作用下的最大步长)
max_sensor=3.5; %传感器节点移动的最大步长(在传感器节点作用下的最大步长)
%----------------------
kp=1;
xp{kp,1}=x;
for t=1:maxiter
    %k3=1;
    F=0;
    tx_old=x; %前一次迭代的传感器节点位置
    for i=1:N
        k2=1;
        flag1=0;
        for k=1:(k1-1)
            dik=sqrt((x(i,1)-no_cover(k,1))^2+(x(i,2)-no_cover(k,2))^2); %x为传感器节点的位置,no_cover是未被网络覆盖的网格点的位置
            if ((dik>r)&(dik<=R))   %R为通信半径,r为感知半径
                 F(k2,1)=no_cover(k,1)-x(i,1);
                 F(k2,2)=no_cover(k,2)-x(i,2);
                k2=k2+1;
                %flag1=flag1+1;
            end
        end
        Fx=0;
        Fy=0;
        for j=1:(k2-1)
           Fx=Fx+F(j,1); %水平力
           Fy=Fy+F(j,2); %垂直力
        end
        Fxy=sqrt(Fx^2+Fy^2);
        if Fxy==0
            x(i,1)=x(i,1)+0;
            x(i,2)=x(i,2)+0;
        else
            x(i,1)=x(i,1)+Fx/Fxy*max_step*exp(-1/Fxy);%节点位置更新
            x(i,2)=x(i,2)+Fy/Fxy*max_step*exp(-1/Fxy);
        end
        %if flag1~=0
        %   k3=k3+1;
       %end
    end
    %-----------------------传感器之间的斥力-------
    Fs=0;
    for i=1:N
        ks=1;
        for j=1:N
            if i~=j
                ds=sqrt((x(i,1)-x(j,1))^2+(x(i,2)-x(j,2))^2);
                if ds<(1.4*r)
                    Fs(ks,1)=x(i,1)-x(j,1);
                    Fs(ks,2)=x(i,2)-x(j,2);
                    ks=ks+1;
                end
            end
        end
        Fsx=0;
        Fsy=0;
        for j=1:(ks-1)
           Fsx=Fsx+Fs(j,1); %水平力
           Fsy=Fsy+Fs(j,2); %垂直力
        end
        Fsxy=sqrt(Fsx^2+Fsy^2);
        if Fsxy==0
            x(i,1)=x(i,1)+0;
            x(i,2)=x(i,2)+0;
        else
            x(i,1)=x(i,1)+Fsx/Fsxy*max_sensor*exp(-1/Fsxy);
            x(i,2)=x(i,2)+Fsy/Fsxy*max_sensor*exp(-1/Fsxy);
        end
        %-----------------不能出边界------
        if x(i,1)<XMIN
            x(i,1)=XMIN;
        end
        if x(i,1)>XMAX
            x(i,1)=XMAX;
        end
         if x(i,2)<YMIN
            x(i,2)=YMIN;
        end
        if x(i,2)>YMAX
            x(i,2)=YMAX;
        end
    end
    %-----------------------------------------------------------
    [no_cover,summ,k1]=compute_cover(xx1,yy1,x,r); %计算覆盖覆盖范围
    %------------------
    for i=1:N
        hold on
        plot([tx_old(i,1) x(i,1)],[tx_old(i,2) x(i,2)],'-'); %画传感器节点的运动轨迹图
    end
    %-------------------
    
    if t==(kp*round(maxiter/4)) %为了得到总迭代步长中每四分之一时的情况
        kp=kp+1;
        xp{kp,1}=x;
    end
    q(t+1,1)=summ/K;
    if q(t+1,1)>=1
       break;
    end
end
function [no_cover,summ,k1]=compute_cover(xx1,yy1,x,r)
%% (xx1,yy1)分别为网格点的横纵坐标;
%% x为传感器节点的坐标, r为传感器节点的感知半径
%% no_cover存储没有被覆盖的格点位置,k1为求被覆盖的格点,summ为被覆盖的格点数;
no_cover=0; %存储没有被覆盖的格点位置
flag=0;
k1=1;
summ=0;
K=length(xx1);
[N,N1]=size(x);
for k=1:K   %K为总的网格点数目
    flag=0;
    for j=1:N   %寻找不在覆盖区域内的网格点,条件为:x-r<=xx1<=x+r
        D=sqrt((xx1(k,1)-x(j,1))^2+(yy1(k,1)-x(j,2))^2); %网格点与传感器节点间的距离
        if (D<=r)  %在圆盘范围内的网格点作为可覆盖的格点
           flag=flag+1;
        end
        if flag==1
            summ=summ+flag;
            break;
        end
    end
    if flag==0
        no_cover(k1,1)=xx1(k,1); %存储没有被覆盖的格点位置
        no_cover(k1,2)=yy1(k,1);
        %text(no_cover(k1,1)+1,no_cover(k1,2),['\fontsize{8}\rm',num2str(k1)]);
        k1=k1+1;  
    end
end
k1=k1-1;

三、运行结果

在这里插入图片描述 在这里插入图片描述 在这里插入图片描述 在这里插入图片描述

四、备注

版本:2014a