具有多个配电发电机 (DG) 的径向配电系统附matlab代码.zip

%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); f=[1:NB]'; for ii=1:NB g=find(F(:,:)==ii); h(ii)=length(g); end k(:,1)=f; k(:,2)=h'; 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)); for ii=1:N a=fb(ii,1); b=tb(ii,1); A(ii,a)=-1; A(ii,b)=1; end A(:,1)=[]; bibc=zeros(size(LD,1)); bibc=(inv(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); Z=complex(rj,xj); 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) 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)); Ps=NBD(:,5)-NBD(:,7); PS=((Ps)/bdh); % schedule powers and Qs=NBD(:,6)-NBD(:,8); QS=((Qs)/bdh); 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, power factor is at unity, i.e., PFDG = 1, a = 0. From below eq, the optimal size of DG at % % each bus i for minimizing losses can be given by reduced equation %PDGi=PDi-1/aii[(bii*QDi)+sum(j=iandj~=i)(aij*Pj-bij*Qj)]% for Xi and Yi inthe esystem % PDGi=PDi-1/aii[(bii*QDi)+sum(j=iandj~=i)(aij*Pj-bij*Qj)]% for Xi and Yi inthe esystem pdg=BD(:,5); qdg=BD(:,6); apf=(tan(acos(0.8225))); 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 nad xi POD=BD(:,7)/bdh; QOD=BD(:,8)/bdh; Pi=((pdg)-(POD)); Qi=((((pdg))*apf)-(QOD)); 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 int he 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); 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(index56); Ploss=(Plos_type3)*bdh; plot(STABILITYINDEX,'DisplayName','STABILITYINDEX');hold all; plot(VBN,'DisplayName','VBN');hold all; plot(STABILITYINDEX1,'DisplayName','STABILITYINDEX1');hold all; plot(VBN1,'DisplayName','VBN1');hold off; Pin=(Pi)*bdh; Qin=Qi*bdh; fpdgi=val56*bdh; for ii=1:NB PDfg(ii)=abs(pdg_type3(ii))/sqrt(((abs(pdg_type3(ii)))^2)+((abs(qdg_type3(ii)))^2)); end Pdgi=(abs(pdg_type3)*bdh); Qdgi=abs(qdg_type3)*bdh; %% cost of energy loss with out DG Eloss_1dg=sum(val56)*bdh; EC=4.63 ; % http://www.aperc.gov.in/TariffOrders/Orders/2013/RST&termsandconditions13-14.pdf % Closs in rupeessthere is continurs load int he system as the load curve is constant with year 365 days 24/7 Closs_1dg=Eloss_1dg*(EC*8760) %% DG cost F_Pdgi=Pdgi(index56); ai=0;bi=20;ci=0.25; cost=ai*(a*(F_Pdgi)^2)+bi*(F_Pdgi)+ci %% output file disp('#########################################################################################'); disp('---------