一、简介

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算法的流程和伪代码

二、源代码

clear
clc
close all
%% 参数初始化
c1=2.05;
c2=2.05;
maxgen=5000;
sizepop=30;
k=0.6;
% wV=1.1;
wP=1.1;
v=5;
popmax=30;
popmin=-30;

%      pso_option = struct('c1',1.5,'c2',1.7,'maxgen',200,'sizepop',20, ...
%                    'k',0.6,'wV',1,'wP',1,'v',5, ...
%                     'popcmax',100,'popcmin',0.1,'popgmax',10^3,'popgmin',10^(-2));
% c1:初始为1.5,pso参数局部搜索能力
% c2:初始为1.7,pso参数全局搜索能力
% maxgen:初始为200,最大进化数量
% sizepop:初始为20,种群最大数量
% k:初始为0.6(k belongs to [0.1,1.0]),速率和x的关系(V = kX)
% wV:初始为1(wV best belongs to [0.8,1.2]),速率更新公式中速度前面的弹性系数
% wP:初始为1,种群更新公式中速度前面的弹性系数
% v:初始为5,SVM Cross Validation参数
% popcmax:初始为100,SVM 参数c的变化的最大值.
% popcmin:初始为0.1,SVM 参数c的变化的最小值.
% popgmax:初始为1000,SVM 参数g的变化的最大值.
% popgmin:初始为0.01,SVM 参数c的变化的最小值.
D=10;  %%%维数
Vmax =k*popmax;
Vmin = -Vmax ;

eps =1E-5;

%% 产生初始粒子和速度
pop=zeros(sizepop,D);
V=zeros(sizepop,D);
fitness=zeros(sizepop,1);
for i=1:sizepop
    
    % 随机产生种群和速度
    pop(i,:) = (popmax-popmin)*rand(1,D)+popmin;
    
    V(i,:)=Vmax*rands(1,D);
    
    % 计算初始适应度
    fitness(i)=myfunc_fit1(pop(i,:));
    
end
Xd_ave0=repmat(sum(pop)/sizepop,sizepop,1);
D_t0=sum((sum((pop-Xd_ave0).^2,2)).^0.5)/sizepop/(popmax-popmin);
wV=1/(1+exp(-12*(D_t0-0.5)));

D_min=D_t0*0.2;

% 找极值和极值点
[global_fitness bestindex]=min(fitness); % 全局极值
local_fitness=fitness;   % 个体极值初始化

global_x=pop(bestindex,:);   % 全局极值点
local_x=pop;    % 个体极值点初始化

% 每一代种群的平均适应度
avgfitness_gen = zeros(maxgen,1);
fit_gen=zeros(maxgen,1);

%% 迭代寻优
for i=1:maxgen
    
    for j=1:sizepop
        
        %速度更新
        V(j,:) = wV*V(j,:) + c1*rand*(local_x(j,:) - pop(j,:)) + c2*rand*(global_x - pop(j,:));
        if find(V(j,:) > Vmax)
            V_maxflag=find(V(j,:) > Vmax);
            V(j,V_maxflag) = Vmax;
        end
        if find(V(j,1) < Vmin)
            V_minflag=find(V(j,1) < Vmin);
            V(j,V_minflag) = Vmin;
        end
        
        %种群更新
        pop(j,:)=pop(j,:) + wP*V(j,:);
        if find(pop(j,:) > popmax)
            pop_maxflag=find(pop(j,:) > popmax);
            pop(j,pop_maxflag) = popmax;
        end
        if find(pop(j,:) < popmin)
            pop_minflag=find(pop(j,:) < popmin);
            pop(j,pop_minflag) = popmin;
        end
        
        % 自适应粒子变异
        if rand>0.5
            k=ceil(2*rand);
            pop(j,k) = (popmax-popmin)*rand + popmin;
        end
        
        %适应度值
        fitness(j)=myfunc_fit1(pop(j,:));
        
        %个体最优更新
        if fitness(j) < local_fitness(j)
            local_x(j,:) = pop(j,:);
            local_fitness(j) = fitness(j);
        end
        
%         if fitness(j) == local_fitness(j) && length(pop(j,:) < local_x(j,:))
%             local_flag=find(pop(j,:) < local_x(j,:));
%             local_x(j,local_flag) = pop(j,local_flag);
%             local_fitness(j) = fitness(j);
%         end
        
        %群体最优更新
        if fitness(j) < global_fitness
            global_x = pop(j,:);
            global_fitness = fitness(j);
        end
        function C=PSO_FUNC(X)
