粒子群算法应用于重心法选址问题-基于MATLAB实现

朱佩棋（代码版）

已于 2022-11-20 14:55:37 修改

阅读量3.1k

点赞数 7

分类专栏：启发式算法文章标签：启发式算法算法 matlab

于 2022-10-05 16:44:54 首次发布

本文链接：https://blog.csdn.net/m0_73970831/article/details/127174071

版权

启发式算法专栏收录该内容

7 篇文章 2 订阅

订阅专栏

本文是个人学习过程中记录一下笔记，借鉴了各位大佬的内容，仅供学习使用，如果侵一定删！

1、根据求平面中物体重心的方法求配送中心坐标点（模拟重心法）

2、精确重心法（基于MATLAB实现迭代过程）

3、粒子群算法（Particle Swarm Optimization）的基本原理

选址问题：设在某计划区域内有n个需求点，它们各自的坐标为 $(x_{i},y_{_{i}})$ ，各点的需求量为 $w_{i}$ $_{}$ ，配送中心到需求点的运费率是 $r_{i}$ ，设该配送中心的坐标是（x,y）。

1、根据求平面中物体重心的方法求配送中心坐标点（模拟重心法）

代入数值，实际求得的值，即为所求得配送中心位置的坐标。

2、精确重心法（基于MATLAB实现迭代过程）

假设配送中心到各需求点的运费为 $r_{i}$ ，总的运费用为 D，则有：

要想使求得总费用D最小的方法是根据一阶偏导数为零的原理求解的，计算公式如下：

如果从两个方程式的右边完全消除 $x_{0}$ 和 $y_{0}$ ，计算将变得很复杂，计算量也很大。因此，可以采用迭代的方法进行计算，迭代过程如下：

可用模拟重心公式计算得出的结果作为初始解代入
利用方程式（1-4）和（2-1）计算该重心点与需求点的总运费
把分别代入方程式（1-4）和（2-2）中，计算配送中心的改善坐标

计算相对应的总的运输费用C,并与上一次总运费进行比较，如果，则说明为最优解；如果，则说明计算结果得到改善，并且有待更进一步优化，于是返回第三步进行迭代计算。

利用MATLAB对上面的迭代计算进行演示，本例的基础是共有33个需求点，选择一个配送中心使得配送费用最低。为了简便计算，将运费率和需求量均设置为1。

%%精确重心法迭代计算过程
%% 输入各需求点坐标
clc;clear all; %清屏
KnownPoint=[23.58961883	34.50001618;14.04321388	70.08361889;128.643594	56.97054047
    49.94811626	54.10749291;64.12058384	40.96510475;182.5604304	17.99404973
    135.2960124	14.13726843;170.6347982	51.45892853;84.3409894	99.43543245
    180.8850119	17.10261255;15.78952671	61.96977021;30.75658012	33.64912783
    179.6510086	66.71733882;13.76432092	44.74775586;65.34921128	26.62973917
    40.3213122	12.65279543;187.4120223	25.49512307;137.8455491	74.17467689
    85.64481108	22.74382151;121.6057362	81.17053312;68.78676629	79.27644001
    113.1183539	57.04987996;40.18263065	99.81085283;37.24769422	34.50411952
    158.8958357	57.79362356;80.95579814	52.59104834;164.6554548	69.93213962
    185.7131125	49.19363862;140.2637087	89.34146483]; %各需求点坐标
w=ones(1,size(KnownPoint,1)); %各需求点的需求量
r=ones(1,size(KnownPoint,1)); %配送中心到需求点的运费率
num=1;
%% 计算初始重心坐标及运费
[x0,y0,C0,distance]=chushi(KnownPoint,w,r); %用重心公式计算得出的结果作为初始解
disp(['初始重心点坐标为:(',num2str(x0),',',num2str(y0),')']); %输出初始重心点坐标
LsatCost=C0; %上一次计算的运费
NextCost=0; %迭代后计算的运费
%% 开始迭代，根据初始重心点计算迭代后的改善坐标和运费
while 1
    [xi,yi,C,dis_diedai]=diedai(KnownPoint,w,r,distance,num);
    disp(['第',num2str(num),'次迭代重心点坐标为:(',num2str(xi),',',num2str(yi),')']);
    NextCost=C;
    num=num+1;
    distance=dis_diedai;
    %%%判断是否满足终止条件%%%
    if (NextCost>LsatCost) %把下一次迭代费用与上次费用作比较
        break;
    end
    LsatCost=NextCost;
