基于Matlab模拟分布式电源的33节点配电网位置容量优化

%FOR A 33 bus test ssystem 12.62KV AND 420 KVA SYSTEM WITH PER UNIT WIHT A %PF OF 0.85 tic clc clear all LD=Ldat33; BD=Bdat33; bdh=max(BD(:,7));%%最大有功负荷 pu=bdh/(((12.62)^2)*1000);%%标幺值 bus=BD(:,1);%节点编号 rj=LD(:,4)*pu;%电阻 xj=LD(:,5)*pu;%电抗 F=LD(:,2:3);%支路端点矩阵 M=max(LD(:,2:3));%支路终点 NB=max(M);%%NB=33 f=[1:NB]'; for ii=1:NB g=find(F(:,:)==ii);%%find找到F矩阵中与ii相等的数的序号组成的矩阵 h(ii)=length(g);%%计算g矩阵包含的数字个数并组成h矩阵 end%%找相邻节点的个数 k(:,1)=f; k(:,2)=h';%k矩阵第一列为1-33，第二列为每个节点相邻节点的个数 CG=unique(LD(:,2:3));%%各节点编号 cbus=min(CG);%%最小编号节点 %% Nanch matrix formation NBD=BD; NBD(cbus,:)=[]; REAL_POWER=BD(:,7);%有功负荷 REACTIVE_POWER=BD(:,8);%无功负荷 fb=LD(:,2);%起始节点矩阵 tb=LD(:,3);%终止节点矩阵 N=length(LD(:,1));%两节点之间线路条数N=32 for ii=1:N a=fb(ii,1);%起始节点 b=tb(ii,1);%终止节点 A(ii,a)=-1; A(ii,b)=1;%%A矩阵中，在电流流过的方向上，相邻两条线路在后面的为-1 end A(:,1)=[]; bibc=zeros(size(LD,1));%定义bibc为32*32维矩阵 bibc=(inv(A)');%求A'的逆，回代矩阵 %% fbs load flow Ps=NBD(:,5)-NBD(:,7); %节点需要满足的有功负荷 PS=(abs(Ps)/bdh); % schedule powers and /节点需要的功率 Qs=NBD(:,6)-NBD(:,8); %节点需要满足的无功负荷 QS=(abs(Qs)/bdh); % initial values /占节点总负荷的百分比 S=complex(PS,QS);%视在功率S=PS+i QS Z=complex(rj,xj);%阻抗Z=r+i x Vo=NBD(:,3);%初始电压（标幺值） VBN=1; ZD=diag(Z);%对角矩阵 Vb=complex(NBD(:,3)); Vprev=Vb; toler=1;%精准度 iter=1;%循环次数 ZB=(ZD*(bibc))'; while (toler>0.0001)%精准到0.0001 for ii=1:N IB(ii)=conj(complex(PS(ii),QS(ii))/Vo(ii)); end B=bibc*IB'; %计 LI=bibc*abs(IB)';%算 BCBV=ZB*bibc*IB';%前推矩阵 for ii=1:N Vo(ii)=Vb(ii)-BCBV(ii);%得到新的节点电压 end iter=iter+1; %循环加一 toler=max(abs(abs(Vo)-abs(Vprev))); Vprev=Vo; end del1=angle(Vo);%%电压的相角弧度值 Vbus=abs(Vo); VBN=[VBN;Vbus]; del2=0; del2=[del2;del1]; iter [val1, index1]=min(VBN);%%最小电压及节点 del=(180/pi)*del2;%角度值 %% loss and stability index br=max(LD(:,1));%%编号=32 Ps=NBD(:,5)-NBD(:,7); %需要补偿的有功 PS=((Ps)/bdh); % schedule powers and Qs=NBD(:,6)-NBD(:,8); %需要补偿的无功 QS=((Qs)/bdh); % initial values pm2=bibc*(PS);%% qm2=bibc*(QS);%%支路功率？ LP=[]; QP=[]; for ii=1:N lpp=rj(ii)*((((pm2(ii))^2)+(qm2(ii))^2))/Vbus(ii)^2; LP=[LP;lpp]; qpp=xj(ii)*((((pm2(ii))^2)+(qm2(ii))^2))/Vbus(ii)^2; QP=[QP;qpp]; %#ok<*AGROW> end for ii=1:N Pfg(ii)=abs(PS(ii))/realsqrt((PS(ii)^2)+((QS(ii)^2))); %#ok<SAGROW>%%功率因数 end PM2=max(abs(pm2)); PM2=[PM2;pm2]; QM2=max(qm2); QM2=[QM2;qm2]; for ii=1:N p=fb(ii);%起点 q=tb(ii);%终点 R(q)=(Vbus(p)^4)-(4.0*(PM2(q)*xj(p)-QM2(q)*rj(p))^2)-(4.0*(PM2(q)*rj(p)+QM2(q)*xj(p))*(Vbus(p)^2)); %电压值 end R(:,1)=1; STABILITYINDEX=(R)'; [val2, index2]=min(STABILITYINDEX);%%电压偏差最小值及节点 Pd=PS*bdh; Qd=QS*bdh;%节点需要的功率 %% optimum allocation of distribution genarators分布式电源位置优化 % % Type DG: For Type 1 DG, 功率因数统一, i.e., PFDG = 1, a = 0. From below eq, the 最佳容量 of DG at % % each bus i for 减少损耗 can be given by 简化公式 %PDGi=PDi-1/aii[(bii*QDi)+sum(j=iandj~=i)(aij*Pj-bij*Qj)]% for Xi and Yi in the esystem % PDGi=PDi-1/aii[(bii*QDi)+sum(j=iandj~=i)(aij*Pj-bij*Qj)]% for Xi and Yi in the esystem pdg=BD(:,5); qdg=BD(:,6);%电源功率 apf=(tan(acos(0.8225)));%%Q/P ybus=ybuspg_33;%%导纳阵 fff=sparse(ybus);%%生成稀疏矩阵 Zbus=inv(ybus);%%逆 r=(real(Zbus)); %实部 x=(imag(Zbus)); %虚部 F=LD(:,2:3); %PS and QS are specified powers %% for the calculation of aij aij=zeros(NB,NB); for ii=1:NB for jj=1:NB aij(ii,jj)=(r(ii,jj)*(cos(del2(ii)-del2(jj))/(VBN(ii)*VBN(jj)))); cij(ii,jj)=(x(ii,jj)*(cos(del2(ii)-del2(jj))/(VBN(ii)*VBN(jj)))); end end %% for the calculation of bij bij=zeros(NB,NB); for jj=1:NB for ii=1:NB bij(ii,jj)=(r(ii,jj)*(sin(del2(ii)-del2(jj))/(VBN(ii)*VBN(jj)))); dij(ii,jj)=(x(ii,jj)*(sin(del2(ii)-del2(jj))/(VBN(ii)*VBN(jj)))); end end %% for the calculation of xi and yi POD=BD(:,7)/bdh;%有功负荷/最大有功负荷 QOD=BD(:,8)/bdh; %无功负荷/最大无功负荷 Pi=((pdg)-(POD)); Qi=((((pdg))*apf)-(QOD)); %也可以写Qi=（（qdg）-（QOD）），前面定义过qdg XI=[]; for ii=1:NB xi=0; for jj=1:NB if jj~=ii xi=xi+(aij(ii,jj)*Pi(jj))-(bij(ii,jj)*Qi(jj)); end end XI=[XI;xi]; end YI=[]; for ii=1:NB yi=0; for jj=1:NB if jj~=ii yi=yi+(aij(ii,jj)*Qi(jj))+(bij(ii,jj)*Pi(jj)); end end YI=[YI;yi]; end %% without dg losses plosWG=[]; pdg=zeros(NB,1); qdg=zeros(NB,1); Pi=((pdg)-(POD)); Qi=((((pdg))*apf)-(QOD)); Pj=((pdg)-(POD)); Qj=(((pdg)*apf)-(QOD)); for ii=1:NB ploss=0; for jj=1:NB ploss=ploss+(aij(ii,jj)*((Pi(ii)*Pj(jj))+(Qi(ii)*Qj(jj)))+(bij(ii,jj)*((Qi(ii)*Pj(jj))-(Pi(ii)*Qj(jj))))); end plosWG=[plosWG;ploss]; end [val3, index3]=min(plosWG);%%网损最小值和节点 %% cost of energy loss with out DG Eloss=sum(plosWG)*bdh; EC=4.63 ; % http://www.aperc.gov.in/TariffOrders/Orders/2013/RST&termsandconditions13-14.pdf % Closs in rupeessthere is continurs load in the system as the load curve is constant with year 365 days 24/7 Closs=Eloss*(EC*8760); %% for the calculation of pdgi and plos,qlos for type_3 pdg_type3=BD(:,4);%%当作pq节点类型 for dg=1:NB pdg_type3(dg)=((aij(dg,dg)*(POD(dg)+(QOD(dg)*apf)))-(YI(dg)*apf)-XI(dg))/(((aij(dg,dg)*apf^2)+aij(dg,dg))); end for dg=1:NB qdg_type3(dg)=apf*pdg_type3(dg); end data2=zeros(NB,1); Sdgi1=zeros(N,1); %% after placing dg for kk=1:NB clear data1 data1(kk)=pdg_type3(kk); clear data2 data2=zeros(NB,1); data2(kk,1)=data1(kk); Plos_type3=[]; Pi=((data2)-(POD)); Qi=((((data2))*apf)-(QOD)); Pj=((data2)-(POD)); Qj=(((data2)*apf)-(QOD)); for m=1:NB ploss_type3=0; for n=1:NB ploss_type3=ploss_type3+((aij(m,n)*((Pi(m)*Pj(n))+(Qi(m)*Qj(n))))+(bij(m,n)*((Qi(m)*Pj(n))-(Pi(m)*Qj(n))))); end Plos_type3=[Plos_type3;ploss_type3]; end Sum_Plos_type3(kk)=sum(Plos_type3); end [val56,index56]=min(Sum_Plos_type3); [val57,index57]=max(Sum_Plos_type3);%%有功网损最大值 MfdP2=zeros(NB,1); MfdP2(index56,1)=pdg_type3(index56); MfdQ2=zeros(NB,1); MfdQ2(index56,1)=qdg_type3(index56); MfdQ2(cbus,:)=[]; MfdP2(cbus,:)=[]; Ps=MfdP2-NBD(:,7); PS1=(abs(Ps)/bdh); % schedule powers and Qs1=MfdQ2-NBD(:,8); QS1=(abs(Qs1)/bdh);% initial values %% dg placement voltage Vbus1=Vbus; S1=complex(PS1,QS1); Z1=complex(rj,xj); Vo2=Vbus1; VBN1=1; ZD1=diag(Z1); SDG=complex(MfdP2,MfdQ2); Vb1=ones(N,1); Vprev1=Vb1; toler=1;iter=1; ZB1=(ZD1*(bibc))'; while max(abs(Vo2))<=1.000&&max(abs(Vo2>=0.9))&&(toler>0.001) for jj=1:N IB(jj)=conj((S1(jj))/Vo2(jj)); Idg(jj)=conj((SDG(jj))/Vo2(jj)); Ii1(jj)=-Idg(jj)+IB(jj); end B=bibc*Ii1'; LI1=bibc*(Ii1)'; BCBV1=ZB1*bibc*Ii1'; for ff=1:N Vo2(ff)=Vb1(ff)-BCBV1(ff); end iter=iter+1; toler=max(abs(abs(Vo2)-abs(Vprev1))); Vprev1=Vo2; end del11=angle(Vo2); Vbus1=abs(Vo2); VBN1=[VBN1;Vbus1]; del21=0; del21=[del21;del11]; [val17, index17]=min(VBN1); del2=(180/pi)*del21; Ploss_type3=[]; %% PM2=max(pm2); PM2=[PM2;pm2]; QM2=max(qm2); QM2=[QM2;qm2]; for ii=1:N p=fb(ii);q=tb(ii); R1(q)=(Vbus1(p)^4)-(4.0*(PM2(q)*xj(p)-QM2(q)*rj(p))^2)-(4.0*(PM2(q)*rj(p)+QM2(q)*xj(p))*(Vbus1(p)^2)); end R1(:,1)=1; STABILITYINDEX1=(R1)'; [val21, index21]=min(STABILITYINDEX1); Pdgi=(abs(pdg_type3)*bdh); Qdgi=abs(qdg_type3)*bdh; Sdgi=abs(complex(Pdgi,Qdgi')); Pdgi_type3=abs(pdg_type3)*bdh; Qdgi_type3=abs(qdg_type3)*bdh; Real_TYPE_3=Pdgi_type3(inde