MATLAB风光柴储微网粒子群算法

本程序实现了风光柴储微网中的粒子群优化（PSO）算法，用于优化微网的能源调度问题。具体来说，程序考虑了光伏发电、风力发电、柴油机发电（柴储），并使用粒子群算法来优化这些能源的调度，以满足负载需求。以下是对代码的详细说明：

1. PSO算法参数设置

最大迭代次数 maxgen：设定粒子群算法的最大迭代次数为2000次。
种群规模 N：设定粒子群中的粒子数量为100。
粒子速度限制 v_max, v_min：粒子速度的上限和下限，分别为2和-2。
惯性权重 w_max, w_min：控制粒子群在搜索空间中探索的幅度。
学习因子 c1, c2：分别表示个体认知学习因子和社会学习因子。

2. 初始化

负载 P_load：表示微网的负荷需求。
光伏发电 solar：光伏电池板的发电量，按时间分布。
风力发电 pv：风力发电机组的发电量，按时间分布。
发电总量 P_pv：光伏和风力发电总和。

粒子群的初始化过程中，粒子位置 x 和速度 v 被随机生成。粒子的位置对应了微网中不同能源的发电量（包括风力、光伏、柴油发电机组等），这些位置的初始值是根据各能源的最大发电能力进行约束的。

3. 适应度计算

每个粒子的适应度是通过调用 fitness 函数来计算的。适应度函数衡量了粒子方案的优劣，旨在通过最小化成本、最大化系统效率或实现其他优化目标。

4. 主循环

在每次迭代中，粒子的位置和速度会根据以下公式进行更新：

速度更新：结合了惯性项、个体认知项和社会认知项。
位置更新：通过粒子的当前速度更新位置。

每个粒子的边界都会根据各个能源的最大值和最小值进行限制，并采用反转速度的方式进行边界处理，以确保粒子不会越过边界。

5. 粒子群更新

每个粒子在更新其位置和速度后，会重新评估其适应度。如果当前粒子的适应度比之前的最优适应度更好，则更新该粒子的最佳位置。全局最优解（即整个群体中最优的粒子）会在每次迭代后更新。

6. 结果可视化

程序会绘制多个图表展示优化结果：

图1：展示已消纳的功率与总光伏发电功率的关系。
图2：展示已消纳的功率与总风电发电功率的关系。
图3：展示光伏发电的已消纳功率。
图4：展示风电发电的已消纳功率。

这些图表帮助分析不同能源的调度效果，以及微网在负荷需求下的能源使用情况。

7 主程序代码

clc;
clear;
%rng default;
%pso参数设置
format long;
maxgen=2000;%迭代次数
N=100;%种群规模
%c1=2;%自我学习因子
%c2=2;%群体学习因子
v_max=2; v_min=-2;%个体速度
w_max=0.9;w_min=0.1;
dmaxp=10;
dminp=0;
wminp=-10;
wmaxp=10;
s=1;
P_load=[52 50 50 51 56 63 70 75 76 80 78 74 72 72 76 80 85 88 90 87 78 71 65 56 ];%负载
%P_load=[152 150 150 151 156 163 170 175 176 180 178 174 172 172 176 180 185 188 190 187 178 171 165 156 ];
solar=[0 0 0 0 1 3 5 10 15 16 17 20 20 20 18 15 15 10 10 5 0 0 0 0 ]%光伏
pv=[75 65 60 55 50 57 60 66 50 58 60 60 57 50 61 58 55 65 56 60 72 66 75 76 ];%风力
P_pv=pv+solar;
%初始化种群
for i=1:N;for j=1:96;%4个维度v(i,j)=0.0;if j<=24;x(i,j)=rand()*pv(j);elseif j>24&&j<=48;x(i,j)=rand()*solar(j-24);elseif    j>48&&j<=72;x(i,j)=dminp+rand()*(dmaxp-dminp);elseif   j>72&&j<=96;  x(i,j)=wminp+rand()*(wmaxp-wminp);endend
end%计算各个粒子的适应度，并初始化Pi和Pg
for i=1:Np(i)=fitness(x(i,:),s);y(i,:)=x(i,:);%每个粒子的个体寻优值
end
Pbest=fitness(x(1,:),s);
pg=x(1,:);%Pg为全局最优
for i=2:Nif fitness(x(i,:),s)<fitness(pg,s)Pbest=fitness(x(i,:),s);pg=x(i,:);%全局最优更新end
end
%进入主循环
for t=1:maxgenfor i=1:Nw=w_max-(w_max-w_min)*t/maxgen;%惯性权重更新c1=(0.5-2.5)*t/maxgen+2.5; %认知c2=(2.5-0.5)*t/maxgen+0.5; %社会认识 
v(i,:)=w*v(i,:)+c1*rand()*(y(i,:)-x(i,:))+c2*rand()*(pg-x(i,:));
for m=1:96if(v(i,m)>v_max)v(i,m)=v_max;elseif(v(i,m)<-v_max)v(i,m)=-v_max;end
endx(i,:)=x(i,:)+v(i,:);%对粒子边界处理for n=1:96if n<25if x(i,n)<0x(i,n)=0;v(i,n)=-v(i,n); elseif x(i,n)>pv(n)x(i,n)=pv(n);v(i,n)=-v(i,n); %反转运动方向endelseif n>24&&n<49if x(i,n)<0x(i,n)=0;v(i,n)=-v(i,n);  elseif  x(i,n)>solar(n-24)x(i,n)=solar(n-24);v(i,n)=-v(i,n);                     endelseif n>49&&n<73if x(i,n)<dminpx(i,n)=dminp;v(i,n)=-v(i,n); elseif x(i,n)>dmaxpx(i,n)=dmaxp;v(i,n)=-v(i,n); endelse if x(i,n)<wminpx(i,n)=wminp;v(i,n)=-v(i,n); elseif x(i,n)>wmaxpx(i,n)=wmaxp;v(i,n)=-v(i,n); endendendend %对粒子进行评价，寻找最优值if fitness(x(i,:),t)<p(i)p(i)=fitness(x(i,:),t);y(i,:)=x(i,:);end if p(i)<PbestPbest=p(i);pg=y(i,:);s=t;end
endfor m=1:24pg1(m)=pg(m);
end
for m=25:48pg2(m-24)=pg(m);
end
for m=49:72pg3(m-48)=pg(m);
end
for m=73:96pg4(m-72)=pg(m);
end  figure(1);
plot( pg1,'r-')
hold on
plot( pv,'-')
xlim([1 24])
grid 
legend('已消纳功率','总光伏功率');
xlabel('时间');
ylabel('功率');
title('风机发电出力')figure(2);
plot( pg2,'g-')
hold on
plot( solar,'-')
xlim([1 24])
grid
legend('已消纳功率','总风电功率');
xlabel('时间');
ylabel('功率');
title('光伏发电出力')figure(3);
plot( pg3,'b-')
xlim([1 24])
grid       figure(4);
plot( pg4,'k-')
xlim([1 24])
grid   %plot (P_load,'-')
grid on