潮流计算范例

%以下部分为输入原始数据（到标示‘///’标志为止）。 g(5,1)=0.000; b(5,1)=-31.746; g(5,2)=0.000; b(5,2)=0.000; g(1,2)=0.830; b(1,2)=-3.112;

g(1,3)=0.755; b(1,3)=-2.642; g(2,3)=0.624; b(2,3)=-3.900; g(3,4)=0.000; b(3,4)=0.000; g(1,4)=0.000; b(1,4)=0.000; g(2,4)=0.000; b(2,4)=-63.492; g(1,5)=0.000; b(1,5)=-31.746; g(2,5)=0.000; b(2,5)=0.000; g(2,1)=0.830; b(2,1)=-3.112; g(3,1)=0.755; b(3,1)=-2.642; g(3,2)=0.624; b(3,2)=-3.900; g(4,3)=0.000; b(4,3)=0.000; g(4,2)=0.000; b(4,2)=-63.492; g(4,1)=0.000; b(4,1)=0.000; g(1,1)=0.000; b(1,1)=1.762; g(2,2)=0.000; b(2,2)=3.524; g(3,3)=0.000; b(3,3)=0.250; g(4,4)=0.000; b(4,4)=-3.175; g(5,5)=0.000; b(5,5)=-1.587;

d(5,1)=-j*1.587;d(1,5)=j*1.512;d(4,2)=-j*3.175; d(2,4)=j*3.023;d(3,2)=j*0.250;d(2,3)=j*0.250; d(1,2)=j*0.250;d(2,1)=j*0.250;

%求取节点导纳矩阵。 for m=1:5

for n=1:5 if m==n

G(m,m)=g(m,1)+g(m,2)+g(m,3)+g(m,4)+g(m,5);

B(m,m)=b(m,1)+b(m,2)+b(m,3)+b(m,4)+b(m,5); else

G(m,n)=-g(m,n); B(m,n)=-b(m,n); end end end

Y=G+j*B;

%//////////////////////////////////////////////////////////////

%//设定节点起始节点电压，并给出已知功率值。 %//下面将题中节点 2、3、4、5、1分别替换为节点1、2、3、4、5即节点4为PV节点，节点5为平衡节点。 delt(1)=0;delt(2)=0;delt(3)=0;delt(4)=0; u(1)=1.0;u(2)=1.0;u(3)=1.0;

p(1)=-3.7;q(1)=-1.3;p(2)=-2;q(2)=-1; p(3)=-1.6;q(3)=-0.8;p(4)=5;q(4)=0.00;

%置迭代次数k=0，并计算节点功率的不平衡量 ， 。 k=0;precision=1; k,delt,u

N1=4;

while precision>0.00001

u(4)=1.05; delt(5)=0;u(5)=1.05; for m=1:N1

for n=1:N1+1

pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n))); end

pp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt); end

%///////////////////////////////////////////////////////

%求取雅可比矩阵元素 （m=n 时）。 for m=1:N1

for n=1:N1+1

h0(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));

n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

l0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n))); end

H(m,m)=sum(h0)-u(m)^2*(G(m,m)*sin(delt(m)-delt(m

))-B(m,m)*cos(delt(m)-delt(m)));

N(m,m)=sum(n0)-2*u(m)^2*G(m,m)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))

+B(m,m)*sin(delt(m)-delt(m)));

J(m,m)=sum(j0)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));

L(m,m)=sum(l0)+2*u(m)^2*B(m,m)+u(m)^2*(G(m,m)*sin(delt(m)-delt(m))

-B(m,m)*cos(delt(m)-delt(m))); end

for m=1:N1-1

JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m); JJ(2*m,2*m-1)=J(m,m);JJ(2*m,2*m)=L(m,m); end

for m=N1:N1

JJ(2*m-1,2*m-1)=H(m,m); end

%///////////////////////////////////////////////////////

%求取雅可比矩阵元素 （m不等于n时） for m=1:N1

for n=1:N1 if m==n else

H(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B

(m,n)*cos(delt(m)-delt(n)));

