需求响应前风光储能粒子群算法优化调度附matlab代码.zip

clc; clear; %% 参数区 %参数初始化 c1=3; %个体学习因子 c2=3; %社会学习因子 w_max=0.9; %设置最大惯性权重为0.9 w_min=0.4; %设置最小惯性权重为0.4 max_die_dai=250; %迭代次数设置为200 size_zhong_qun=100; %粒子种群规模设置为100 WT_N=randi(10,size_zhong_qun,1); %初始化粒子风机的位置 PV_N=randi(100,size_zhong_qun,1); %初始化粒子光伏的位置 DE_N=randi(7,size_zhong_qun,1); %初始化粒子柴油发电机的位置 SB_N=randi(150,size_zhong_qun,1)+100; %初始化粒子蓄电池的位置 X=round([WT_N PV_N DE_N SB_N]); v=3.*rands(size_zhong_qun,4); %初始化粒子的飞翔速度 %定义适应度函数 f=zeros(size_zhong_qun,1); for i=1:size_zhong_qun f(i)=fitness(X(i,1),X(i,2),X(i,3),X(i,4)); end v_max=3; %速度最大值为3 v_min=-3; %速度最小值为-3 WT_N_max=12; %种群WT中个体位置的最大值 WT_N_min=1; %种群WT中个体位置的最小值 PV_N_max=300; %种群PV中个体位置的最大值 PV_N_min=1; %种群PV中个体位置的最小值 DE_N_max=7; %种群DE中个体位置的最大值 DE_N_min=1; %种群DE中个体位置的最小值 SB_N_max=800; %种群SB中个体位置的最大值 SB_N_min=1; %种群SB中个体位置的最小值 %群体中个体的最佳位置随着迭代的进行而不断发生变化 personal_best_local_x=X; %种群中每个个体的最佳位置，随着迭代的进行，数值将会发生改变 personal_best_value_x=f; %种群中每个个体的最佳适应度值 [global_best_value_x,i]=min(f); %计算出当前的全局最佳值和最佳值所对应的位置 global_best_local_x=personal_best_local_x(i,:); %当前的全局最佳位置 %% %迭代寻优 for i = 1:max_die_dai disp(['当前迭代次数为：',num2str(i),' 剩余迭代次数：',num2str(max_die_dai-i)]); %更新惯性权重因子 w=w_max-((w_max-w_min)/max_die_dai)*i; %粒子群wt的更新 for j = 1:size_zhong_qun % 速度更新 v(j,:) = w*v(j,:) + c1*rand*(personal_best_local_x(j,:)- X(j,:)) + c2*rand*(global_best_local_x - X(j,:)); if(v(j,:)> v_max) v(j,:) = v_max; else if(v(j,:)< v_min) v(j,:) = v_min; end end % 种群更新 X(j,:) = round(X(j,:) + v(j,:)); if( X(j,1)> WT_N_max) X(j,1) = WT_N_max; else if(X(j,1)< WT_N_min) X(j,1)= WT_N_min; end end if( X(j,2)> PV_N_max) X(j,2) = PV_N_max; else if(X(j,2)< PV_N_min) X(j,2)= PV_N_min; end end if( X(j,3)> DE_N_max) X(j,3) = DE_N_max; else if(X(j,3)< DE_N_min) X(j,3)= DE_N_min; end end if( X(j,4)> SB_N_max) X(j,4) = SB_N_max; else if(X(j,4)< SB_N_min) X(j,4)= SB_N_min; end end % 适应度值更新 f_update(j) = fitness(X(j,1),X(j,2),X(j,3),X(j,4)); end for j = 1:size_zhong_qun % 个体最优更新 if (f_update(j)<personal_best_value_x(j,:)) personal_best_local_x(j,:) = X(j,:); %更新粒子群x中个体最优的位置 personal_best_value_x(j,:) = f_update(j); %更新个体的最优值 end % 群体最优更新 if (f_update(j) < global_best_value_x) global_best_local_x = X(j,:); %更新粒子群中全局最优位置 global_best_value_x = f_update(j); %更新粒子群中全局最优值 end end yy1(i,:) = global_best_value_x; %将每次迭代产生的全局最优值储存在yy1中 end %% %画图 plot(yy1) grid minor; title('目标函数f的最优适应度','fontsize',12); xlabel('迭代次数','fontsize',12);ylabel('适应度','fontsize',12); %% 以下部分是为了绘制结果图和显示中间变量 num_wt=global_best_local_x(1); %风机数量 num_pv=global_best_local_x(2); %光伏板数量 num_g=global_best_local_x(3); %柴油发电机数量 num_sb=global_best_local_x(4); %蓄电池数量 load DATE; Load=DATE(3,:);%全年负荷数据 Speed_WT=DATE(1,:);%全年风速数据 Solar_PV=DATE(2,:);%全年光照强度数据 T_c=DATE(4,:);%光伏板表面温度 %% 参数区 %1.风机 single_WT=35;%单台风机的额定功率/kW v_ci=3;%切入风速 v_co=25;%切出风速 v_n=11;%额定风速 P_wt=zeros(1,8760);%单台风机的年功率曲线 u_wt=3.86;%风机投资成本（台/万元） k_wt=0.02;%风机维护费用（万元/年） Y_wt=20;%风机的使用年限 R_wt=3;%风机单台置换费用（万元/台） %2.光伏 single_PV=1;%单个光伏板额定功率/kW T_r=25;%太阳能电池温度为25℃ namd=-0.0047;%温度条件系数 G_n=1;%光谱辐照度为1000W/m2 P_pv=zeros(1,8760);%单片光伏板的年功率曲线 u_pv=0.8;%光伏投资成本（台/万元） k_pv=0.002;%光伏维护费用（万元/年） Y_pv=20;%光伏的使用年限 R_pv=0.7;%光伏单台置换费用（万元/台） %3.柴油发电机 a=0.08415;%截距系数/(L/kWh) b=0.246;%斜率/(L/kWh) P_deo=50;%额定功率 P_deo_min=10;%最小运行功率 u_de=2.09;%柴油发电机投资成本（台/万元） k_de=0.05;%柴油发电机维护费用（万元/年） Y_de=10;%柴油发电机的使用年限 R_de=1.8;%柴油发电机单台置换费用（万元/台） K_fuel=8.38;%柴油的价格8.38（元/L） V_CO2=0.210;%CO2治理成本（元/kg） V_SO2=14.842;%SO2治理成本（元/kg） V_NOx=62.964;%NOx治理成本（元/kg） kr_CO2=649;%CO2污染物排放系数（g/kw·h） kr_SO2=0.206;%SO2污染物排放系数（g/kw·h） kr_NOx=9.890;%NOx污染物排放系数（g/kw·h） %4.蓄电池 E_sb=2;%蓄电池容量 P_cha_max=0.2*num_sb*E_sb;%蓄电池最大充电功率 P_dis_max=0.2*num_sb*E_sb;%蓄电池最大放电功率 c_cha=0.85;%充电效率 c_dis=0.85;%放电效率 Soc_max=0.9;%荷电状态上限 Soc_min=0.2;%荷电状态下限 Soc_0=0.5;%荷电状态初始值 u_sb=0.16;%蓄电池投资成本（台/万元） k_sb=0.002;%蓄电池维护费用（万元/年） Y_sb=1.36;%蓄电池的使用年限 R_sb=0.09;%蓄电池单台置换费用（万元/台） %5.其他参数 r=0.0475;%成本回收系数 Y=20;%运行年限 %% 变量区 P_deo_rel=zeros(1,8760);%柴油发电机实际功率 P_cha=zeros(1,8760);%蓄电池充电功率 P_dis=zeros(1,8760);%蓄电池放电功率 P_bre=zeros(1,8760);%负荷中断量 P_was=zeros(1,8760);%弃风弃光量 SOC=zeros(1,8760);%蓄电池荷电状态 %% 计算区 %根据风速求单台风机的年功率曲线 for i=1:8760 if Speed_WT(i)>=0 && Speed_WT(i)<=v_ci P_wt(i)=0; elseif Speed_WT(i)>=v_ci && Speed_WT(i)<=v_n P_wt(i)=single_WT*(Speed_WT(i)-v_ci)/(v_n-v_ci); elseif Speed_WT(i)>=v_n && Speed_WT(i)<=v_co P_wt(i)=single_WT; else P_wt(i)=0; end end %根据光照强度和温度求单块光伏板的年功率曲线 for i=1:8760 P_pv(i)=(single_PV*Solar_PV(i)/G_n)*(1+namd*(T_c(i)-T_r)); end %% 功率分配 P_exc=Load-(num_wt*P_wt+num_pv*P_pv); for i=1:8760 if i==1 if P_exc(i)>=0 %负荷功率缺额 if P_exc(i)<=P_dis_max %case1:净负荷小于蓄电池最大功率且soc允许，由蓄电池补充缺额 P_dis(i)=P_exc(i); SOC(i)=Soc_0+((P_cha(i)*c_cha)-(P_dis(i)/c_dis))/(num_sb*E_sb);%SOC if SOC(i)<Soc_min SOC(i)=Soc_0; P_dis(i)=0; if P_exc(i)<num_g*P_deo_min P_deo_rel(i)=num_g*P_deo_min; P_cha(i)=num_g*P_deo_min-P_exc(i); else P_deo_rel(i)=P_exc(i); end end elseif P_exc(i)>