%%改进粒子群算法
clear all
clc
tic
%%参数设置
maxgen=600; %迭代次数
sizepop=2000; %种群大小
c=1.49445; %速度更新参数，c1=c2
wmax=0.9; %惯性权重
wmin=0.4;
popmax=[2.1.*ones(1,24) 8]; %
popmin=[zeros(1,24) 1]; %
vmax=[0.1.*ones(1,24) 1]; %速度上限
vmin=[-0.1.*ones(1,24) -1]; %速度下限
dim=25; %变量维数
P_buy=[0.300000000000000,0.300000000000000,0.300000000000000,0.300000000000000,0.300000000000000,0.300000000000000,0.600000000000000,0.600000000000000,0.600000000000000,0.600000000000000,0.600000000000000,0.600000000000000,0.600000000000000,0.450000000000000,0.450000000000000,0.450000000000000,0.450000000000000,0.600000000000000,0.600000000000000,0.600000000000000,0.600000000000000,0.600000000000000,0.300000000000000,0.300000000000000];
Pload=50.*[0.170000000000000,0.170000000000000,0.170000000000000,0.170000000000000,0.270000000000000,4.25000000000000,0.170000000000000,0.170000000000000,0.170000000000000,0.170000000000000,2.35000000000000,0.370000000000000,0.170000000000000,0.170000000000000,0.170000000000000,0.170000000000000,0.170000000000000,5.56000000000000,3.37000000000000,0.370000000000000,0.370000000000000,0.370000000000000,0.170000000000000,0.170000000000000];
Pv=[0,0,0,0,0,0.300000000000000,0.800000000000000,1.50000000000000,2,2.30000000000000,2.60000000000000,2.30000000000000,2.20000000000000,2,1.50000000000000,1,0.500000000000000,0,0,0,0,0,0,0];
Tout=[21,20,19,18,18,19,20,22,24,25,27,28,30,31,31,29,28,27,26,25,24,23,22,22];%%%%%%%
wh=(wmax+wmin)/2;
%%%%%%%%%
%%初始化种群
for i=1:sizepop
pop(i,1:dim)=rand(1,dim).*(popmax-popmin)+popmin;
pop(i,25)=round(pop(i,25));
for vi=1:25
if pop(i,vi)>popmax(vi)
pop(i,vi)=popmax(vi);
elseif pop(i,vi)<popmin(vi)
pop(i,vi)=popmin(vi);
end
end
v(i,1:dim)=0.2*(2.*rand(1,dim)-1);
%计算适应度值
fitness(i) =mg_fit1(pop(i,:));
end
%%%%
detaf=0;
f=zeros(sizepop,1);
%zuixiaosherngjing deta
shengjing=sum((popmax-popmin).^2);
d=zeros(sizepop,sizepop);
sd=zeros(sizepop,sizepop);
sdd=zeros(sizepop,1);
fshengjing=zeros(sizepop,1);
for i=1:sizepop
fh=mean(fitness);
%%%%%%
f(i)=fitness(i)-fh;
if f(i)<1
f(i)=1;
end
detaf=detaf+abs(fitness(i)-fh)/f(i)/sizepop;
%%%%%%%%%%%zuixiaoshengjing
for j=1:sizepop
for k=1:12
d(i,j)=d(i,j)+(pop(i,k)-pop(j,k))^2;
end
sd(i,j)=1-(d(i,j)/shengjing)^2;
sdd(i)=sdd(i)+sd(i,j);
end
fshengjing(i)=1/sdd(i);
end
pop=[pop fshengjing];
bpop=sortrows(pop,26);
eb=1;
[best index]=min(fitness);
gtsite=pop(:,1:25); %个体最佳位置
gtfit=fitness; %个体最佳适应度
popsite=pop(index,1:25); %全局最佳位置
popfit=fitness(index); %全局最佳适应度
trace(1)=popfit;
pop=bpop(sizepop/2+1:sizepop,1:25);
% figure(1)
% scatter3(pop(:,1),pop(:,dim(1)+1),pop(:,dim(1)+dim(2)+1),'r');
% title('粒子分布图');
lame=zeros(maxgen+1,1);
lame(1)=rand;
for j=1:maxgen
if detaf<=eb
w=wmax-(wmax-wmin)*(1/(1+exp(-j/maxgen)))+wh*(4*lame(j)*(1-lame(j)));
else
w=wmax-(wmax-wmin)*(1/(1+exp(-j/maxgen)));
end
%w=wmax-(wmax-wmin)/maxgen*j;
for i=1:sizepop/2
%速度更新
v(i,:)=w*v(i,:)+c*rand*(gtsite(i,:)-pop(i,:))+c*rand*(popsite-pop(i,:));
%限制微粒速度
for vi=1:25
if v(i,vi)>vmax(vi)
v(i,vi)=vmax(vi);
elseif v(i,vi)<vmin(vi)
v(i,vi)=vmin(vi);
end
end
%位置更新
pop(i,:)=pop(i,:)+v(i,:);
pop(i,25)=round(pop(i,25));
%限制微粒位置
for vi=1:25
if v(i,vi)>vmax(vi)
v(i,vi)=vmax(vi);
elseif v(i,vi)<vmin(vi)
v(i,vi)=vmin(vi);
end
end
for vi=1:25
if pop(i,vi)>popmax(vi)
pop(i,vi)=popmax(vi);
elseif pop(i,vi)<popmin(vi)
pop(i,vi)=popmin(vi);
end
end
end
for i=1:sizepop/2
%计算适应度
fitness=mg_fit1(pop(i,:));
%个体适应度及位置更新
if fitness<gtfit(i)
gtfit(i)=fitness;
gtsite(i,:)=pop(i,:);
end
%全局适应度及位置更新
if fitness<popfit
popfit=fitness;
popsite=pop(i,:);
end
end
trace(j+1)=popfit; %记录每次迭代的最佳值
end
figure;
plot(trace,'b-*')
xlabel('迭代次数');
ylabel('适应度');
figure;
plot(Pv,'r-o')
xlabel('时间');
ylabel('典型日单位光伏发电量');
%%%%%%%%%%%%%%%%结果
kongtiao=popsite(1:24);
sx=1:24;
figure;
plot(sx,kongtiao,'r')
xlabel('时间','Fontname','宋体','Fontsize',7.5)
ylabel('功率','Fontname','宋体','Fontsize',7.5)
yyaxis right
stairs(sx,P_buy,'k--');
ylim(0:1.2);
xlabel('时间','Fontname','宋体','Fontsize',7.5)
ylabel('电价','Fontname','宋体','Fontsize',7.5)
legend('空调','电价');
Tin0=24;
tl=23;th=26;
%温度
alfa=0.2;beta=2;
for i=1:24
if i==1
Tin(i)=Tin0+alfa*(Tout(i)-Tin0)-sign(Tout(i)-Tin0)*beta*kongtiao(i);
else
Tin(i)=Tin(i-1)+alfa*(Tout(i)-Tin(i-1))-sign(Tout(i)-Tin0)*beta*kongtiao(i);
end
end
figure;
plot(Tin,'r-')
hold on
plot(Tout,'b-')
plot(tl.*ones(1,24),'--')
plot(th.*ones(1,24),'--')
xlabel('时间','Fontname','宋体','Fontsize',7.5)
ylabel('温度','Fontname','宋体','Fontsize',7.5)
yyaxis right
plot(kongtiao);
legend('室内温度','室外温度','下限温度','上限温度','空调功率');
xlabel('时间','Fontname','宋体','Fontsize',7.5)
ylabel('功率','Fontname','宋体','Fontsize',7.5)
- 1
- 2
前往页