J(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n))); N(m,n)=-J(m,n);L(m,n)=H(m,n); end end end

for m=1:N1-1 %求前六行六列元素(非对角) for n=1:N1-1 if m==n else

JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n); JJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n); end end end

for m=N1

for n=1:N1-1 %求取第七行的元素

JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n); end end

for n=N1

for m=1:N1-1 %求的第七列元素

JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m,2*n-1)=J(m,n); end end

%///////////////////////////////////////////////

////////

%解修正方程式，由 ， 和 计算电压修正量 和 。 for m=1:N1-1

PP(2*m-1)=pp(m);PP(2*m)=qq(m); end

for m=N1

PP(2*m-1)=pp(m); end

uu=-inv(JJ)*PP';precision=max(abs(uu));

%///////////////////////////////////////////////////////

%若结果不收敛，执行下列语句。 for n=1:N1-1

delt(n)=delt(n)+uu(2*n-1); u(n)=u(n)+uu(2*n); end

for n=N1

delt(n)=delt(n)+uu(2*n-1); end

k=k+1; k,delt,u end

%///////////////////////////////////////////////////////

%若结果收敛，计算各节点电压，平衡节点功率，PV节点功率和线路功率 。

for n=1:N1+1

U(n)=u(n)*(cos(delt(n))+j*sin(delt(n))); end

for m=1:N1+1

I(m)=Y(5,m)*U(m); end

S5=U(5)*sum(conj(I)) %平衡节点功率 %PV节点功率 for n=1:N1+1

q4(n)=u(4)*u(n)*(G(4,n)*sin(delt(4)-delt(n))-B(4,n)*cos(delt(4)-delt(n))); end

Q4=sum(q4) %线路功率 for m=1:N1+1

for n=1:N1+1

S(m,n)=U(m)*(conj(U(m))*conj(d(m,n))+(conj(U(m))-conj(U(n)))*conj(-Y(m,n)));

end end

%///////////////////////////////////////////////////////

%显示运行结果（至结束）。 Y

JJ S B pp qq uu U k

PQ分解法潮流计算程序

G(1,1)=0.004817;B(1,1)=-0.016114;

G(1,2)=-0.004569;B(1,2)=0.007030;G(1,3)=0;

B(1,3)=0;G(1,4)=0;B(1,4)=0;G(1,5)=-0.000248;B(1,5)=0.009084;

G(2,1)=-0.004569;B(2,1)=0.007030;

G(2,2)=0.009138;B(2,2)=-0.01406;G(2,3)=-0.004569;

B(2,3)=0.007030;G(2,4)=0;B(2,4)=0;G(2,5)=0;B(2,5)=0;

G(3,1)=0;B(3,1)=0;

G(3,2)=-0.004569;B(3,2)=0.007030;G(3,3)=0.006079

;B(3,3)=-0.050450;

G(3,4)=-0.001510;B(3,4)=0.043420;G(3,5)=0;B(3,5)=0;

B(4,1)=0;G(4,2)=0;B(4,2)=0;G(4,3)=-0.001510;B(4,3)=0.043420;

G(4,4)=0.007255;B(4,4)=-0.074090;G(4,5)=-0.005745;B(4,5)=0.030670;G(4,1)=0;

G(5,1)=-0.000248;B(5,1)=0.009084;G(5,2)=0;B(5,2)=0;

G(5,3)=0;B(5,3)=0;G(5,4)=-0.005745;B(5,4)=0.030670;G(5,5)=0.005993;B(5,5)=-0.039754; Y=G+j*B

B1=B(1:4,1:4); B11=B1;

delt(1)=0; delt(2)=0; delt(3)=0; delt(4)=0;

u(1)=220; u(2)=220; u(3)=220; u(4)=220;

p(1)=40; q(1)=30; p(2)=0; q(2)=10;

p(3)=-180; q(3)=-100; p(4)=-50; q(4)=-30; k=0;precision=1; k,delt,u N1=4;

while precision>0.00001

delt(5)=0;u(5)=242; for m=1:N1

