一、简介

1 粒子群算法的概念
粒子群优化算法(PSO：Particle swarm optimization) 是一种进化计算技术（evolutionary computation）。源于对鸟群捕食的行为研究。粒子群优化算法的基本思想：是通过群体中个体之间的协作和信息共享来寻找最优解．
PSO的优势：在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。

2 粒子群算法分析
2.1基本思想
粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子在搜索空间中单独的搜寻最优解，并将其记为当前个体极值，并将个体极值与整个粒子群里的其他粒子共享，找到最优的那个个体极值作为整个粒子群的当前全局最优解，粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。下面的动图很形象地展示了PSO算法的过程：

2 更新规则
PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest，gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。

公式(1)的第一部分称为【记忆项】，表示上次速度大小和方向的影响；公式(1)的第二部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；公式(1)的第三部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。以上面两个公式为基础，形成了PSO的标准形式。

公式(2)和 公式(3)被视为标准PSO算法。
3 PSO算法的流程和伪代码

二、源代码

%%  optimization   +  PSO 求解   选址_路径优化问题
%% 时间：2021年4月13日:12:09
clc;clear;close all;
%% 初始化种群
N = 300;                         % 初始种群个数
d = 6;                           % 空间维数
ger =60;                         % 最大迭代次数     
limit =zeros(2,2);
limit(:,1)=[-46000  63000 ];  
limit(:,2)=[-26000  70000 ] ;      
% 设置位置参数限制
vlimit(1) = -10000;               % 设置速度限制
vlimit(2)=10000 ;
wmax = 0.8;
wmin = 0.4;                        % 惯性权重
c1 = 1;                            % 自我学习因子
c2 = 2;                            % 群体学习因子 

model=CreateModel();               %调用模型函数
for i = 1:d
  if mod(i,2)==0
      x(:,i) = limit(2,1) + (limit(2,2) - limit(2,1)) * rand(N, 1);%初始种群的位置
  else
      x(:,i) = limit(1,1) + (limit(1,2) - limit(1,1)) * rand(N, 1); 
  end
end
v = rand(N, d);                  % 初始种群的速度
xm = x;                          % 每个个体的历史最佳位置
ym = zeros(1, d);                % 种群的历史最佳位置
fxm = 1./zeros(N,1);             % 每个个体的历史最佳适应度
fym = inf;                       % 种群历史最佳适应度   inf无穷大
%% 群体更新
iter = 1;
record = zeros(ger, 2);          % 记录器
average=record;
while iter <= ger
    fx=zeros(N,1);
    for i=1:N 
    fx(i)= TourLength(model,x(i,:)) ; % 个体当前适应度   
    end
    intermediate_x=zeros(size(xm));
    intermediate_x(1:N,:) = xm;
    intermediate_x(N + 1 : N + N,1 : d) =  x;

  for i=1:N*2
     intermediate_x(i,d+3)=150;
     intermediate_x(i,d+1:d+2)=Tour(model, intermediate_x(i,:)) ; % 个体当前适应度   
  end
    intermediate_x = ...
        non_domination_sort_mod(intermediate_x, 2, d);
    intermediate_x = replace_x(intermediate_x, 2, d, N);

    record(iter,:) = intermediate_x(1,d+1:d+2);%最小值记
    average(iter,1)=sum(intermediate_x(:,d+1))/N;
    average(iter,2)=sum(intermediate_x(:,d+2))/N;
    ym=intermediate_x(1,1:d);
    fym=record(iter,:);
    
    disp(['第',num2str(iter),'次迭代''最小值：',num2str(record(iter,:)),'抢修中心坐标：',num2str(ym)]);
    iter = iter+1;
    xm=intermediate_x(:,1:d);
    
    w=wmax-(wmax-wmin)*(ger-iter)/ger  ;%权重更新
    v = v * w + c1 * rand * (xm - x) + c2 * rand * (repmat(ym, N, 1) - x);% 速度更新
    % 边界速度处理
    for ii=1:N
        for jj=1:d       
     if v(ii,jj)>vlimit(2)  v(ii,jj)= vlimit(2);end
      if v(ii,jj)<vlimit(1)  v(ii,jj)= vlimit(1);end
        end
    end
    x = x + v;% 位置更新
    % 边界位置处理
        for ii=1:N
        for jj=1:d         
             if mod(jj,2)==0 
                 if x(ii,jj)>limit(2,2)  x(ii,jj)= limit(2,2);end
                   if x(ii,jj)<limit(2,1)  x(ii,jj)= limit(2,1);end
             else
             if  x(ii,jj)>limit(1,2)  x(ii,jj)= limit(1,2);end
              if x(ii,jj)<limit(1,1)  x(ii,jj)= limit(1,1);end
             end
        end
        end
