牛拉法潮流计算(2)

j=0; for i=1:n

if i~=isb&B2(i,6)==2 h=h+1; for j=1:n

OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3))); end end end for i=1:n

if i~=isb&B2(i,6)==3 h=h+1; for j=1:n

OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3))); end end end % OrgS

%创建DetaS h=0; for i=1:n

if i~=isb&B2(i,6)==2 h=h+1;

DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1); end end t=0; for i=1:n

if i~=isb&B2(i,6)==3 h=h+1; t=t+1;

% DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1);

DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2; end end

% DetaS %创建I

i=zeros(n-1,1); h=0; for i=1:n if i~=isb h=h+1;

I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3)); end end % I

%创建Jacbi

Jacbi=zeros(2*n-2); h=0; k=0; for i=1:n

if B2(i,6)==2 h=h+1; for j=1:n if j~=isb k=k+1; if i==j

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1)); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1)); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); else

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3)); Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3)); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end

if k==(n-1) k=0; end end end end end k=0; for i=1:n

if B2(i,6)==3 h=h+1; for j=1:n if j~=isb k=k+1; if i==j

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1)); Jacbi(2*h,2*k-1)=2*imag(B2(i,3)); Jacbi(2*h,2*k)=2*real(B2(i,3)); else

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3)); Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3)); Jacbi(2*h,2*k-1)=0; Jacbi(2*h,2*k)=0; end

if k==(n-1) k=0; end end end end end % Jacbi

DetaU=zeros(2*n-2,1); DetaU=inv(Jacbi)*DetaS; % DetaU

%修正节点电压 j=0; for i=1:n

if B2(i,6)==2 j=j+1;

B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end end for i=1:n

if B2(i,6)==3 j=j+1;

B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end end % B2

Times=Times+1; %迭代次数加1 end

disp('迭代次数为：'); disp(Times);

disp('收敛时电压修正量为：：'); disp(DetaU);

for k=1:n

E(k)=B2(k,3);

e(k)=real(E(k)); f(k)=imag(E(k));

V(k)=sqrt(e(k)^2+f(k)^2);

sida(k)=atan(f(k)./e(k))*180./pi; end

%=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：'); disp(E); %显示各节点的实际电压标幺值E用复数表示 disp('-----------------------------------------------------');

disp('各节点的电压大小V为(节点号从小到大排列)：'); disp(V); %显示各节点的电压大小V的模值 disp('-----------------------------------------------------');

disp('各节点的电压相角sida为(节点号从小到大排列)：'); disp(sida); %显示各节点的电压相角

for p=1:n C(p)=0; for q=1:n

C(p)=C(p)+conj(Y(p,q))*conj(E(q)); %计算各节点的注入电流的共轭值 end

S(p)=E(p)*C(p); %计算各节点的功率 S = 电压 X 注入电流的共轭值 end

disp('各节点的功率S为(节点号从小到大排列)：'); disp(S); %显示各节点的注入功率 Sline=zeros(n1,5);

disp('-----------------------------------------------------');

disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：'); for i=1:n1

p=B1(i,1); q=B1(i,2);

Sline(i,1)=B1(i,1); Sline(i,2)=B1(i,2); if B1(i,6)==0

Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5)))); Siz(i)=Si(p,q); else

Si(p,q)=E(p)*(conj(E(p))*((1-B1(i,5))/B1(i,3))+(conj(E(p))-conj(E(q)))*(B1(i,5)/B1(i,3))); Siz(i)=Si(p,q); end

SSi(p,q)=Si(p,q); Sline(i,3)=Siz(i);

ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))]; disp(ZF); end

disp('-----------------------------------------------------');

disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：'); for i=1:n1

p=B1(i,1);q=B1(i,2); if B1(i,6)==0

Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5)))); Sjy(i)=Sj(q,p); else

Sj(q,p)=E(q)*(conj(E(q))*((B1(i,5)*(B1(i,5)-1))/B1(i,3))+(conj(E(q))-conj(E(p)))*(B1(i,5)/B1(i,3))); Sjy(i)=Sj(q,p); end

SSj(q,p)=Sj(q,p); Sline(i,4)=Sjy(i);

ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))]; disp(ZF); end

disp('-----------------------------------------------------');

disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：'); for i=1:n1

p=B1(i,1); q=B1(i,2);

DS(i)=Si(p,q)+Sj(q,p); DDS(i)=DS(i); Sline(i,5)=DS(i);

ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))]; disp(ZF); end

disp('-----------------------------------------------------');

disp('各支路首端编号 末端编号 首端功率 末端功率 线路损耗'); disp(Sline);