for n=1:N1+1

pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n))); end

pp(m)=p(m)-sum(pt); qq(m)=q(m)-sum(qt); end

for m=1:N1

PU(m)=pp(m)/u(m);

PU1(m)=pp(m)/(u(m)^2); QU(m)=qq(m)/u(m); end

ddelt=-inv(B1)*PU1';

precision=max(abs(PU)); for n=1:N1

delt(n)=delt(n)+ddelt(n); end

uu=-inv(B11)*QU';

precision=max(abs(QU));

for n=1:N1

u(n)=u(n)+uu(n); end

k=k+1;

k,delt,u end

for n=1:N1+1

U(n)=u(n)*(cos(delt(n))+j*sin(delt(n))); end

for m=1:N1+1

I(m)=Y(5,m)*U(m); end

S5=U(5)*sum(conj(I)) for m=1:N1+1

for n=1:N1+1

S(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n));

end end S

x=1:k;

subplot(4,1,1); plot(x,u(1),'rp') grid on

title('迭代过程中各节点的电压幅值') subplot(4,1,2); plot(x,u(2),'rp') grid on

subplot(4,1,3); plot(x,u(3),'rp') grid on

ylabel('电压标幺值') subplot(4,1,4); plot(x,u(4),'rp') grid on

xlabel('迭代次数')

x=1:k;

subplot(4,1,1);

plot(x,delt(1),'rp') grid on

title('迭代过程中各节点的电压幅值') subplot(4,1,2); plot(x,u(2),'rp') grid on

subplot(4,1,3); plot(x,u(3),'rp') grid on

ylabel('')

subplot(4,1,4); plot(x,u(4),'rp')

grid on

xlabel('迭代次数')

clc; clear;

g(5,1)=5.0000; b(5,1)=-15.0000; g(5,2)=1.2500; b(5,2)=-3.7500; g(1,2)=1.6667; b(1,2)=-5.0000; g(1,3)=3.333; b(1,3)=-6.6667;

g(2,3)=10.0000; b(2,3)=-30.0000; g(3,4)=1.2500;b(3,4)=-3.7500; g(1,4)=2.7500; b(1,4)=-8.2500;

g(1,5)=5.0000; b(1,5)=-15.0000; g(2,5)=1.2500; b(2,5)=-3.7500; g(2,1)=1.6667; b(2,1)=-5.0000; g(3,1)=3.333; b(3,1)=-6.6667;

g(3,2)=10.0000; b(3,2)=-30.0000; g(4,3)=1.2500;b(4,3)=-3.7500; g(4,1)=2.7500; b(4,1)=-8.2500; g(1,1)=0.2750;

b(1,1)=-0.8250; %g(1,6)?ú±íg(1,0)??b(1,6)?ú±íb?¨1,0??

g(4,4)=-0.2500;

b(4,4)=0.7500; %g(4,6)?ú±íg?¨4,0????b?¨4??6???ú±íb?¨4,0??

%?ó?????????????ó?? for m=1:5

for n=1:5 if m==n

G(m,m)=g(m,1)+g(m,2)+g(m,3)+g(m,4)+g(m,5);

B(m,m)=b(m,1)+b(m,2)+b(m,3)+b(m,4)+b(m,5); else

G(m,n)=-g(m,n); B(m,n)=-b(m,n); end end end

Y=G+j*B

%////////////////////////////////////////////////////////////// %//?è?¨???????????????????????????????????? %//?????????????? 2??3??4??5??1??±???????????1??2??3??4??5??????4??PV??????????5????????????

e(1)=1.0;e(2)=1.0;e(3)=1.0;e(4)=1.10; e(5)=1.06;f(5)=0.00;

f(1)=0;f(2)=0;f(3)=0;f(4)=0;

u(1)=1.0;u(2)=1.0;u(3)=1.0;u(4)=1.0; uu(1)=0;uu(2)=0;uu(3)=0;uu(4)=0;

p(1)=0.20;q(1)=0.20;p(2)=-0.45;q(2)=-0.15;

