节点牛顿拉夫逊法法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;
end
end
end
%disp(Jac);%求解包含平衡节点在内的雅可比矩阵；
Jac1=Jac(3:2*n,3:2*n);%平衡节点不计入雅可比矩阵内
for k=1:2*n-2
m=0;
for i=k+1:2*n-2
m=Jac1(i,k)./Jac1(k,k);
w1(i)=w1(i)-m*w1(k);%对应于线性方程的增光矩阵做初等变换
for j=k+1:2*n-2
Jac1(i,j)=Jac1(i,j)-m*Jac1(k,j);
%通过矩阵Jac1的第i行-m*对应第一行的元素，从而通过矩阵的初等变换将第i行的第一个元素是化为0
end
end
end%将雅可比矩阵Jac1转变为上三角矩阵（所有0在下面）
E(2*n-2)=w1(2*n-2)./Jac1(2*n-2,2*n-2);
for i=2*n-3:-1:1
%i=s:d:f，s表示起始值，d表示增量，f表示终点值（这里2*n-3为Jac1的最后一行或一列;Jac1为(2n-2）*(2n-2）阶矩阵；Jac为2n*2n阶矩阵)
c=0;
for k=i+1:2*n-2
c=c+Jac1(i,k)*E(k);
E(i)=(w1(i)-c)./Jac1(i,i);
end
end
%disp(E);%求得电压实部和虚部的修正量（△V）
for i=1
:n-1
e(i+1)=e(i+1)-E(2*i-1);
f(i+1)=f(i+1)-E(2*i);%（1）这里，故减号。

（2）由于e(i)和f(i)的第一个元素包含了平衡节点，而在计算E(i)中不包含平衡节点，所以为i+1，即从2开始。

end
%disp(e);
%disp(f);
for i=1:n-1
b(2*i-1)=abs(E(2*i-1));
b(2*i)=abs(E(2*i));
end
KB=max(b);
%disp(KB);
if KB<pr
t=1;
else
t=0;
end
end
%disp(e);
%disp(f);
for i=1:n
fz(i)=sqrt(e(i)^2+f(i)^2);
end
disp(fz);%牛顿拉夫逊法求解正常状态下的潮流，这里只有电压的幅值。