粒子群算法求解多物流中心选址问题（p-中值问题）

多物流中心选址问题（p-中值问题）介绍

p-中值问题模型是指在给定m个备选物流中心位置中，选取p个位置，并为每个需求点指定相应的物流中心，使得各个物流中心到需求点的运输费用之和最低(如图所示)。

基本假设

① 各需求区域和备选物流中心均当作几何上的点看待，各需求点的需求量和位置以及备选物流中心的位置已知； 
② 需求点到物流中心的距离按直线距离计算； 
③ 需求点到物流中心的运费与两点之间的距离成正比； 
④ 各需求点的需求均可得到物流中心的满足。

模型建立

中值问题这是解决这样的一个问题：假设有N个设施候选点和M个需求点，从N中选择P个设施点，使得所有M到最近设施的距离之和最小。题目所问的p-中值问题即是如此。

粒子群算法

粒子群算法（ Particle Swarm Optimization, PSO）最早是由Eberhart和Kennedy于1995年提出，它的基本概念源于对鸟群觅食行为的研究。PSO算法就从鸟类这种生物种群的觅食行为特性中得到启发并用于求解优化问题。

用一种粒子来模拟上述的鸟类个体，每个粒子可视为N维搜索空间中的一个搜索个体，粒子的当前位置即为对应优化问题的一个候选解，粒子的飞行过程即为该个体的搜索过程．粒子的飞行速度可根据粒子历史最优位置和种群历史最优位置进行动态调整．粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子单独搜寻的最优解叫做个体极值，粒子群中最优的个体极值作为当前全局最优解。不断迭代，更新速度和位置。最终得到满足终止条件的最优解。

参数设置

%% Problem Definition

model=SelectModel();                        % Select Model

CostFunction=@(xhat) MyCost(xhat,model);    % Cost Function

VarSize=[model.N model.N];      % Decision Variables Matrix Size

nVar=prod(VarSize);             % Number of Decision Variables

VarMin=0;          % Lower Bound of Decision Variables
VarMax=1;          % Upper Bound of Decision Variables


%% PSO Parameters

MaxIt=250;      % Maximum Number of Iterations

nPop=150;        % Population Size (Swarm Size)

w=1;                % Inertia Weight
wdamp=0.99;         % Inertia Weight Damping Ratio
c1=1.5;             % Personal Learning Coefficient
c2=2.0;             % Global Learning Coefficient

粒子群代码

%% Initialization

empty_particle.Position=[];
empty_particle.Cost=[];
empty_particle.Sol=[];
empty_particle.Velocity=[];
empty_particle.Best.Position=[];
empty_particle.Best.Cost=[];
empty_particle.Best.Sol=[];

particle=repmat(empty_particle,nPop,1);

BestSol.Cost=inf;

for i=1:nPop
    
    % Initialize Position
    particle(i).Position=unifrnd(VarMin,VarMax,VarSize);
    
    % Initialize Velocity
    particle(i).Velocity=zeros(VarSize);
    
    % Evaluation
    [particle(i).Cost, particle(i).Sol]=CostFunction(particle(i).Position);
    
    % Update Personal Best
    particle(i).Best.Position=particle(i).Position;
    particle(i).Best.Cost=particle(i).Cost;
    particle(i).Best.Sol=particle(i).Sol;
    
    % Update Global Best
    if particle(i).Best.Cost<BestSol.Cost
        
        BestSol=particle(i).Best;
        
    end
    
end

BestCost=zeros(MaxIt,1);


%% PSO Main Loop

for it=1:MaxIt
    
    for i=1:nPop
        
        % Update Velocity
        particle(i).Velocity = w*particle(i).Velocity ...
            +c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ...
            +c2*rand(VarSize).*(BestSol.Position-particle(i).Position);
        
        % Apply Velocity Limits
        particle(i).Velocity = max(particle(i).Velocity,VelMin);
        particle(i).Velocity = min(particle(i).Velocity,VelMax);
        
        % Update Position
        particle(i).Position = particle(i).Position + particle(i).Velocity;
        
        % Velocity Mirror Effect
        IsOutside=(particle(i).Position<VarMin | particle(i).Position>VarMax);
        particle(i).Velocity(IsOutside)=-particle(i).Velocity(IsOutside);
        
        % Apply Position Limits
        particle(i).Position = max(particle(i).Position,VarMin);
        particle(i).Position = min(particle(i).Position,VarMax);
        
        % Evaluation
        [particle(i).Cost, particle(i).Sol] = CostFunction(particle(i).Position);
        
        % Mutation
        NewParticle=particle(i);
        NewParticle.Position=Mutate(particle(i).Position, model);
        [NewParticle.Cost, NewParticle.Sol]=CostFunction(NewParticle.Position);
        if NewParticle.Cost<=particle(i).Cost || rand < 0.1
            particle(i)=NewParticle;
        end
        
        % Update Personal Best
        if particle(i).Cost<particle(i).Best.Cost
            
            particle(i).Best.Position=particle(i).Position;
            particle(i).Best.Cost=particle(i).Cost;
            particle(i).Best.Sol=particle(i).Sol;
            
            % Update Global Best
            if particle(i).Best.Cost<BestSol.Cost
                
                BestSol=particle(i).Best;
                
            end
            
        end
        
    end
    
    % Local Search based on Mutation
    for k=1:3
        NewParticle=BestSol;
        NewParticle.Position=Mutate(BestSol.Position, model);
        [NewParticle.Cost, NewParticle.Sol]=CostFunction(NewParticle.Position);
        if NewParticle.Cost<=BestSol.Cost
            BestSol=NewParticle;
        end
    end

    BestCost(it)=BestSol.Cost;
    
    disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
    
    w=w*wdamp;
    
    % Plot Best Solution
    figure(1);
    PlotSolution(BestSol.Sol,model);
    pause(0.01);
    
end

%% Results

figure;
plot(BestCost,'LineWidth',2);
xlabel('Iteration');
ylabel('Best Cost');

结果展示

欢迎加入群智能讨论群

QQ群：688369315