function microgird_dispatch_2()
%第一行 DG,第二行MT,第三行FC
tic;
w_start=0.9;
w_end=0.4;
c1=2;
c2=2;
Vmax=200;
PopSize=20;
MaxIter=300;
iter=0;
T=24;
N=4; %可再生能源的种类数
% 输入原始数据,包括各时段负荷大小,发电机有功输出上下限,发电机耗量成本系数,各时段风电场预测的平均输出功率
pmax(1)=7;pmin(1)=0;
pmax(2)=4;pmin(2)=0;
pmax(3)=4;pmin(3)=0;
% a=[1000,970,700,680,450,370,480,660,665,670];
% b=[16.19,17.26,16.60,16.50,19.70,22.26,27.74,25.92,27.27,27.79];
% c=[0.00048,0.00031,0.002,0.00211,0.00398,0.00712,0.00079,0.00413,0.00222,0.00173];
% 柴油发电机的耗量参数
a=0.0071;b=0.2333;c=0.4333;
% 微型燃气轮机,P_mt表示燃气轮机发出的功率,Xl_mt表示燃气轮机的效率
% Xl_mt=0.0753*(P_mt/65)^3-0.3095*(P_mt/65)^2+0.1068;
% 燃料电池
Xl_fc=0.4;
% Cost_fc=Price_fc*P_fc/Xl_fc;
xx=1:T;
%负荷数据
% PL=[11,4,4,4,5,6,7.5,8,7.5,6,5,6,6.5,7,8,5,8.5,11,13,13,10,10,6.5,5];
PL=[3,4,4,4,5,6,6.5,7,7.5,8.5,9,10,10.5,10,9,8.5,9,10,11,11.5,10,9,5.5,5];
%PL=PL*10;
plot(xx,PL,'-*r');
xlabel('t/h');ylabel('P/kW');
axis([1,24,0,15]);
% 光伏和风力发电24小时的预测数据
PV=[0,0,0,0,0.0067,0.0133,0.01661,0.0298,0.043,0.0446,0.0529,0.05605,0.0543,0.0558,0.046,0.0394,0.0164,0.0131,0.0099,0.0049,0.0033,0,0,0];
WT=[0.2762,0.2762,0.3077,0.2824,0.3209,0.3275,0.2701,0.2824,0.2348,0.2521,0.2701,0.3275,0.3275,0.3479,0.3209,0.3077,0.3209,0.3342,0.3618,0.334,0.321,0.307,0.3,0.308];
PV_WT=PV(1:T)+WT(1:T);
% figure(1);
% plot(xx,PV,'-*r',xx,WT,'-b^');
% xlabel('t/h');ylabel('P/kW');
PGTTT=zeros(N,T);
PGTTT(N,:)=PV_WT;
TTT=T;
%_________________________________________________________________________
%_________________________________________________________________________
for iii=1:TTT
T=1;
B=ones(N,T,PopSize);
% B=zeros(N,T,PopSize);
% for i=1:PopSize
% B(:,:,i)=Y; %生成一个(N,T,PopSize)维的矩阵,生成了PopSize个11*24的矩阵,构成了A
% end
%随机初始化
Z=zeros(N,T,PopSize);
Z(N,:,:)=PV_WT(iii);
for jj=1:PopSize;
X=zeros(N,T);
for i=2:3
for j=1:T
X(i,j)=rand(1)*(pmax(i)-pmin(i))+pmin(i);
% X(i,j)=Y(i,j)*(rand(1)*(pmax(i)-pmin(i))+pmin(i));
end
end
X(N,:)=PV_WT(iii);
% X(1,:)=PL(1:T)-sum(X(2:4,:));
X(1,:)=PL(iii)-sum(X(2:4,:));
k= X(1,:); % 对第一行平衡机组越限的调整,取其限值
change=k>pmax(1);
k(find(change))=pmax(1);
change=k<pmin(1);
k(find(change))=pmin(1);
X(1,:)=k;
Z(:,:,jj)=X;
end
Z;
M=Z; % 用在之后第一次迭代更新的时候,如果某个粒子不满足要求,则用初始化时的粒子替代
%适应值函数
Cost=zeros(PopSize,1);
for jj=1:PopSize
Cost(jj)=shiyingzhi_fadianchengben_2(Z(:,:,jj),iii);
% Cost(jj)=shiyingzhi_pollution_2(Z(:,:,jj),iii);
% Cost(jj)=shiyingzhi_zonghexiaoyi_2(Z(:,:,jj),iii);
end
V=rand(4,T,PopSize);%微粒速度随机初始化
%设定当前位置为粒子的最好位置,并记录其最好值。
PBest=Z;
FPBest=Cost;
%找到初始微粒群体的最好微粒
[Fgbest,r]=min(Cost);
CF=Fgbest;%记录当前全局最优值
Best=Z(:,:,r);%用于保存最优粒子的位置
FBest=Fgbest;
HengZB=1:MaxIter+1;
ZongZB=zeros(1,MaxIter);
%-------------------------------------------------------------------------------------------------------
iter=0;
%循环
while(iter<=MaxIter)
iter=iter+1;
% 更新惯性权重的值
w_now=((w_start-w_end)*(MaxIter-iter)/MaxIter)+w_end;
A=Z(:,:,r);
for i=1:PopSize
A(:,:,i)=Z(:,:,r);%生成一个(11,24,PopSize)维的矩阵,生成了PopSize个11*24的矩阵,构成了A
end
%生成随机数
R1=rand(N,T,PopSize);
R2=rand(N,T,PopSize);
%速度更新
V=w_now*V+c1*R1.*(PBest-Z)+c2*R2.*(A-Z);
%对进化后速度小于最小速度的微粒进行处理
changeRows=V>Vmax;
V(find(changeRows))=Vmax;
changeRows=V<-Vmax;
V(find(changeRows))=-Vmax;
%微粒位置进行更新
Z=Z+1.0*V.*B;
%判断更新后的粒子是否满足约束条件,对不满足约束的粒子进行处理,使其满足约束,如果不满足则重新生成粒子的速度继而更新其位置,直到满足约束条件,如果重复
%次数超过规定的次数,则用原来可行的该粒子代替。
for i=1:PopSize
for j=2:3
k=Z(j,:,i);
change=k>pmax(j);
k(find(change))=pmax(j);
change=k<pmin(j);
k(find(change))=pmin(j);
Z(j,:,i)=k;
end
Z(N,:,i)=PV_WT(iii);
% Z(1,:,i)=PL(1:T)-sum(Z(2:4,:,i));
Z(1,:,i)=PL(iii)-sum(Z(2:4,:,i));
k= Z(1,:,i); % 对第一行平衡机组越限的调整,取其限值
change=k>pmax(1);
k(find(change))=pmax(1);
change=k<pmin(1);
k(find(change))=pmin(1);
Z(1,:,i)=k;
end
%重新计算新位置的适应度值
Cost=zeros(PopSize,1);
for jj=1:PopSize
Cost(jj)=shiyingzhi_fadianchengben_2(Z(:,:,jj),iii);
% Cost(jj)=shiyingzhi_pollution_2(Z(:,:,jj),iii);
% Cost(jj)=shiyingzhi_zonghexiaoyi_2(Z(:,:,jj),iii);
end
%更新每个微粒的最好位置
d=Cost;
P=d<FPBest;
FPBest(find(P))=d(find(P));%适应值更换
PBest(:,:,find(P))=Z(:,:,find(P));%粒子位置更换
[Fgbest,g]=min(FPBest);%保存最好适应值
Best=PBest(:,:,g);%自己加的一句,因为如果下边的if语句如果不成立的话,Best就不更新了,不合逻辑
if Fgbest<CF%如果本次适应值好于上次则保存
[FBbest,h]=min(FPBest);%最好适应值为FBest
Best=PBest(:,:,h);
end
CF=Fgbest;
ZongZB(iter)=CF;
end
CF;
Fgbest;
Best;
PGTTT(:,iii)=Best;
cccc=sum(Best);
figure(5);
plot(HengZB,ZongZB);
xlabel('迭代次数');ylabel('全局最优值/$');
end
aaaaaa=sum(PGTTT);
Cost=shiyingzhi_fadianchengben_all_2(PGTTT)
% Cost=shiyingzhi_pollution_all_2(PGTTT)
% Cost=shiyingzhi_zonghexiaoyi_all_2(PGTTT)
%__________________________________________________________________________
%__________________________________________________________________________
% figure(2);
% plot(HengZB,ZongZB);
PGTTT % 输出最优的各类型微电源结果
figure(3);
plot(xx,PGTTT(1,:),'-ksquare',xx,PGTTT(2,:),'-r*',xx,PGTTT(3,:),'-b^'); % 绘制24小时的各微电源的功率曲线图
xlabel('t/h');ylabel('P/kW');
axis([1,24,-1,8]);
figure(4);
plot(xx,sum(Best),':',xx,PL(1:T),'*',xx,sum(Best)-PL(1:T),'x');
toc;