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)>
