EI太阳能学报复现（粒子群算法PSO）——风-水电（抽水蓄能）联合优化运行分析

%% 清空环境 clc clear tic %记录运行开始时间 %% 基本遗传算法求解 %输入参数列表 %P_v 风电场的可利用风能 %P_g 风电场装机容量 %eta_p 水泵抽水效率 %eta_h 水力发效率 %t 每个时段长度 %C 上网电价 %C_p 抽水费用 %P_DL 舍弃的功率 %E 水库储能 %输出参数列表 %P_w 风机功率 %P_h 水力发电功率 %P_p 水泵抽水功率 %P 风电场输出功率(整个系统输出功率) %% 模型参数设置 load('P_v.mat') %输入已知参数一天内风能功率P_v n=length(P_v); %确定时间尺度 [C,C_p]=price(n); %确定上网电价 P_gmax=12; %风电场最大装机容量 P_gmin=0; P_hmax=3; P_hmin=0; P_pmax=3; P_pmin=0; E_max=24; P_min=3; P_max=8; %电网限制的功率，限制风电场输送到电网的效率 eta_p=0.8; %水泵抽水效率 eta_g=0.75; eta_h=eta_g/eta_p; %水力发电效率 t=1; %每个时段长 %% PSO算法中参数初始化 %粒子群算法中的两个参数——加速度因数，∈[0,4] c1 =2; c2 =2; w_start=0.9; w_end=0.4; %惯性因子,影响全局和局部搜索能力 Vmax=3; Vmin=-3; %速度的最大最小值 maxgen=300; %进化次数 sizepop=500; %种群规模,种群粒子数，每个粒子的维数为24×3=72个 NVAR=24; %变量个数，即维度Dim %% 产生初始粒子 %-----------------flag指示此时为发电或抽水状态-----------------------------% % flag=1 表示此时为水力发电状态 即Ph>0 % flag=0 表示此时为抽水蓄能状态 即Pp>0 %{ flag=2*ones(NVAR,1); for i=1:NVAR if P_v(i)>P_max flag(i)=0; else if P_v(i)<P_min flag(i)=1; end end end %} %{ for i=1:NVAR if flag(i)==0 P_h(:,i)=0; end end %} %-------------确定水泵抽水功率P_p的取值范围-------------------------------% Ppmin=P_hmin*ones(1,NVAR); %[0,3] Ppmin=[3 3 3 3 3 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; %Ppmin=[2 2 2 2 2 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; Ppmax=P_hmax*ones(1,NVAR); Ppmax=[3 3 3 3 3 0 3 0 0 0 0 0 0 0 2 0 0 0 2 0 0 0 0 0]; %Ppmax=[2 2 2 2 2 0 2 0 0 0 0 0 0 0 2 0 0 0 2 0 0 0 0 0]; RangePp=ones(sizepop,1)*(Ppmax-Ppmin); %矩阵取值范围 P_p=rand(sizepop,NVAR).*RangePp+ones(sizepop,1)*Ppmin; %初始化粒子群 %{ for i=1:NVAR if flag(i)==1 P_p(:,i)=0; end end %} VStepPp=rand(sizepop,NVAR)*(Vmax-Vmin)+Vmin; %初始化速度 %-----------------确定风机功率P_w的取值范围-------------------------------% Pwmax=min(P_max,P_v)'; %[0,min(8,Pv)] Pwmin=P_gmin*ones(1,NVAR); RangePw=ones(sizepop,1)*(Pwmax-Pwmin); %矩阵取值范围 P_w=rand(sizepop,NVAR).*RangePw+ones(sizepop,1)*Pwmin; %初始化粒子群 VStepPw=rand(sizepop,NVAR)*(Vmax-Vmin)+Vmin; %初始化速度 %-------------确定水力发电功率P_h的取值范围-------------------------------% Phmax=P_hmax*ones(1,NVAR); %[0,3] %Phmax=[0 0 0 0 0 1 0 3 3 0 2 3 0 3 0 2 0 0 0 2 0 0 0 0]; Phmax=[0 0 0 0 0 1 0 2 2 0 2 2 0 2 0 2 0 0 0 2 0 0 0 0]; Phmin=P_hmin*ones(1,NVAR); %Phmin=[0 0 0 0 0 0 0 3 3 0 0 3 0 3 0 0 0 0 0 0 0 0 0 0]; Phmin=[0 0 0 0 0 0 0 2 2 0 0 2 0 2 0 0 0 0 0 0 0 0 0 0]; RangePh=ones(sizepop,1)*(Phmax-Phmin); %矩阵取值范围 P_h=rand(sizepop,NVAR).*RangePh+ones(sizepop,1)*Phmin; %初始化粒子群 %{ for i=1:NVAR if flag(i)==0 P_h(:,i)=0; end end %} VStepPh=rand(sizepop,NVAR)*(Vmax-Vmin)+Vmin; %初始化速度 Vary=[P_p P_w P_h]; V=[VStepPp VStepPw VStepPh]; popmax=[Ppmax Pwmax Phmax]; popmin=[Ppmin Pwmin Phmin]; %% 优化 F=fun(P_w,P_h,P_p,sizepop,NVAR,C,C_p); %计算目标函数（适应度值） E=zeros(sizepop,NVAR+1); %---------对不符合条件的解（粒子）加上惩罚因子----------------------------% %判断哪些解需要加入惩罚因子 P=P_h+P_w; M=1000; %惩罚因子 for j=1:sizepop % 约束条件（2） for nvar=1:NVAR if P(j,nvar)>P_max %flag=1时表示不满足约束条件 F(j)=F(j)-M; %加上惩罚因子 end if P(j,nvar)<P_min %flag=1时表示不满足约束条件 F(j)=F(j)-M; %加上惩罚因子 end % 约束条件（3） if P_p(j,nvar)+P_w(j,nvar)>P_gmax F(j)=F(j)-M; end if P_p(j,nvar)+P_w(j,nvar)<P_gmin F(j)=F(j)-M; end % 约束条件（4） if P_h(j,nvar)>min(P_hmax,E(j,nvar)*eta_h/t) F(j)=F(j)-M; end E(j,nvar+1)=E(j,nvar)+t*(eta_p*P_p(j,nvar)-P_h(j,nvar)/eta_h); %下一时刻水库储能 %E(j,nvar+1)=max(0,E(j,nvar+1)); %若E小于0，应加以约束 if P_h(j,nvar)<P_hmin F(j)=F(j)-M; end % 约束条件（6） if E(j,nvar)<0 F(j)=F(j)-M; end if E(j,nvar)>E_max F(j)=F(j)-M; end % 约束条件（5） if P_p(j,nvar)>P_pmax F(j)=F(j)-M; end if P_p(j,nvar)<P_pmin F(j)=F(j)-M; end % 约束条件（7） if P_v(nvar)-P_w(j,nvar)-P_p(j,nvar)<0 F(j)=F(j)-M; end %附加约束条件（1） %{ if P_h(j,nvar)*P_p(j,nvar)~=0 F(j)=F(j)-M*(P_h(j,nvar)*P_p(j,nvar)); end %} end end % 个体极值和群体极值（初始情况） [bestfitness,I]=max(F); %找出最大的惩罚函数：bestfitness为惩罚函数值；I为序号数 zbest=Vary(I,:); %全局最佳 gbest=Vary; %个体最佳 fitnessgbest=F; %个体最佳适应度值 fitnesszbest=bestfitness; %全局最佳适应度值 %% 迭代寻优 w=zeros(1,maxgen); for i=1:maxgen %迭代次数 w(i)=w_start*(w_start-w_end)*(maxgen-i)/maxgen; %线性递减惯性权重 for j=1:sizepop %--------------------------速度更新---------------------------------------% V(j,:) =w(i)* V(j,:) + c1*rand*(gbest(j,:) - Vary(j,:)) + c2*rand*(zbest - Vary(j,:)); V(j,find(V(j,:)>Vmax))=Vmax; V(j,find(V(j,:)<Vmin))=Vmin; %确保V值范围 %----------------------------种群更新-------------------------------------% Vary(j,:)=Vary(j,:)+V(j,:); for k=1:72 if Vary(j,k)>popmax(k) Vary(j,k)=popmax(1,k); end if Vary(j,k)<popmin(1,k) Vary(j,k)=popmin(1,k); %通过此约束，使种群更新后仍满足条件 end end %------------------------------适应度值-----------------------------------% P_p=Vary(:,1:24); %{ for n=1:NVAR if flag(n)==1 P_p(:,n)=0; end end %} P_w=Vary(:,25:48); P_h=Vary(:,49:72); P=P_h+P_w; %{ for n=1:NVAR if flag(n)==0 P_h(:,n)=0; end end %} F(j)=funsel(P_w,P_h,P_p,NVAR,C,C_p,j); %新粒子适应度值 end %---------对不符合条件的解（粒子）加上惩罚因子----------------------------% %判断哪些解需要加入惩罚因子 for j=1:sizepop % 约束条件（2） for nvar=1:NVAR if P(j,nvar)>P_max %flag=1时表示不满足约束条件 F(j)=F(j)-M; end if P(j,nvar)<P_min F(j)=F(j)-M; end % 约束条件（3） if P_p(j,nvar)+P_w(j,nvar)>P_gmax F(j)=F(j)-M; end if P_p(j,nvar)+P_w(j,nvar)<P_gmin F(j)=F(j)-M; end % 约束条件（4） if P_h(j,nvar)>min(P_hmax,E(j,nvar)*eta_h/t) F(j)=F(j)-M; end E(j,nvar+1)=E(j,nvar)+t*(eta_p*P_p(j,nvar)-P_h(j,nvar)/eta_h); %下一时刻水库储能 %E(j,nvar+1)=max(0,E(j,nvar+1)); %若E小于0，应加以约束 if P_h(j,nvar)<P_hmin F(j)=F(j)-M; end % 约束条件（6） if E(j,nvar)<0 F(j)=F(