clear
clc
format long;
%初始化条件****************************************
%微型燃气轮机最大功率
MTMaxPower=5;
%微型燃气轮机最小功率
MTMinPower=1;
%电网输入微网最大功率
GridMaxImportPower=250;
%电网输入微网最小功率
GridMinImportPower=10;
%储能最大放电功率
StorageMaxDischargingPower=150;
%储能最大充电功率
StorageMaxChargingPower=-100;
Max_Dt=300;%最大迭代次数300
D=72;%搜索空间维数（未知数个数）
N=600;%粒子个数600
w_max=0.9;
w_min=0.4;
v_max=2;
s=1;
%初始化种群个体（位置和速度）***********************
for i=1:N
for j=1:72
% v(i,j)=randn;
v(i,j)=0.0;
if j<25
x(i,j)=MTMinPower+rand()*(MTMaxPower-MTMinPower);
elseif j>24&&j<49
x(i,j)=GridMinImportPower+rand()*(GridMaxImportPower-GridMinImportPower);
elseif j>48&&j<73
x(i,j)=StorageMaxChargingPower+rand()*(StorageMaxDischargingPower-StorageMaxChargingPower);
end
end
end
%计算各个粒子的适应度，并初始化Pi和Pg****************
for i=1:N
p(i)=fitness(x(i,:),s);
y(i,:)=x(i,:);%每个粒子的个体寻优值
end
Pbest=fitness(x(1,:),s);
pg=x(1,:);%Pg为全局最优
for i=2:N
if fitness(x(i,:),s)<fitness(pg,s)
Pbest=fitness(x(i,:),s);
pg=x(i,:);%全局最优更新
end
end
%进入主循环*****************************************
for t=1:Max_Dt
for i=1:N
w=w_max-(w_max-w_min)*t/Max_Dt;%惯性权重更新
c1=(0.5-2.5)*t/Max_Dt+2.5; %认知
c2=(2.5-0.5)*t/Max_Dt+0.5; %社会认识
% w=0.7;
% c1=2.05; %认知
% c2=2.05; %社会认识
v(i,:)=w*v(i,:)+c1*rand()*(y(i,:)-x(i,:))+c2*rand()*(pg-x(i,:));
for m=1:72
if(v(i,m)>v_max)
v(i,m)=v_max;
elseif(v(i,m)<-v_max)
v(i,m)=-v_max;
end
end
x(i,:)=x(i,:)+v(i,:);
%对粒子边界处理*****************************
for n=1:72
if n<25
if x(i,n)<MTMinPower
x(i,n)=MTMinPower;
v(i,n)=-v(i,n);
elseif x(i,n)>MTMaxPower
x(i,n)=MTMaxPower;
v(i,n)=-v(i,n);
else
delt(i,n)=0;
end
elseif n>24&&n<49
if x(i,n)<GridMinImportPower
x(i,n)=GridMinImportPower;
v(i,n)=-v(i,n);
elseif x(i,n)>GridMaxImportPower
x(i,n)=GridMaxImportPower;
v(i,n)=-v(i,n);
else
delt(i,n)=0;
end
else
if x(i,n)<StorageMaxChargingPower
x(i,n)=StorageMaxChargingPower;
v(i,n)=-v(i,n);
elseif x(i,n)>StorageMaxDischargingPower
x(i,n)=StorageMaxDischargingPower;
v(i,n)=-v(i,n);
else
delt(i,n)=0;
end
end
end
%对粒子进行评价，寻找最优值******************
if fitness(x(i,:),t)<p(i)
p(i)=fitness(x(i,:),t);
y(i,:)=x(i,:);
end
if p(i)<Pbest
Pbest=p(i);
pg=y(i,:);
s=t;
end
end
%Pbest(t)=fitness(pg,s);
end
disp('*************************************************************')
disp('函数的全局最优位置为：')
Solution=pg'
for m=1:24
pg1(m)=pg(m);
end
for m=25:48
pg2(m-24)=pg(m);
end
for m=49:72
pg3(m-48)=pg(m);
end
figure
subplot(311)
plot( pg1,'-*')
xlim([1 24])
grid
title('MT运行计划')
subplot(312)
plot( pg2,'-*')
xlim([1 24])
grid
title('GRID运行计划')
subplot(313)
plot( pg3,'-*')
xlim([1 24])
grid
title('BA运行计划')
disp('最后得到的优化极值为：')
Result=economic(pg)
Resulta=fitness(pg,s)
s
disp('*************************************************************')
