MICROGRID-IN-GRID-MODE.zip_dispatch_储能微电网_多目标 调度_电网 储能_电网调度

function microgird_dispatch_battery() % 并网运行模式 % 应用启发式规则确定下一个时段对蓄电池充电还是放电 or 利用模糊控制 % 第一行表示 从电网吸收的功率，可以从0到任意值，2,3,4行依次为FC，MT，DG tic; w_start=0.9; w_end=0.4; c1=2; c2=2; Vmax=100; PopSize=30; MaxIter=500; iter=0; T=24; N=4; %可再生能源的种类数 % 输入原始数据,包括各时段负荷大小，发电机有功输出上下限，发电机耗量成本系数，各时段风电场预测的平均输出功率 pmax(1)=40;pmin(1)=0; % DG % pmax(1)=0;pmin(1)=0; pmax(2)=40;pmin(2)=10; % FC pmax(3)=40;pmin(3)=10; % MT % 柴油发电机的耗量参数 a=0.4333;b=0.2333;c=0.0071; % 微型燃气轮机,P_mt表示燃气轮机发出的功率，Xl_mt表示燃气轮机的效率 % Xl_mt=0.0753*(P_mt/65)^3-0.3095*(P_mt/65)^2+0.4174*PGT(3)/65+0.1068; % 燃料电池 % Xl_fc=0.4; Xl_fc=0.8; % 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=[40,35,30,35,45,60,65,70,75,85,90,100,105,100,90,85,90,100,110,110,100,90,55,50]; PL=[40,35,30,35,45,60,65,70,70,80,90,100,102,100,90,85,90,110,120,115,110,90,55,50];% 这组负荷有几个时段即使所有的微电源都满发都不能满足，选用蓄电池补足 % PL_chu=PL; figure(6); plot(xx,PL,'-*r'); axis([1,24,-20,120]); % xp=sum(PL); % xpxp=xp*0.2; % PL=[3,4,4,4,5,6]; % P_battery=[1 1 1 0.5 -1 -1 -1 -0.5 1 1 1 0.5 -1 -1 -1 -0.5 1 1 1 0.5 -1 -1 -1 -0.5]; % 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=PV*10; WT=WT*10; % PV_WT=PV(1:T)+WT(1:T); % PV=[0,0,0,0,0.7,0.1,0.1,0.3,0.4,0.4,0.5,0.6,0.5,0.5,0.4,0.4,0.2,0.1,0.1,0,0.3,0,0,0]; % WT=[2.8,2.8,3,2.8,3.2,3.2,2.7,2.8,2.3,2.5,2.7,3.3,3.3,3.5,3.2,3.1,3.2,3.3,3.6,3.3,3.2,3.1,3,3.1]; PV_WT=PV(1:T)+WT(1:T); % PLall=sum(pmax); % Pcd=PL-PLall-PV_WT; % k= Pcd; % change=k<=0; % k(find(change))=0; % Pcd=k; % PL=PL-Pcd; % figure(1); % plot(xx,PV,'-*r',xx,WT,'-b^',xx,PV_WT,'-or'); % xlabel('t/h');ylabel('P/kW'); PGTTT=zeros(N,T); PGTTT(N,:)=PV_WT; TTT=T; Pbattery_start=zeros(1,T); Pbattery_end=zeros(1,T); Pcd=zeros(1,T); % %蓄电池来回充放电 % Pcd=[0.3 0.3 0.3 0.15 -0.3 -0.3 -0.3 -0.15 0.3 0.3 0.3 0.15 -0.3 -0.3 -0.3 -0.15 0.3 0.3 0.3 0.15 -0.3 -0.3 -0.3 -0.15]; % Pcd=[7 7 7 3.5 -7 -7 -7 -3.5 7 7 7 3.5 -7 -7 -7 -3.5 7 7 7 3.5 -7 -7 -7 -3.5]; % 蓄电池容量为50，初始容量为25 % 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]; % Pcd=[0, 0, 0, -1, 0, 0, 0, 0, 0.5,-0.5, 0, 0, 0.5, 0, 0, 0, -1, 0, 0, 1.5, 0.5, 0,-1.2, 0]; % Pcd=[-10,-10,-5,0,0,0, 0,0,0,0,0,0, 0,0,0,0,5,5, 5,5,5,0,0,0]; % PL=PL-Pcd; %_________________________________________________________________________ %_________________________________________________________________________ for iii=1:TTT %——————————————————— % %增加模糊控制策略，决定蓄电池的充放电 % %蓄电池剩余容量 Pbattery_all=25; Pbattery_start(1)=Pbattery_all; %Pcd表示charge或者discharge的功率 if iii~=1 Pbattery_start(iii)=Pbattery_end(iii-1); end fismat=readfis('fuzzy_2inputs_2'); Pcd(iii)=evalfis([Pbattery_start(iii) PL(iii)],fismat); if iii==22 Pcd(22)=(Pbattery_end(21)-Pbattery_all)/3; end if iii==23 Pcd(23)=Pcd(22); end if iii==24 Pcd(24)=Pcd(22); end % % Pcd(19)=6.2830; %这两个时段除了蓄电池外，其他的微电源的总出力不能满足负荷需求，强制发出这么多功率。维持功率平衡 % % Pcd(20)=1.6110; Pbattery_end(iii)=Pbattery_start(iii)-Pcd(iii); % PL(iii)=PL(iii)-Pcd(iii); %——————————————————— 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); % 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,:)); % for j=2:4 % k=X(j,:); % change=k>pmax(j); % k(find(change))=pmax(j); % change=k<pmin(j); % % k(find(change))=pmin(j); % k(find(change))=0; % X(j,:)=k; % end X(1,:)=PL(iii)-sum(X(2:N,:)); k= X(1,:); % 对第一行平衡机组越限的调整，取其限值 change=k>pmax(1); k(find(change))=pmax(1); change=k<pmin(1); % k(find(change))=pmin(1); k(find(change))=0; X(1,:)=k; Z(:,:,jj)=X; end Z; M=Z; % 用在之后第一次迭代更新的时候，如果某个粒子不满足要求，则用初始化时的粒子替代 %适应值函数 Cost=zeros(PopSize,1); for jj=1:PopSize Cost(jj)=shiyingzhi(Z(:,:,jj),iii,PL(iii)); end V=rand(N,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); k(find(change))=0; 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:N,:,i)); k= Z(1,:,i); % 对第一行平衡机组越限的调整，取其限值 change=k>pmax(1); k(find(change))=pmax(1); change=k<pmin(1); % k(find(change))=pmin(1); k(find(change))=0; Z(1,:,i)=k; end %重新计算新位置的适应度值 Cost=zeros(PopSize,1); for jj=1:PopSize Cost(jj)=shiyingzhi(Z(:,:,jj),iii,PL(iii)); end %更新每个微粒的最好位置 d=Cost; P=d<FPBest; FPBest(find(P))=d(find(P));%适应值更换 PBest(:,:,find(P))=Z(:,:,find(P));%粒子位置更换 [Fgbest,g]=mi