end
   %%%%%%%%计算初始坐标函数（模拟重心法）%%%%%%%%
%模拟重心法函数：xo=sum(ri*wi*xi)/sum(ri*wi)
%               yo=sum(ri*wi*yi)/sum(ri*wi)
%% 模拟重心公式计算初始重心坐标
function [x0,y0,C0,distance]=chushi(KnownPoint,w,r) %用重心公式计算得出的结果作为初始解代入
%%初始化参数
n=size(KnownPoint,1);
z1=zeros(n,1);
z2=zeros(n,1);
z3=zeros(n,1);
distance=zeros(n,1);
C=zeros(n,1);
%%计算初始位置坐标
for i=1:n
    z1(i)=r(i)*w(i)*KnownPoint(i,1);
    z2(i)=r(i)*w(i);
    z3(i)=r(i)*w(i)*KnownPoint(i,2);
end
x0=sum(z1)/sum(z2);
y0=sum(z3)/sum(z2);
%%计算各需求点到重心的距离和相应的总的运输费用
for i=1:n
    distance(i)=sqrt((x0-KnownPoint(i,1))^2+(y0-KnownPoint(i,2))^2);
    C(i)=r(i)*w(i)*distance(i);
    disp(['第',num2str(i),'个地点的初始距离：',num2str(distance(i))]);
end
C0=sum(C); %计算初始重心点到各需求点的总运费
disp(['初始运费：',num2str(C0)]);
end
   %%%%%%%%%%%精确重心法函数%%%%%%%%%
%目标函数x=sum(ri*wi*xi)/di除以sum(ri*wi)/di
%       y=sum(ri*wi*yi)/di除以sum(ri*wi)/di
%% 开始迭代
function [xi,yi,C,dis_diedai]=diedai(KnownPoint,w,r,distance,num)
%%初始化参数
n=size(KnownPoint,1);
z1=zeros(n,1);
z2=zeros(n,1);
z3=zeros(n,1);
dis_diedai=zeros(n,1);
TotalCost=zeros(n,1);
%%计算配送中心的改善坐标
for i=1:n
    z1(i)=r(i)*w(i)*KnownPoint(i,1)/distance(i);
    z2(i)=r(i)*w(i)/distance(i);
    z3(i)=r(i)*w(i)*KnownPoint(i,2)/distance(i);
end
xi=sum(z1)/sum(z2);%采用迭代的方法进行计算配送中心的改善地点
yi=sum(z3)/sum(z2);
%%计算各需求点到改善重心点的距离以及运费
for j=1:n
    dis_diedai(j)=sqrt((xi-KnownPoint(j,1))^2+(yi-KnownPoint(j,2))^2);
    TotalCost(j)=r(j)*w(j)*dis_diedai(j);
    disp(['第',num2str(j),'个地点第',num2str(num),'次迭代距离：',num2str(dis_diedai(j))]);
end
C=sum(TotalCost); %累加改善后配送中心到各需求点的总运费
disp(['第',num2str(num),'次迭代运费：',num2str(C)]);
end

3、粒子群算法（Particle Swarm Optimization）的基本原理

粒子群算法（particle swarm optimization，PSO）是由Kennedy和Eberhart提出的全局优化算法，基于群体智能理论的优化算法，通过群体中粒子间的合作与竞争产生的群体智能优化搜索。PSO 采取速度-位移模型，可以动态跟踪当前的搜索情况从而调整搜索策略，具有“自我学习”和“向他人学习”的优点，从而能在较少迭代次数内寻找最优解。

3.1PSO算法流程

初始化PSO参数

初始化粒子群算法的基本参数：问题维度（dimension）、种群规模（population）、迭代次数（interaction）、粒子位置Xi、粒子速度Vi 、个体历史最优位置（pbest）、群体全局最优位置（gbest）、位置限制范围（Xlimit）、速度限制范围（Vlimit）、学习因子（自我学习因子C1；社会学习因子C2）、惯性权重w

2.产生粒子群，根据以下公式初始化每个粒子的初始位置和速度
 %进入for循环，随机生成一个种群
for i=1:pop 
    for d=1:Dimension
        x(i,d)=Lb(d)+(Ub(d)-Lb(d))*rand(1);
        V(i,d)=Vmin(d)+(Vmax(d)-Vmin(d))*rand(1);
    end
