功率计算程序

clear

n=4;

n1=1;

n2=2;

isb=4;

pr=0.00001;

K=[0 0 0 0;0 0 0 0.9625;0 0 0 0;0 0 0 0];

C=[0 0.02+0.06i 0.01+0.03i inf; 0.02+0.06i 0 0.03+0.07i 0.0+0.05i;0.01+0.03i 0.03+0.07i 0 0.02+0.05i;inf 0.0+0.05i 0.02+0.05i 0];

y=[0 0.01i 0.01i 0;0.01i 0 0 0;0.01i 0 0 0;0 0 0 0];

U=[1+0i 1+0i 1.02+0i 1.05+0i];

S=[-0.4-0.3i -0.3-0.2i 0.4 0];

Z=zeros(1,n);N=zeros(n1+n2,n2);L=zeros(n2,n2);QT1=zeros(1,n1+n2); for m=1:n

for R=1:n

C(m,m)=C(m,m)+y(m,R);

if K(m,R)~=0

C(m,m)=C(m,m)+1/(C(m,R) /( K(m,R) * (K(m,R)-1))) ;

C(R,R)=C(R,R)+1/(C(m,R)/(1-K(m,R)));

C(m,R)=C(m,R)/K(m,R);

C(R,m)=C(m,R);

end

end

end

for m=1:n

for R=1:n

if m~=R

Z(m)=Z(m)+1/C(m,R);

end

end

end

for m=1:n

for R=1:n

if m==R

Y(m,m)=C(m,m)+Z(m);

else

Y(m,R)=-1/C(m,R);

end

end

end

disp('结点导纳矩阵:');

disp(Y);

disp('迭代中的雅克比矩阵:');

G=real(Y);

B=imag(Y);

O=angle(U);

U1=abs(U);

k=0;

PR=1;

P=real(S);

Q=imag(S);

while PR>pr

for m=1:n2

UD(m)=U1(m);

end

for m=1:n1+n2

for R=1:n

PT(R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R )));

end

PT1(m)=sum(PT);

PP(m)=P(m)-PT1(m);

PP1(k+1,m)=PP(m);

end

for m=1:n2

for R=1:n

QT(R)=U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R )));

end

QT1(m)=sum(QT);

QQ(m)=Q(m)-QT1(m);

QQ1(k+1,m)=QQ(m);

end

PR1=max(abs(PP));

PR2=max(abs(QQ));

PR=max(PR1,PR2);

for m=1:n1+n2

for R=1:n1+n2

if m==R

H(m,m)=U1(m)^2*B(m,m)+QT1(m);

else

H(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O( R)));

end

end

end

for m=1:n1+n2

for R=1:n2

if m==R

N(m,m)=-U1(m)^2*G(m,m)-PT1(m);

else

N(m,R)=-U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O (R)));

end

end

end

for m=1:n2

for R=1:n1+n2

if m==R

J(m,m)=U1(m)^2*G(m,m)-PT1(m);

else

J(m,R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O( R)));

end

end

end

for m=1:n2

for R=1:n2

if m==R

L(m,m)=U1(m)^2*B(m,m)-QT1(m);

else

L(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O( R)));

end

end

end

JJ=[H N;J L];

disp(JJ);

PQ=[PP';QQ'];

DA=-inv(JJ)*PQ;

DA1=DA';

for m=1:n1+n2

OO(m)=DA1(m);

end

for m=n:n1+n2+n2

UU1(m-n1-n2)=DA1(m);

end

UD2=diag(UD);

UU=UU1*UD2;

for m=1:n1+n2

O(m)=O(m)+OO(m);

end

for m=1:n2

U1(m)=U1(m)+UU(m);

end

for m=1:n1+n2

o(k+1,m)=180/pi*O(m);

end

for m=1:n2

u(k+1,m)=U1(m);

end

k=k+1;

end

for m=1:n

b(m)=U1(m)*cos(O(m));

c(m)=U1(m)*sin(O(m));

end

U=b+i*c;

for R=1:n

PH1(R)=U(isb)*conj(Y(isb,R))*conj(U(R)); end

PH=sum(PH1);

for m=1:n

for R=1:n

if m~=R

C1(m,R)=1/C(m,R);

else

C1(m,m)=C(m,m);

end

end

end

for m=1:n

for R=1:n

SS(m,R)=U1(m)^2*conj(C1(m,m))+U(m)*(conj(U(m))-conj(U(R)))*conj (C1(m,R));

end

end

SS

disp('迭代中的△P：');disp(PP1);

disp('迭代中的△Q：');disp(QQ1);

disp('迭代中相角:');disp(o);

disp('迭代中电压的模:');disp(u);

disp('平衡结点的功率:');disp(PH);

disp('全部线路功率分布:');disp(SS);