clc;
clear;
close all;
warning off;
addpath(genpath(pwd));
t=[4 5;16 25.8;10 45;20 55;30 65;35 55;29 31;37 26;47 27;30 31.3;31 17;14 7;35.6 13.8;26.7 22.5;21 39;38 42;5 26;28 53;20 13;10 60;26 31;54 38;7 58;12 36;30 2] %24个点,第25个点事origin
save t.mat t
load t.mat
value=[1 1 1 2 3 2 1 3 3 2 2 2 2 2 1 2 3 3 1 1 2 1 1 1]; %24个目标的价值
value=value/100;
time=zeros(1,25); %侦察UAV时间数组，里面放的是飞机走的航程，除以速度便是时间，设速度为‘1’
attacktime=zeros(1,25); %打击UAV时间数组
%把侦察过的任务对应无人机走过的航程存到该任务的一个矩阵里，当做时间，然后打击任务如果选定某任务，check一下时间是否合格，合格的话，可以打击，并存入禁忌表，不合格的话，选次概率的
%注：目的是把所有目标执行完所有任务，所以每次迭代最后所有无人机收获的总价值都一样，都是所有目标的价值之和；所以本程序考虑优先执行价值大的目标，防止无人机飞很久、打很久后攻打效率变低
%的情况出现
%%计算城市间相互距离
n=size(t,1);
D=zeros(n,n);
for i=1:n
for j=1:n
if i~=j
D(i,j)=sqrt(sum((t(i,:)-t(j,:)).^2));
else
D(i,j)=1e-2;
end
end
end
%%初始化参数
m=10; %蚂蚁个数
alpha=1; %信息素重要程度因子
beta=1; %启发函数重要程度因子
gama=2;
rho=0.3; %信息素挥发因子
Q=1.0; %总量
eta=1./D; %启发函数
tau=ones(n,n)+7.1192e-005;%信息素矩阵
iter=1; %迭代次数初始值
iter_max=80; %迭代次数最大值
length_best=zeros(iter_max,1);%每次迭代最佳路径长度(应该是一次比一次小)
length_ave=zeros(iter_max,1); %每次迭代路径平均长度
%%迭代寻找最佳路径
while iter<=iter_max
whta=cell(8,1);
lieend=zeros(8,1);
for zu=1:8
city_index=1:25; %城市来标号
table=[];
start=zeros(4,1);
temp=randperm(24);
for i=1:4
start(i)=temp(i);
end
table(:,1)=start;
j=2;
while (j<=30)
for i=1:4
if i==1 %UAV1只负责“侦察”任务
if table(1,(j-1))~=25
table1=table(1,:);
table1=[table1;table(3:4,:)];
tabu1=table1(:); %UAV1的禁忌表出来了 %25如果也在tabu1里的话，那么
allow_index1=~ismember(city_index,tabu1); %【走过的变成0，能走的为1】【若tabu=(1 4)则allow_index=(0 1 1 0 1 1 1...)】【注意：allow_index与city_index同维】
allow1=city_index(allow_index1); %把还能走的序号摘出来了（待访问的城市集合）
P1=allow1;
%计算城市的转移概率
if numel(allow1)~=0
for k=1:max(size(allow1))-1
P1(k)=(tau(table(1,(j-1)),allow1(k))^alpha)*(eta(table(1,(j-1)),allow1(k))^beta)*10000+7.1192e-004;
end
P1(max(size(allow1)))=7.1192e-005;
P1=P1/sum(P1);
[d1,ind1]=sort(P1,2,'descend');%从大到小排序是d1,对应的原序号是ind1
target1=allow1(ind1(1));
%轮盘赌法选择下一个城市
%pc1=cumsum(P1); % （p1 p1+p2 p1+p2+p3 p1+p2+p3+p4 ....）【p1<->allow(1) p2<->allow(2) ...】
%target_index1=find(pc1>=rand);
%target1=allow1(target_index1(1)); %这次返回的是allow数组中城市的真正序号
table(1,j)=target1; %把选好这个点放到路径表里面
rr=D(25,table(1,1));
time(table(1,1))=rr;
if j>2
for c=2:(j-1)
rr=rr+D(table(1,c-1),table(1,c));
end
end
rrr=rr+D(table(1,j-1),target1);%rrr就是UAV1到该点时走过的航程
time(target1)=rrr;
else
table(1,j)=25;
end
end
if table(1,(j-1))==25
table(1,j)=25;
end
end
if i==2 %UAV2只负责“打击”任务
if (table(2,(j-1))~=25)
table(2,1)=table(1,1); %设定它第一次打击的是UAV1侦察过的目标
ta2=table(1:(4*(j-1)+1)); %当前元素之前所有的元素
tabu21=[];
tabu22=[];
tabu2=[];
for y=1:24
if sum(ta2==y)==2
tabu21=[tabu21;y];
end
end %出现过两次的放在tabu21里
tabu22=setdiff(1:24,ta2); %一次都没出现的放在tabu22里
tabu2=[tabu21',tabu22]; %tabu2出来了
allow_index2=~ismember(city_index,tabu2); %【走过的变成0，能走的为1】【若tabu=(1 4)则allow_index=(0 1 1 0 1 1 1...)】【注意：allow_index与city_index同维】
allow2=city_index(allow_index2); %把还能走的序号摘出来了（待访问的城市集合）
P2=allow2;
%计算城市的转移概率
for k=1:(length(allow2)-1)
P2(k)=tau(table(2,(j-1)),allow2(k))*eta(table(2,(j-1)),allow2(k))*value(allow2(k))*10000;
end
P2(max(size(allow2)))=7.1192e-005;
P2=P2/sum(P2);
[d2,ind2]=sort(P2,2,'descend');%从大到小排序是d1,对应的原序号是ind1
target2=allow2(ind2(1)); %target2=d1(1);
%轮盘赌法选择下一个城市
%pc2=cumsum(P2); % （p1 p1+p2 p1+p2+p3 p1+p2+p3+p4 ....）【p1<->allow(1) p2<->allow(2) ...】
%target_index2=find(pc2>=rand); %选中那个概率较大的选中的点，返回的是allow数组中的序号
%target2=allow2(target_index2(1)); %这次返回的是allow数组中城市的真正序号
%table(2,j)=target2; %把选好这个点放到路径表里面
oo=D(25,table(2,1));
attacktime(table(2,1))=oo;
if j>2
for c=2:(j-1)
oo=oo+D(table(2,c-1),table(2,c));
end
end
ooo=oo+D(table(2,j-1),target2);%ooo就是UAV2到该点时走过的航程
if numel(d2)>5
u=2;
while (ooo>time(target2)+20 & u<6)
target2=allow2(ind2(u));
ooo=oo+D(table(2,(j-1)),target2);
u=u+1;
end
end
table(2,j)=target2;
attacktime(target2)=ooo;
end
if table(2,(j-1))==25
table(2,j)=25;
end
end
if i==3 %UAV3是“察打”任务
if table(3,(j-1))~=25
ta3=table(1:(4*(j-1)+2));
tabu3=[];
tabu3c=[];
for y=1:24
if sum(ta3==y)==2
tabu3=[tabu3;y];
end
end %出现两次的放在tabu3里
for y=1:24
if sum(ta3==y)==1
tabu3c=[tabu3c;y];
end
end %tabu3c是待打的任务，已侦查完的任务
allow_index3=~ismember(city_index,tabu3); %【走过的变成0，能走的为1】【若tabu=(1 4)则allow_index=(0 1 1 0 1 1 1...)】【注意：allow_index与city_index同维】
allow3=city_index(allow_index3); %把还能走的序号摘出来了（待访问的城市集合）
P3=allow3;
%计算城市的转移概率
for k=1:(length(allow3)-1)
%if ismember(allow3(k),tabu3c)==1
h=table(3,(j-1))
P3(k)=(tau(table(3,j-1),allow3(k))^alpha)*(eta(table(3,(j-1)),allow3(k))^beta)*value(allow3(k))*10000+7.1192e-005;%这是要打的，需要价值
%else
%P3(k)=(tau(table(3,(j-1)),allow3(k))^alpha)*(eta(table(3,(j-1)),allow3(k))^beta)*100+7.1192e-005;%这些是待侦察的，没有价值
%end
end
P3(max(size(allow3)))=7.1192e-009;
P3=P3/sum(P3);
[d3,ind3]=sort(P3,2,'descend');%从大到小排序是d1,对应的原序号是ind1
target3=allow3(ind3(1));
%轮盘赌法选择下一个城市
%pc3=cumsum(P3); % （p1 p1+p2 p1+p2+p3 p1+p2+p3+p4 ....）【p1<->allow(1) p2<->allow(2) ...】
%target_index3=find(pc3>=rand); %选中那个概率较大的选中的点，返回的是allow数组中的序号
%target3=allow3(target_index3(1)); %这次返回的是allow数组中城市的真正序号
%table(3,j)=target3; %把选好这个点放到路径表里面
ww=D(25,table(3,1));
time(table(3,1)