end
3.开始迭代，计算每个粒子的适应度
4.更新粒子的历史最优值和全局最优值

本式以最小化问题为例：

if fit（xi）<fit（pbesti）

pesti=xi

if fit（pbesti）< fit (gbesti)

gbesti=pbesti

5.根据以下公式更新粒子的速度和位置（pso的重点公式）

速度更新: $v^{k+1}= wv^{k}+c_{1}rand()\left (pbest- x^{k} \right )+c_{2}rand()\left ( gbest- x^{k}\right )$

位置更新： $x^{k+1}=x^{k}+v^{k+1}$

5.边界处理

在迭代过程中，如果粒子i的位置超过了位置和速度的边界，就将该粒子i的位置和速度设置为边界值。if Vi > Vmax,则Vi=Vmax；if Vi<-Vmax，则Vi =- Vmax

6.判断是否达到终止条件（迭代次数），若达到最大迭代次数，则输出全局最优解

3.2 PSO的流程图和伪代码

（2）伪代码

%%本例以求最小值为目标

procedure PSO

for each particle（i）

Initialize position(Xi) and velocity(Vi) for each particle

Evaluate particle(i) and set pbest=Xi

end for

gbest=min{pbest}

while not stop

for i=1to pop

Update the position and velocity of each particle

Evaluate particle（i）

if fitness（Xi）<fitness(pbest)

pbest=Xi

if fitness(gbest)< fitness(pbest)

gbest=pbest;

end for

end while

print gbest

end procedure

3.3全局模式PSO与局部模式 PSO

全局模式（Global Version PSO）:每个粒子的运动轨迹受群体中所有粒子状态的影响，粒子追寻自身极值和种群全局极值。全局模式的收敛速度相对较快，但是

局部模式（Local Version PSO）:每个粒子只收自身的认识和临近的粒子状态影响，不是被种群中所有粒子的状态影响，只追寻自身极值和临近粒子中的局部极值。局部模式具有较高的鲁棒性，但收敛速度相对较慢。

局部模式的速度更新公式：

$v^{k+1}= wv^{k}+c_{1}rand()\left (pbesti- x^{k} \right )+c_{2}rand()\left ( nbesti- x^{k}\right )$

3.4粒子群算法的特点

粒子群算法和其他进化算法都是用种群表示一组解空间中的个体集合，采用随机初始化种群方法，使用适应度评价个体，然后进行一定的随机搜素，故不能保障一定能找到最优解；

具有一定的选择性。PSO通过不同代种群之间的竞争实现种群的进化过程，具有一定的选择机制；

算法具有并行性。搜索过程从一个解集开始的，不容易陷入局部最小值，这个并行性在计算机上易于实现，提高了算法的性能和效率。

算法收敛更快。PSO进化过程中同时具有记忆位置和速度更新，在全局模式中只有全局最优粒子提供信息给其他粒子，整个搜索过程是跟随当前最优解的过程，因此所有粒子很可能较快收敛于最优解。

4.PSO算法的MATLAB编码过程

%%     I.粒子群算法应用重心法选址问题
clear all;%清除工作区所有变量
close all;%关闭所有图形窗口
clc;%清屏
%% II.输入已知仓库地点坐标（x,y)
KnownPoint=[23.58961883	34.50001618
    14.04321388	70.08361889
    128.643594	56.97054047
    49.94811626	54.10749291
    64.12058384	40.96510475
    182.5604304	17.99404973
    135.2960124	14.13726843
    170.6347982	51.45892853
    84.3409894	99.43543245
    180.8850119	17.10261255
    15.78952671	61.96977021
    30.75658012	33.64912783
    179.6510086	66.71733882
    13.76432092	44.74775586
    65.34921128	26.62973917
    40.3213122	12.65279543
    187.4120223	25.49512307
    137.8455491	74.17467689
    85.64481108	22.74382151
    121.6057362	81.17053312
    68.78676629	79.27644001
    113.1183539	57.04987996
    40.18263065	99.81085283
    37.24769422	34.50411952
    158.8958357	57.79362356
    80.95579814	52.59104834
    164.6554548	69.93213962
    185.7131125	49.19363862
    140.2637087	89.34146483];%需求点坐标
