matlab 任务分配粒子群算法

任务模型基于此篇论文

Xiaowei Jiang, Qiang Zhou, Ying Ye

Method of task assignment for UAV based on particle swarm optimization in logistics[P]

clear;clc;
%%
%    [ x    ，  y   ， need  ，ai   ， bi   ，service_time]
global T
T = [[0.00  , 0.00  , 0.00  , 6.00  , 12.0  ,   0  ]      %0分布矩阵
     [43.6  , 4.80  , 0.50  , 7.00  , 9.00  , 15/60]      %1
     [27.2  , 22.5  , 0.65  , 7.00  , 9.00  , 30/60]      %2
     [39.9  , 46.8  , 0.30  , 7.00  , 8+2/3 , 15/60]      %3
     [42.8  , 32.3  , 0.23  , 7+1/3 , 8.00  , 30/60]      %4
     [15.8  , 37.1  , 0.50  , 8.00  , 9.00  , 30/60]      %5
     [6.7   , 25.4  , 0.55  , 7.00  , 9+1/6 , 30/60]      %6
     [22.0  , 27.2  , 0.42  , 7.00  , 8.5   , 15/60]      %7
     [34.4  , 16.1  , 0.51  , 7.00  , 8.00  , 15/60]      %8
     [43.2  , 29.3  , 0.40  , 7.00  , 9.00  , 15/60]      %9
     [45.4  , 13.4  , 0.35  , 8.00  , 9+1/6 , 15/60]      %10
     ];
global demand   %uav总的搬运能力
demand = 2;  
global L_max    %uav行驶的最大距离
L_max = 1000000;
global L_c      %uav行驶的单位代价 yuan/km
L_c = 1;
global v        %uav飞行速度    km/h
v = 50;
global y        %迟到惩罚代价   yuan/h
y = 1.5;
n = 10;                         %n个任务
m = 3;                          %m个无人机 
dis = ones(n + 1, n + 1);           %记录无人机与任务两两之间欧式距离
for i = 1:n + 1
    for j = 1:n + 1
        dis(i, j) = ((T(i, 1) - T(j, 1))^2 + (T(i, 2) - T(j, 2))^2)^0.5;
    end
end
 
NP = 100;                   % 种群规模
w  = 0.5;                   % 惯性权重
c1 = 1.5;                   % 自我认知学习因子
c2 = 1.5;                   % 社会认知学习因子
g = ones(NP, n);            % 粒子群
best = 0;                   % 最佳迭代次数
num_max = 500;             % 最大进化代数
 
% 粒子初始化
for i = 1:NP
   % 粒子随机生成排序
    g(i, :) = randi([1,m],1,n);
end
 
gb = ones(1, 1);   %存储每次迭代得全局最优，用于plot        
num = 1;                % 进化代数
cost = zeros(NP,1);     % 存放个体的适应值
global_best_fit = 10^5;    % 存放全局最优值
global_best = g(1, :);        % 存放全局对应最优值的解
individual_best = g;    % 存放个体最优解
individual_best_fit = cost;   % 存放个体的最优适应值

for i = 1:NP            % 计算初始种群中的最优值和最优值对应的解
    individual_best_fit(i) = fitness(g(i, :), dis, n, m);
    if(individual_best_fit(i) <= global_best_fit)
        global_best_fit = individual_best_fit(i);
        global_best=g(i, :);
    end
end
doc save
%%
for time = 1:num_max
    for i = 1:NP
    % 惯性操作
        position_i = randperm(n,1);   %随机生成i
        position_j = randperm(n,1);   %随机生成j
    %置换粒子i,j位置处内容
        tmp = g(i, :);
        tmp(position_i) = g(i, position_j);
        tmp(position_j) = g(i, position_i);
    % 适应值更新
        tmp_cost = fitness(tmp, dis, n, m); % 计算新的适应度
        
        % 适应度值小于个体最优值对个体最优值进行更新；
        if tmp_cost < individual_best_fit(i)
            individual_best(i, :) = tmp;
            individual_best_fit(i) = tmp_cost;
        else        % 随机数小于w，也进行粒子更新
            if rand() < w
                g(i, :) = tmp;
            end
        end
        
        % 小于种群最优值对种群最优值进行更新
        if tmp_cost < global_best_fit
            global_best = tmp;
            global_best_fit = tmp_cost;
            best = time;
        end
    end
    
    % 自我认知
    K1 = randperm(NP,1);
    if fitness(individual_best(K1, :), dis, n, m) > global_best_fit
        g(K1, :) = individual_best(K1, :);  %第i个元素的内容替换为个体最优记录的第i个元素
    end
    if rand < c1
        g(K1, :) = individual_best(K1, :);  %第i个元素的内容替换为个体最优记录的第i个元素
    end
%     if fitness(individual_best(K1, :), dis, n, m) < global_best_fit
%         A=fitness(individual_best(K1, :), dis, n, m);
%         global_best_fit = fitness(individual_best(K1, :), dis, n, m);
%     end
    
    % 社会认知
    K2 = randperm(NP,1);
    if fitness(individual_best(K2, :), dis, n, m) > global_best_fit
        g(K2, :) = global_best;  %第i个元素的内容替换为个体最优记录的第i个元素
    end
    if rand < c2
        g(K2, :) = global_best;         %粒子个体第i个元素的内容替换为种群最优记录的第i个元素
    end
