function shexing%%目標在監控區內運動軌跡為蛇形曲線時的定位%%M=100;                                                                       %蒙特卡羅實驗次數N0=50;%%初始化監控區域的長、寬、無線傳感器個數和目標在監控區內的運動時間%%a=100;%input('請輸入監控區域的長\n');                                       %初始化監控區域的長b=100;%input('請輸入監控區域的寬\n');                                       %初始化監控區域的寬s=100;%input('請輸入無線傳感器個數\n');                                     %初始化監控區域中的無線傳感器個數t=50;%input('請輸入目標在監控區內的運動時間（單位：秒）\n');                  %選擇觀察某一監控時刻（單位：秒）dn=0;%input('請輸入聯盟中損壞節點個數\n');%初始化目標軌跡n1=10;n2=25;n3=a-25;n4=a-10;dt=(a-20)/t;for i=1:t    if i*dt+n1<=n2        Xt(i)=i*dt+n1;        Yt(i)=b/2;    elseif i*dt+n1<=n3&i*dt+n1>n2        Xt(i)=i*dt+n1;        Yt(i)=(0.8*b/2)*sin((Xt(i)-n2)*2*pi/(a-50))+b/2;    elseif i*dt+n1<=n4&i*dt+n1>n3        Xt(i)=i*dt+n1;        Yt(i)=b/2;    endendZt=zeros(1,t);                                                             %目標的垂直高度%%初始化無線傳感器的監控范圍和傳輸范圍%%Ls=20;                                                                     %無線傳感器的監測范圍Lt=20;                                                                     %無線傳感器的傳輸范圍alfa=2;%%初始化目標定位坐標的數組%%Xr=zeros(1,t);                                                             %目標定位的橫坐標Yr=zeros(1,t);                                                             %目標定位的縱坐標Xrd=zeros(1,t);                                                            %傳感器受損情況下的目標定位的橫坐標Yrd=zeros(1,t);                                                            %傳感器受損情況下的目標定位的縱坐標Xm=zeros(1,t);                                                             %蒙特卡羅實驗得到的橫坐標Ym=zeros(1,t);                                                             %蒙特卡羅實驗得到的縱坐標Xmd=zeros(1,t);                                                            %傳感器受損情況下的蒙特卡羅實驗得到的橫坐標Ymd=zeros(1,t);                                                            %傳感器受損情況下的蒙特卡羅實驗得到的縱坐標D0=zeros(1,t);                                                             %距離誤差D0d=zeros(1,t);                                                            %距離誤差D1=zeros(1,t);                                                             %距離誤差D1d=zeros(1,t);                                                            %距離誤差%初始化一些傳感器信息p=0.0;                                                                     %無線傳感器信息素在單位時間內的揮發后剩余量de_mark=zeros(1,s);                                                        %傳感器探測到目標的標志位D=zeros(1,s);                                                              %初始化各傳感器一個位置常量E0=10000;                                                                  %傳感器初始能量for kk=1:N0    Xm=zeros(1,t);    Ym=zeros(1,t);    Xmd=zeros(1,t);    Ymd=zeros(1,t);    for k=1:M        X=rand(1,s)*a;                                                         %隨機產生各傳感器橫坐標        Y=rand(1,s)*b;                                                         %隨機產生各傳感器縱坐標        neigh=zeros(s,s);                                                      %傳感器的鄰近節點表,儲存節點序號        [D,neigh]=dec(s,X,Y,D,Lt);                                             %計算位置常量        w0=eps;                                                                %無線傳感器的目標存在概率之和        wd=eps;                                                                %無線傳感器的目標存在概率之和        N=zeros(1,s);                                                          %初始化各傳感器的信息素值        E=ones(1,s)*E0;                                                        %傳感器能量        for i=1:t            N=p*N;                                                             %傳感器信息素單位時間內的衰減            w=zeros(1,s);                                                      %無線傳感器在某一時刻的目標存在概率            de_mark=zeros(1,s);                                                %傳感器探測到目標的標志位            dp=zeros(s,s);                                                     %傳感器信息素增量表            S=zeros(1,s);                                                      %初始化各傳感器的工作狀態。1代表活躍狀態，0代表睡眠狀態。            Sd=zeros(1,s);                                                     %傳感器受損情況下的各傳感器的工作狀態。            for j=1:s                                                          %傳感器信息素增量                d=sqrt((X(j)-Xt(i))^2+(Y(j)-Yt(i))^2+Zt(i)^2)+eps;                if d<=Ls                    de_mark(j)=1;                    N(j)=N(j)+(Ls/(d))^alfa;                end                if de_mark(j)==1                                               %如果傳感器探測到目標則向外傳播信息素                    for m=1:s                        if m~=j&neigh(j,m)==1                            dp(j,m)=dp(j,m)+rat(j,m,X,Y,Lt,alfa)*N(j);                            E(j)=E(j)-1;                                       %規定傳感器廣播一次信息素消耗1單位能量，其余行為不計消耗                        end                    end                end            end            [S,Sd,N]=judg(N,de_mark,D,E,E0,s,Ls,alfa,dn,neigh,dp);            for j=1:s %目標定位                if S(j)==1                    w(j)=2*atan(N(j))/pi;                    Xr(i)=Xr(i)+X(j)*w(j);                    Yr(i)=Yr(i)+Y(j)*w(j);                    w0=w0+w(j);                end                if Sd(j)==1                    w(j)=2*atan(N(j))/pi;                    Xrd(i)=Xrd(i)+X(j)*w(j);                    Yrd(i)=Yrd(i)+Y(j)*w(j);                    wd=wd+w(j);                end            end            Xr(i)=Xr(i)/w0;            Yr(i)=Yr(i)/w0;            w0=eps;            Xrd(i)=Xrd(i)/wd;            Yrd(i)=Yrd(i)/wd;            wd=eps;        end        Xm=Xm+Xr;        Ym=Ym+Yr;        Xr=zeros(1,t);                                                         %目標定位的橫坐標        Yr=zeros(1,t);                                                         %目標定位的縱坐標        Xmd=Xmd+Xrd;        Ymd=Ymd+Yrd;        Xrd=zeros(1,t);                                                        %目標定位的橫坐標        Yrd=zeros(1,t);                                                        %目標定位的縱坐標    end    Xm=Xm/M;    Ym=Ym/M;    Xmd=Xmd/M;    Ymd=Ymd/M;    D0=zeros(1,t);    D0d=zeros(1,t);    for i=1:t        D0(i)=sqrt((Xt(i)-Xm(i))^2+(Yt(i)-Ym(i))^2);        D0d(i)=sqrt((Xt(i)-Xmd(i))^2+(Yt(i)-Ymd(i))^2);    end        D1=D1+D0;    D1d=D1d+D0;endD1=D1/N0;D1d=D1d/N0;figureplot(Xt,Yt)axis([0,a,0,b])xlabel('單位：米')ylabel('單位：米')legend('目標軌跡')figureplot(Xt,Yt)axis([0,a,0,b])xlabel('單位：米')ylabel('單位：米')hold onplot(Xm,Ym,'r')legend('目標實際軌跡','定位軌跡')figuredmax=max(D1);ddmax=max(D1d);plot(D1,'-')axis([1,t-1,0,ddmax*1.1])xlabel('單位：秒')ylabel('定位誤差/米')if dn>0    hold on    plot(D1d,'-.')endclear allfunction r=rat(i,j,X,Y,Lt,alfa)d=sqrt((X(i)-X(j))^2+(Y(i)-Y(j))^2);r=(1-d/Lt)^alfa;function [Q,neig]=dec(s,X,Y,D,Lt)neig=zeros(s,s);for i=1:s    h0=0;    k=0;    h=zeros(1,s);    for m=1:s        d=sqrt((X(i)-X(m))^2+(Y(i)-Y(m))^2);        if d<=2*Lt           h(m)=d;           k=k+1;        end        if d<=Lt           neig(i,m)=1;        end    end    for m=1:s    h0=h0+h(m);    end    D(i)=h0/(2*(k-1)*Lt);endQ=D;function [Q0,Qd,N1]=judg(N,de_mark,D,E,E0,s,Ls,alfa,dn,neigh,dp)Q0=zeros(1,s);Qd=zeros(1,s);N1=zeros(1,s);count=0;for i=1:s    if D(i)*N(i)*E(i)/E0>=(Ls/(Ls-2))^alfa&de_mark(i)==1        Q0(i)=1;        count=count+1;    endendQd=Q0;if dn>0    if count<dn        Qd=zeros(1,s);    else        for i=1:dn            for j=1:s                if Qd(j)==1;                    Qd(j)=0;                    for m=1:s                        if neigh(j,s)==1                            Qd(m)=1;                            N(m)=N(m)+dp(j,m);                        end                    end                    break                end            end        end    endendN1=N;