%% III.参数初始化
c1=2; %自我认知学习因子
c2=2; %社会认知学习因子
w=0.5; %慣性权重
MaxInteration=500; %最大迭代次数
Sizepop=30; %种群规模
Dimension=2; %问题维度
Ub =[300 300]; %搜索范围的上界
Lb=[-200 -200];%搜索范围的下界
Vmax =0.009*Ub; %[100 100];%粒子速度的上界
Vmin=0.009*Lb; %[-100 -100];%粒子速度的下界
%%%%给矩阵分配内存%%%%
x=zeros(Sizepop,Dimension);%%初始化种群的位置；
V=zeros(Sizepop,Dimension);%%初始化种群的速度；
ObjectiveValue=zeros(1,Sizepop);%%初始化目标函数值；
FitnessValue=zeros(1,Sizepop);%%初始化适应度值；
PBestPos=zeros(Sizepop,Dimension);%%历史最优位置；
PBestObj=zeros(1,Sizepop);%%历史最优目标函数值；
PBestFit=zeros(1,Sizepop);%%历史最优适应度值；
GBestPos=zeros(1,Dimension);%%全局最优位置；
GBestObj=0;%%全局最优目标函数值；
GBestFit=0;%%全局最优适应度值；
Curvalue=zeros(MaxInteration);%%绘制收敛曲线；
%% IV.产生初始粒子和速度
for i=1:Sizepop  %进入for循环，随机生成一个种群
    for d=1:Dimension
        x(i,d)=Lb(d)+(Ub(d)-Lb(d))*rand(1);
        V(i,d)=Vmin(d)+(Vmax(d)-Vmin(d))*rand(1);
    end
end
%% V.开始迭代
for t=1:MaxInteration
    for i=1:Sizepop
        ObjectiveValue(1,i)=fun(KnownPoint,x(i,1),x(i,2)); %计算目标函数值；
        FitnessValue(1,i)=fitness(ObjectiveValue(1,i)); %计算适应度值；
        if(FitnessValue(1,i)>PBestFit(1,i)) %如果适应度值大于历史最优；
            PBestFit(1,i)=FitnessValue(1,i); %更新历史最优适应度值；
            PBestObj(1,i)=ObjectiveValue(1,i); %更新历史最优目标；
            PBestPos(i,:)=x(i,:); %更新历史最优位置；
        end
        if(FitnessValue(1,i)>GBestFit) %如果适应度值大于全局最优；
            GBestFit=FitnessValue(1,i); %更新全局最优适应度值；
            GBestObj=ObjectiveValue(1,i); %更新全局最优目标；
            GBestPos=x(i,:); %更新全局最优位置；
        end
    end
    for i=1:Sizepop %更新每个粒子的速度和位置
        V(i,:)=w*V(i,:)+c1*rand()*(PBestPos(i,:)-x(i,:))+c2*rand()*(GBestPos-x(i,:)); %速度更新
        for d=1:Dimension %速度越界处理
            if(V(i,d)>Vmax(d))
                V(i,d)=Vmax(d);
            end
            if(V(i,d)<Vmin(d))
                V(i,d)=Vmin(d);
            end
        end
        x(i,:)=x(i,:)+V(i,:); %位置更新
        for d=1:Dimension %位置越界处理
            if(x(i,d)>Ub(d))
                x(i,d)=Ub(d);
            end
            if(x(i,d)<Lb(d))
                x(i,d)=Lb(d);
            end
        end
    end
    Curvalue(t)=GBestObj;
end
%% VI.输出结果
disp ("最低运输费用C="+GBestObj); %输出全局最优目标函数值
disp("配送中心位置坐标:"+"("+GBestPos(1,1)+"," +GBestPos(1,Dimension)+")"); %输出配送中心位置坐标
semilogy(Curvalue); %创建绘图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%目标函数f(x)=((xi-x)^2+(yi-y)^2)^0.5
function TotalCost=fun(KnownPoint,x,y)
TotalCost=0;
j=size(KnownPoint,1);%返回矩阵KnownPoint的行数
for i=1:j %利用for循环对配送中心到各个仓库的费用进行累加求和
Distance= sqrt((KnownPoint(i,1)-x)^2+(KnownPoint(i,2)-y)^2); %配送中心到仓库的距离
TotalCost= TotalCost+Distance; %计算所选配送中心到各仓库的总费用
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 计算目标函数的适应度值
function fitnessValue=fitness(x)
fitnessValue=1./(x+0.0001);
end