%     if fitness(individual_best(K2, :), dis, n, m) < global_best_fit
%         global_best_fit = fitness(individual_best(K1, :), dis, n, m);
%     end
    gb(time) = global_best_fit;
end

fprintf("最佳迭代次数：%d 次 \nThe Minimum Distribution cost  = %f \n", best, global_best_fit);
 
[bestvalue ,route_msg] = fitness(global_best, dis, n, m);
for i = 1:m
    j = 1;
    fprintf("机器人%d：", i);
    while(j <=  10 && route_msg(i, j) > 0)
        fprintf("--> %d ", route_msg(i, j));
        j = j + 1;
    end
    fprintf("\n");
end
%%
%绘图
f1 = figure;
plot(gb);
title('收敛曲线')
xlabel('迭代次数')
ylabel('适应值')
f2 = figure;
plot(T(:,1),T(:,2),"bo",'MarkerFaceColor',[0.4660 0.6740 0.1880],'MarkerSize',10 );
L = size(route_msg);
point_setx = zeros(10);
point_sety = zeros(10);
for i = 1:L(1)
    for j = 1:10
        if route_msg(i,j)==0
            break;
        else
            point_setx(i,j+1) = T(route_msg(i,j)+1,1);
            point_sety(i,j+1) = T(route_msg(i,j)+1,2);            
        end
    end
end
hold on
plot(point_setx(1,:),point_sety(1,:),'-','LineWidth',2);
plot(point_setx(2,:),point_sety(2,:),'--','LineWidth',2);
plot(point_setx(3,:),point_sety(3,:),':','LineWidth',2);

text(T(:,1)+1.5,T(:,2)+1.5,{'起点','1','2','3','4','5','6','7','8','9','10',},'Position',[50 50],'Color',[0.8500 0.3250 0.0980],'FontSize',14);
title('路线图')
xlabel('km')
ylabel('km')
%%
%适应度函数                              %n个任务10%m个无人机4
function [value, route_msg] = fitness(L, dis, n, m)
    global T
    global demand
    global L_max
    global L_c
    global v
    global y
    judge = 0;      %用于判断是否跳出循环
    route_msg = zeros(m, 10);
    value = 0.0;
    for i = 1:m
       %暂时存储全局最短路径
       tmp_best = 10000.0;
       %找到无人机i需要完成的任务j的集合
       mission_set = find(L(1:n) == i);
       total_demand=0;
       for t =  1:length(mission_set)
            total_demand = total_demand+T(mission_set(t)+1,3);
       end
       if total_demand > demand
           value = 10^5;
           break;
       end
       %全排列，计算无人机i完成任务的所有可能路径
       per_mission_range = perms(mission_set);    %将mission_set-向量所有排列列出
       temp_arry=size(per_mission_range);         %确定排列个数
       if temp_arry(2) == 0
           continue;  
       end           
       %计算每个uav每种排列方式得代价 选择最优作为路线
       for j = 1:temp_arry(1)
           T_cost = 0;     %初始时间代价  
           %创建储存到达每一个节点的时间 的向量
           arrive_time = zeros(1,length(mission_set) + 1);
       %计算路程成本
           %初始值设为从仓库出发，机器人第一个任务的距离
           L_temp = dis(1, per_mission_range(j,1)+1);
           %机器人按照任务顺序执行任务过程中得总距离
           for k = 1:length(mission_set) - 1
               L_temp = L_temp + dis(per_mission_range(j,k)+1, per_mission_range(j,k + 1)+1);
           end
           %机器人执行完任务回到仓库得总距离
           L_temp = L_temp + dis(per_mission_range(j,length(mission_set ))+1,1) ;
           if L_temp > L_max
                judge = 1;  
                break;
           end
           L_cost=L_c*L_temp;
       %计算时间成本
           %计算到达每一个节点的时间 
           T_temp = 6 + dis(1, per_mission_range(j,1)+1)/v;
           if T_temp <  T(per_mission_range(j,1)+1,4) %如果提前到达节点，将等待节点的开始时间
               T_temp = T(per_mission_range(j,1)+1,4);
           end
           arrive_time(1) = T_temp;
           for t = 1:length(mission_set)-1
                T_temp = T_temp + dis(per_mission_range(j,t)+1, per_mission_range(j,t + 1)+1)/v+ ...,
                                   T(per_mission_range(j,t)+1,6);
               if T_temp <  T(per_mission_range(j,t+1)+1,4) %如果提前到达节点，将等待到节点的开始时间
                   T_temp = T(per_mission_range(j,t+1)+1,4);
               end
               arrive_time(t+1) = T_temp;
           end
           T_temp = T_temp + dis(per_mission_range(j,length(mission_set))+1,1)/v+ ...,
                                   T(per_mission_range(j,length(mission_set ))+1,6);
           arrive_time(length(mission_set)+1) = T_temp;
           %计算时间代价
           for c = 1:length(arrive_time)
                if c==length(arrive_time)
                    bi=T(1,5);
                else
                    bi = T( per_mission_range (j,c)+1,5);                     
                end
                tmp_num=arrive_time(c) - bi;
                T_cost=T_cost + max([tmp_num,0])*y;
           end
        %计算总成本
           total_cost = L_cost + T_cost;
        
           if total_cost < tmp_best
               %记录下无人机i最短路径及其距离
               tmp_best = total_cost;
               route_msg(i, 1: length(mission_set)) = per_mission_range(j,:);
           end
       end
       if judge == 1
           value = 10^5;
           break;
       end
       value = value + tmp_best;
    end
end