%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('---------
【配电】基于Matlab实现多个配电发电机 (DG) 的径向配电系统 上传.zip （9个子文件） 
ybuspg_33.m 1KB ybuspg_ds33.m 1KB
1.png 93KB
33 budi_test_case.m 15KB Ldat33.m 972B Bdat33.m 812B budi_test_case_33_withdg.fig 48KB budi_disrect_test_case_for_voltage_33.m 8KB Itrdat33.m 30B资源评论
天天Matlab科研工作室
- 粉丝: 4w+
- 资源: 1万+
下载权益
C知道特权
VIP文章
课程特权
VIP享7折，此内容立减5.97元 开通VIP
最新资源
资源上传下载、课程学习等过程中有任何疑问或建议，欢迎提出宝贵意见哦~我们会及时处理！
点击此处反馈
