节点牛顿拉夫逊法法matlab程序

节点牛顿拉夫逊法法matlab程序

clear;

clc;

n=9;%节点数；

nl=9;%支路数；

isb=1;%平衡节点号；

pr=0.00001;%误差精度；

b1=[140.0576i01.051;450.017+0.092i0.158i10;560.039+0.17i 0.358i10;360.0586i01.051;670.0119+0.1008i0.209i10;78 0.0085+0.072i0.149i10;280.0625i01.051;890.032+0.161i 0.306i10;940.01+0.085i0.176i10];

%依次是支路首端；末端，支路阻抗；对地电纳；支

比；折算到哪一侧标志（高压侧为1；低压侧为0）；其为支路参数矩阵；

%关于变比为1.05的问题：，=，=1.05，k*=1.05(以全网平均额定电压为基准电压），上述矩阵均是以标幺值给出的

b2=[001.051.0501;1.6301.051.0503;0.8501.051.0503;001 002;00.9+0.3i1002;001002;01+0.35i1002;001002;0 1.25+0.5i1002];

%节点参数矩阵；依次是节点的发电机功率给定值Ps，Qs

（只有2和3节点的功率给定值不为0，分别为1.63+0.067i和0.85-0.109i）；负荷功率给定值；节点电压初值（除发电机节点为1.05外，其它均为1。即一般将PV节点和平衡节点初始电压设为1.05，其它节点初始电压设为1）；PV节点电压Vs给定值（标幺值，除去损耗之后为1，故给定值的标幺值为1.05）；节点无功补偿设备容量；节点分类标号（平衡1；PQ2；PV3）；

Y=zeros(n);%求导纳阵；

for i=1:nl

if b1(i,6)==1

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

else

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

end%为了保证p为低压侧节点，q为高压侧节点

Y(p,q)=Y(p,q)-1./(b1(i,3)*b1(i,5));

Y(q,p)=Y(p,q);

Y(q,q)=Y(q,q)+1./(b1(i,3)*b1(i,5)^2)+b1(i,4)./2;

Y(p,p)=Y(p,p)+1./b1(i,3)+b1(i,4)./2;

end

%disp('系统的导纳阵为：');

%disp(Y);

%求解导纳矩阵；

G=real(Y);B=imag(Y);%取导纳矩阵的实部和虚部；

for i=1:n

e(i)=real(b2(i,3));（1）

f(i)=imag(b2(i,3));（2）

%（1）和（2）为PQ节点电压的实部和虚部

v(i)=b2(i,4);%PV节点电压

end

%电压赋初值；

for i=1:n

S(i)=b2(i,1)-b2(i,2);

B(i,i)=B(i,i)+b2(i,5);%加入无功补偿设备的导纳；

end

P=real(S);Q=imag(S);%给有功功率和无功功率赋初值

w=zeros(2*n,1);%定义功率和电压变化量组成2n*1的矩阵；Jac=zeros(2*n);%定义雅可比矩阵；

t=0;

while t==0

for i=1:n

if b2(i,6)~=isb

C=0;D=0;

for j=1:n

C=C+G(i,j)*e(j)-B(i,j)*f(j);

D=D+G(i,j)*f(j)+B(i,j)*e(j);

end

if b2(i,6)==2%P，Q节点；

w(2*i)=P(i)-e(i)*C-f(i)*D;

w(2*i-1)=Q(i)-f(i)*C+e(i)*D;

else if b2(i,6)==3%P，V节点；

w(2*i)=P(i)-e(i)*C-f(i)*D;

w(2*i-1)=v(i)^2-(e(i)^2+f(i)^2);

end

end

else

w(2*i-1)=0;

w(2*i)=0;%平衡节点不参与求解雅可比矩阵

end

end

%disp(w);%求解有功无功及电压变化量的矩阵；

w1=w(3:2*n);%除去平衡节点

for

i=1:n

for j=1:n

if b2(i,6)~=isb

if b2(i,6)==2%P，Q节点；

if j~=i

Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i)); Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i));

Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);

Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i); else if j==i

m=0;h=0;

for r=1:n

m=m+G(i,r)*e(r)-B(i,r)*f(r);

h=h+G(i,r)*f(r)+B(i,r)*e(r);

end

Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i); Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i); Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i); Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i); end

end

else if b2(i,6)==3%P，V节点；

if j~=i

Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i)); Jac(2*i-1,2*j)=0;

Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);

Jac(2*i-1,2*j-1)=0;

else if j==i

m=0;h=0;

for r=1:n

m=m+G(i,r)*e(r)-B(i,r)*f(r);

h=h+G(i,r)*f(r)+B(i,r)*e(r);

end

Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i); Jac(2*i-1,2*j)=-2*f(i);

Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i); Jac(2*i-1,2*j-1)=-2*e(i);

end

end

end

end

else

Jac(2*i-1,2*j-1)=0;

Jac(2*i,2*j)=0;

Jac(2*i-1,2*j)=0;

Jac(2*i,2*j-1)=0;