global G_AC T_C hour_num Wind_V
C_W=110;   %%%风力发电
P_W=100;
u_PW=6;
m_WG=20;
r0=0.06;
v_ci=3; %切入风速
v_r=12; %额定风速
v_co=25; %切除风速
P_r=P_W;
V_t=repmat(Wind_V(:,3),52,1);

C_S=0.7;   %%%%光伏发电
P_S=0.2;
u_PS=0.009;
m_PV=25;

C_B=0.5;   %%%蓄电池
u_WB=0.0014;
m_B=5;

sigam_bat=1e-4;   %%%自放电率
N_B=2000;
W_B=0.64;
Wbat_0=0.5*N_B*W_B;
Pbat_max=0.2*N_B*W_B;
Pbat_min=-0.2*N_B*W_B;
Pbat_maxt=Pbat_max;
Pbat_mint=Pbat_min;
Wbat_t=zeros(hour_num,1);
Wbat_t(1)=Wbat_0; %%%剩余电量

Pbat_t=zeros(hour_num,1);
Pbat_t(1)=Pbat_max;


C_d=10;     %%%柴油机发电
u_Pd=0.95;
P=4.62;
Q_d0=0.22;
m_die=10;

P_STC=0.2;
G_STC=1;
K=-0.47;
Tr=298.15;
T_C=T_C+273.15;

P1_t=300;     %%%居民负荷峰值  KW
Pdes_t=200;   %%%海水淡化负荷   KW
P_des=25;      %%%单台海水淡化机组的额定功率   KW
N_des=8;     %%%海水淡化机组总台数
G_des=100/24;    %%%单台机组每小时的产水量  100t/d
Rwater_t=500/24;  %%%岛上全天用水需求   t
Rdes_min=0;
Rdes_max=8*100/24;
eta_c=0.97;

Rdes_t=zeros(hour_num,1);
Rdes_t(1)=Rdes_max;   %%最初蓄水量 

P_PV=zeros(hour_num,1);
P_WG=zeros(hour_num,1);
P_PVM_t=zeros(hour_num,1);
P_WGM_t=zeros(hour_num,1);
P_net_t=zeros(hour_num,1);
Pdie_t=zeros(hour_num,1);
C_f=0;   %%%柴油年成本
yeushu1=0;
yeushu2=0;
yeushu3=0;
yeushu4=0;
yeushu5=0;
yeushu6=0;
yeushu7=0;
yeushu8=0;
for i=1:hour_num
    P_PV(i)=P_STC*G_AC(i).*(1+K*(T_C(i)-Tr))/G_STC;
    a=P_r/(v_r^3-v_ci^3);
    b=v_ci^3/(v_r^3-v_ci^3);
    if (V_t(i)<v_ci)
        P_WG(i)=0;
    elseif (v_ci<V_t(i)<v_r)
        P_WG(i)=a*V_t(i)^3-b*P_r;
    elseif (v_r<V_t(i)<v_co)
        P_WG(i)=P_r;
    else
        P_WG(i)=0;
    end
    
    if(Rdes_t(i)-Rdes_max>=Rwater_t)
        Ndes_mint=0;
    else
       Ndes_mint=(Rwater_t-(Rdes_t(i)-Rdes_max))/G_des;
    end
    if(Rdes_t(i)+N_des*P_des-Rwater_t<=Rdes_max)
        Ndes_maxt=N_des*P_des;
    else
        Ndes_maxt=(Rdes_max+Rwater_t-Rdes_t(i))/P_des;
    end
    Pdes_mint=Ndes_mint*P_des;
    Pdes_maxt=Ndes_maxt*P_des;
    
    P_PVM_t(i)=P_PV(i);
    P_WGM_t(i)=P_WG(i);
    P_net_t(i)=P_WGM_t(i)+P_PVM_t(i)-P1_t;
    
   if(i>=2)
           Pbat_maxt=min([Pbat_max ((1-sigam_bat)*Wbat_t(i)-((1-sigam_bat)*Wbat_t(i)+Pbat_t(i)))*eta_c]);
           Pbat_mint=max([Pbat_min ((1-sigam_bat)*Wbat_t(i)-((1-sigam_bat)*Wbat_t(i)+Pbat_t(i)))/eta_c]);
   end   

三、运行结果

四、备注

完整代码或者代写添加QQ 1564658423