p(3)=-0.40;q(3)=-0.05;p(4)=-0.50;q(4)=0.00; %???ü?ú????k=0??????????????????????????. k=0;precision=1; k,e,f,uu N1=4;

while precision>0.00001

for m=1:N1-1

for n=1:N1+1

pt(n)=e(m)*(G(m,n)*e(n)-B(m,n)*f(n))+f(m)*(G(m,n)*f(n)+B(m,n)*e(n));

qt(n)=f(m)*(G(m,n)*e(n)-B(m,n)*f(n))-e(m)*(G(m,n)*f(n)+B(m,n)*e(n)); end

pp(m)=p(m)-sum(pt); qq(m)=q(m)-sum(qt); end

for m=N1

for n=1:N1+1

pt(n)=e(m)*(G(m,n)*e(n)-B(m,n)*f(n))+f(m)*(G(m,n)*f(n)+B(m,n)*e(n)); end

pp(m)=p(m)-sum(pt);

uu(m)=u(m)*u(m)-(e(m)*e(m)+f(m)*f(m)); end

for m=1:N1-1

for n=1:N1+1

h0(n)=G(m,n)*f(n)+B(m,n)*e(n);

n0(n)=G(m,n)*e(n)-B(m,n)*f(n); j0(n)=n0(n); l0(n)=h0(n); end

H(m,m)=-B(m,m)*e(m)+G(m,m)*f(m)+sum(h0); N(m,m)=G(m,m)*e(m)+B(m,m)*f(n)+sum(n0); J(m,m)=-G(m,m)*e(m)-B(m,m)*f(m)+sum(j0);

L(m,m)=-B(m,m)*e(m)+G(m,m)*f(m)-sum(l0); end

for m=N1

for n=1:N1+1

h0(n)=G(m,n)*f(n)+B(m,n)*e(n); n0(n)=G(m,n)*e(n)-B(m,n)*f(n); end

H(m,m)=-B(m,m)*e(m)+G(m,m)*f(m)+sum(h0); N(m,m)=G(m,m)*e(m)+B(m,m)*f(n)+sum(n0); R(m,m)=2*f(m);

S(m,m)=2*e(m); end

for m=1:3

JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m); JJ(2*m,2*m-1)=J(m,m);JJ(2*m,2*m)=L(m,m); end

m=4;

JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m); JJ(2*m,2*m-1)=R(m,m);JJ(2*m,2*m)=S(m,m);

for m=1:N1-1

for n=1:N1+1 if m==n else

H(m,n)=-B(m,n)*e(m)+G(m,n)*f(m); N(m,n)=G(m,n)*e(m)+B(m,n)*f(m); J(m,n)=-N(m,n);L(m,n)=H(m,n); end end end

for m=N1

for n=1:N1+1 if m==n else

H(m,n)=-B(m,n)*e(m)+G(m,n)*f(m); N(m,n)=G(m,n)*e(m)+B(m,n)*f(m); S(m,n)=0;R(m,n)=0 end end end

for m=1:N1-1 for n=1:N1 if m==n else

JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n); JJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n); end

end end

for m=N1

for n=1:N1

JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n); JJ(2*m,2*n-1)=R(m,n);JJ(2*m,2*n)=S(m,n); end end

for m=1:N1-1

PP(2*m-1)=pp(m);PP(2*m)=qq(m); end

for m=N1

PP(2*m-1)=pp(m);PP(2*m)=uu(m); end

uu=inv(JJ)*PP';

precision =max(abs(uu)); for n=1:N1

f(n)=f(n)+uu(2*n-1); e(n)=e(n)+uu(2*n); end

k=k+1; k,e,f,u end

for n=1:N1+1

U(n)=e(n)+j*f(n); end

for m=1:N1+1

I(m)=Y(5,m)*U(m); end

S5=U(5)*sum(conj(I)) for m=1:N1+1 for n=1:N1+1

S(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n)+conj(U(m))*conj(g(m,m)+j*b(m,m))); end end Y JJ S B pp qq uu