syms x;
fd=(1/(x*0.88*sqrt(2*3.14159)))*exp(-(log(x)-3.2)^2/(2*0.88*0.88));%生成充电时长的概率密度函数
fd=inline('(1/(x*0.88*sqrt(2*3.14159)))*exp(-(log(x)-3.2)^2/(2*0.88*0.88))','x');
T=zeros(1,24);
k=fd(11)/normpdf(11,20,15); %计算辅助正态分布的k
N=1000; %汽车数量
td=[];
for i=1:N
accept=0;
while accept==0
u=rand();%生成判断是否拒绝的分布u
z=normrnd(15,10);
if z<=200 && z>=0 && u<=fd(z)/(k*normpdf(z,20,15))
accept=1;
td=[td;z];
end
end
end
td=30-(td*14.2/100);
%转换为剩余的电池容量 每100公里耗电量14.2Kw，满容量30KW.h
num=50;%种群数设置为50
T(11:13)=-220;
T(19)=-100;
Tmax=200;%迭代次数设置为200
X=(rand(N,24,num)*10-5);%充放电功率约束在-5到5之间
Wmax=0.9;%最大惯性权重
Wmin=0.3;%最小惯性权重
c10=0.3;%自我学习因子初值
c11=0.9;%自我学习因子终值
c20=0.9;%群体学习因子初值
c21=0.3;%群体学习因子终值
v=rand(N,24,num);%初始化种群速度
Xbest=X;%初始化每个个体的历史最佳位置
Ybest=zeros(N,24);%初始化种群的历史最佳位置
fxbest=zeros(num,1);%每个个体的历史最佳适应度
fybest=inf;%种群最佳适应度
iter=1;
Price=[0.365 0.365 0.365 0.365 0.365 0.365 0.365 0.687 0.687 0.687 1.07 1.07 1.07 1.07 0.687 0.687 0.687 0.687 1.07 1.07 1.07 1.07 0.687 0.687];%分时电价
Y1max=(3716-2345)/3716;%原始负荷谷峰差率
Y2max=sum((30-td)*1.07);%一直按照最大功率充放电的成本
T(18:20)=-200;
xl=[2500 2452 2450 2345 2446 2570 2674 3241 3242 3403 3545 3214 3255 3311 3280 3414 3501 3716 3554 3478 3457 3234 2910 2530] ;%原充电负荷
while iter<Tmax
for i=1:num
[X(:,:,i),fx(i)]=fittness(X(:,:,i),td,xl,Price,Y1max,Y2max);
end %计算每个个体的适应度
for i=1:num
if fxbest(i)>fx(i)
fxbest(i)=fx(i);
Xbest(:,:,i)=X(:,:,i);%更新个体历史最佳位置
end
end
if fybest>min(fxbest)
[fybest,nmin]=min(fxbest);
Ybest=X(:,:,nmin);
end
for i=1:num
v(:,:,i)=v(:,:,i)*(Wmax-iter*(Wmax-Wmin)/Tmax)+(c10+(c11-c10)*iter/Tmax)*rand(N,24).*(Xbest(:,:,i)-X(:,:,i))+(c20+(c21-c20)*iter/Tmax)*rand(N,24).*(Ybest-X(:,:,i));%更新速度
X(:,:,i)=X(:,:,i)+v(:,:,i);
end
iter=iter+1;
end
for i=1:24
for j=1:N
if Ybest(j,i)>=0
T(i)=T(i)+Ybest(j,i);
end
end
end
T=T+xl;
plot(xl,'b');
hold on;
plot(T,'k--');
legend('常规负荷','叠加电动汽车充电后负荷','Location','northwest');
axis([0 24 0 5500]);
xlabel('时间/(h)');
ylabel('功率/kW');
title('总负荷曲线');
set(gcf,'color',[1 1 1]);
set(gca,'xtick',0:2:24);