end

Schedule = code(x, model);  %调用
%% 绘制测试结果图
% 故障点
    figure(1);
    x=model.trouble(1,:);
    y=model.trouble(2,:);
    for i=1:39
        xx=x(i);
        yy=y(i);
        [ch1]=plot( xx,yy,'ks',...
           'MarkerFaceColor','k',...
           'MarkerSize',4); hold on
       text(xx , yy, num2str(i) );   %带箭头的标注
       hold on
    end
% 绘制抢修中心
for i=1:3
    xx=ym(i*2-1);
    yy=ym(i*2);
    [ch2]= plot( xx ,yy ,'ko',...
        'MarkerFaceColor','k',...
        'LineWidth',6);hold on
    text(xx , yy, num2str( i) );   %带箭头的标注
    hold on
end

 %% 规划结果
 rand('seed',  0);
 C= unifrnd(  0.1, 0.2,  model.Num_CenterSletectd , 3) ;%unifrnd可以创建随机的连续均匀分布的数组
for i = 1: model.Num_CenterSletectd
     Center = Schedule(i).Center;
     Client =  Schedule(i).Client ;
     xx=model.trouble(1,:);
     yy=model.trouble(2,:);
    for j= Client
        x = [ ym(i*2-1) , xx(j ) ] ;
        y = [ ym(i*2)   , yy(j ) ] ;
 plot(   x ,y , '-' ,'color' ,C(i, : ) , 'linewidth' , 1.5 );  
    hold on
    end
end

% 标注
legend([ch1, ch2], {'故障点'  ,'供电所位置' }); %  'Location','SouthWestOutside'注释放置位置
xlabel('X/m','fontsize',15,'fontname','Times new roman');
ylabel('Y/m','fontsize',15,'fontname','Times new roman');
title('抢修分布图'); 
axis on
set(gcf,'color',[1 1 1]);  %设置背景颜色
%% 绘制供电半径图
figure(2);
plot(model.trouble(1, : ), model.trouble(2, : ),'ks',...
      'MarkerFaceColor','k',...
      'MarkerSize',3);hold on
for i=1:3
    plot(ym(i*2-1),ym(i*2),'ko',...
        'MarkerFaceColor','k',...
        'MarkerSize',6);hold on    %  抢修中心分布图 ，半径
    x=ym(i*2-1);
    y=ym(i*2);
    r=model.point(3,i)*1000;
    theta=0:2*pi/3600:2*pi;
    Circle1=x+r*cos(theta);
    Circle2=y+r*sin(theta);
    plot(Circle1,Circle2,'k-','Linewidth',1);
    axis equal
end
legend('故障点'  ,'供电所位置', '供电半径');
title('抢修分布图');
xlabel('X/m','fontsize',15,'fontname','Times new roman');
ylabel('Y/m','fontsize',15,'fontname','Times new roman');

   for i=1:N 
      xm(i,d+1:d+2)=Tour(model,xm(i,1:d)) ; % 个体当前适应度   
   end
   xm = non_domination_sort_mod(xm, 2, d);
 function f  = replace_chromosome(intermediate_chromosome, M, V,pop)


[N, m] = size(intermediate_chromosome);

% Get the index for the population sort based on the rank
[temp,index] = sort(intermediate_chromosome(:,M + V + 1));

clear temp m

% Now sort the individuals based on the index
for i = 1 : N
    sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);
end

% Find the maximum rank in the current population
max_rank = max(intermediate_chromosome(:,M + V + 1));

% Start adding each front based on rank and crowing distance until the
% whole population is filled.
previous_index = 0;
for i = 1 : max_rank
    % Get the index for current rank i.e the last the last element in the
    % sorted_chromosome with rank i. 
    current_index = max(find(sorted_chromosome(:,M + V + 1) == i));
    % Check to see if the population is filled if all the individuals with
    % rank i is added to the population. 
    if current_index > pop
        % If so then find the number of individuals with in with current
        % rank i.
        remaining = pop - previous_index;
        % Get information about the individuals in the current rank i.
        temp_pop = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        % Sort the individuals with rank i in the descending order based on
        % the crowding distance.
        [temp_sort,temp_sort_index] = ...
            sort(temp_pop(:, M + V + 2),'descend');
        % Start filling individuals into the population in descending order
        % until the population is filled.
        for j = 1 : remaining
            f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);
        end
        return;
    elseif current_index < pop
        % Add all the individuals with rank i into the population.
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
    else
        % Add all the individuals with rank i into the population.
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        return;
    end

三、运行结果

四、备注

版本：2014a