基于MATLAB的牛顿拉夫逊迭代法计算潮流(附加短路计算).doc

这个程序可以适用于三机九节点系统（参数见主程序），本来是想编一个通用各种结构的程序的，但是因为鄙人比较懒，老师留作业时候又没说要通用，就没改完。。。惭愧啊。。。不过大同小异啦。。。有兴趣的慢慢改吧。

使用方法：
按后文中给出的代码建立.m文件放在一个文件夹里面。
先运行Powerflow_main.m计算算例系统的潮流；
然后运行ShortcircuitCalc.m计算算例系统三相短路电流；
程序说明详见各.m文件注释部分，写的已经很详细了，慢慢看吧。

Powerflow_main.m文件代码如下：
clear
%牛顿拉夫逊迭代法计算潮流
format short %规定参数数据显示精度
%节点参数矩阵
%第一列为节点编号
%第二列表示有功注入P
%第三列表示无功注入Q
%第四列表示电压幅值U
%第五列表示电压角度θ
%第六列表示发电机x′
%第七列表示发电机E′
%第八列表示节点类型（2表示平衡节点，1表示PV节点，0表示PQ节点）
Node_p=[ 1, 0, 0, 1.04 , 0, 0.3, 1.137, 2;
2, 1.63, 0, 1.025, 0, 0.3, 1.211, 1;
3, 0.85, 0, 1.025, 0, 0.3, 1.047, 1;
4, 0, 0, 1.0, 0, 0, 0, 0;
5, -1.25, -0.5, 1.0, 0, 0, 0, 0;
6, -0.9, -0.3, 1.0, 0, 0, 0, 0;
7, 0, 0, 1.0, 0, 0, 0, 0;
8, -1, -0.35, 1.0, 0, 0, 0, 0;
9, 0, 0, 1.0, 0, 0, 0, 0
];

count_s=0;
countPV=0;
for k=1:size(Node_p,1)
if Node_p(k,8)==1
countPV=countPV+1;
else if Node_p(k,8)==2
count_s=count_s+1;
end;
end;
end;
countPV;
count_s;
countPQ=size(Node_p,1)-1-countPV;

%显示节点参数
disp('节点参数如下：')
disp(Node_p)
%支路参数
%第一列为首节点，第二列为末节点，第三列表示R，第四列表示X，第五列表示B/2
%第六列表示支路类型（1为变比为1的变压器元件；2为输电线元件;0为接地支路）
Branch_p =[ 1, 4, 0 , 0.0576, 0 , 1;
2, 7, 0 , 0.0625, 0 , 1;
3, 9, 0 , 0.0586, 0 , 1;
4, 5, 0.01 , 0.085 , 0.088 , 2;
4, 6, 0.017 , 0.092 , 0.079 , 2;
5, 7, 0.032 , 0.161 , 0.153 , 2;
6, 9, 0.039 , 0.17 , 0.179 , 2;
7, 8, 0.0085, 0.072 , 0.0745, 2;
8, 9, 0.0119, 0.1008, 0.1045, 2
];
%显示支路参数
disp('支路参数如下：')
disp(Branch_p)

%设置节点初值
U=Node_p(:,4);
e_ang=Node_p(:,5);

P=Node_p(:,2);
Q=Node_p(:,3);
save data.mat

format long%计算结果数据显示精度
%显示节点导纳矩阵
admi();
disp('节点导纳矩阵Y');
sparseY

Kmax=10; %设置最大迭代次数k
accuracy=10^-7;%设置迭代精度
k=0;%迭代次数初始化为零
for k1=1:Kmax
[dP,dQ,y]=getY(U,e_ang);
if max(abs(y))<accuracy
break;
end;
J=jacob(U,e_ang,dP,dQ);
x=-inv(J)*y;

de_ang=[0;x(1:8)];
dU=x(9:14);
u1=U(2:9)*dU.';
u=diag(u1);
U(4:9)=U(4:9)+u;
e_ang=e_ang+de_ang;
k=k+1;
end;

e_ang=e_ang/pi*180;
save result.mat

disp('迭代次数：')
k

disp('4号节点至9号节点电压幅值如下：')
disp(U(4:9))
disp('2号节点至9号节点电压相角如下：')
disp(e_ang(2:9))

admi.m文件代码如下：
%节点导纳矩阵的形成
format long %规定数据格式
NI=size(Branch_p,1);
k_t=1;
Y=zeros(NI);
for m1=1:NI;
I=Branch_p(m1,1);
J=Branch_p(m1,2);
R=Branch_p(m1,3);
X=Branch_p(m1,4);
b=Branch_p(m1,5);
Style=Branch_p(m1,6);
if Style==1 %判断为变压器元件
Y(I,I)=Y(I,I)+1/(R+1j*X);
Y(J,J)=Y(J,J)+1/(R+1j*X)/k_t/k_t;
Y(I,J)=Y(I,J)-1/(R+1j*X)/k_t;
Y(J,I)=Y(J,I)-1/(R+1j*X)/k_t;
else if Style==0 %判断为母线接地支路元件
Y(I,J)=Y(I,J)+1/(R+1j*X);

else %判断为输电线元件
Y(I,I)=Y(I,I)+1j*b+1/(R+1j*X);
Y(J,J)=Y(J,J)+1j*b+1/(R+1j*X);
Y(I,J)=Y(I,J)-1/(R+1j*X);
Y(J,I)=Y(J,I)-1/(R+1j*X);
end;
end;
end;
Y;
G=real(Y);
B=imag(Y);
sparseY=sparse(Y);
save data.mat;

getY.m文件代码如下：（线性方程组常写作AX=Y形式，故此处命名为get_Y）
%计算ΔP和ΔQ的函数
function [dP,dQ,y]=getY(U,e_ang)
load data P Q G B Node_p countPV count_s;
dP=zeros(size(Node_p,1),1); sum1=zeros(size(Node_p,1),1);
for i=1:size(Node_p,1)
for j=1:size(Node_p,1)
sum1(i)=sum1(i)+U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j))+B(i,j)*sin(e_ang(i)-e_ang(j)));
end;
dP(i)=P(i)-U(i)*sum1(i);
end;
dQ=zeros(size(Node_p,1),1);sum2=zeros(size(Node_p,1),1);
for i=1:size(Node_p,1)
for j=1:size(Node_p,1)
sum2(i)=sum2(i)+U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j))-B(i,j)*cos(e_ang(i)-e_ang(j)));
end;
dQ(i)=Q(i)-U(i)*sum2(i);
end;

y=[dP((count_s+1):9);dQ((countPV+count_s+1):9)]; %拼接ΔP和ΔQ构成方程 -J*x=y的向量y

jacob.m文件的代码如下：
%形成雅克比矩阵的函数
function J=jacob(U,e_ang,dP,dQ)
load data B G Q P Node_p countPV count_s;

size_Y=size(Node_p,1);
%H矩阵
H1=zeros(size_Y);
for i=1:size_Y
for j=1:size_Y
if j==i;
H1(i,j)= U(i)*U(i)*B(i,j)+Q(i)-dQ(i);
else
H1(i,j)= -U(i)*U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j))-B(i,j)*cos(e_ang(i)-e_ang(j)));
end;
end
end
H=H1(2:size_Y,2:size_Y);
%N矩阵
N1=zeros(size_Y);
for i=1:size_Y
for j=1:size_Y
if j==i;
N1(i,j)= -U(i)*U(i)*G(i,j)-(P(i)-dP(i));
else
N1(i,j)= -U(i)*U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j))+B(i,j)*sin(e_ang(i)-e_ang(j)));
end;
end;
end;
N=N1(2:size_Y,(countPV+count_s+1):size_Y);
%M矩阵
M1=zeros(size_Y);
for i=1:size_Y
for j=1:size_Y
if j==i;
M1(i,j)=U(i)*U(i)*G(i,j)-(P(j)-dP(j));
else
M1(i,j)=U(i)*U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j))+B(i,j)*sin(e_ang(i)-e_ang(j)));
end;
end;
end;
M=M1((countPV+count_s+1):size_Y,2:size_Y);
%L矩阵
L1=zeros(size_Y);
for i=1:size_Y
for j=1:size_Y
if j==i;
L1(i,j)= U(i)*U(i)*B(i,j)-(Q(i)-dQ(i));
else
L1(i,j)= -U(i)*U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j))-B(i,j)*cos(e_ang(i)-e_ang(j)));
end;
end
end
L=L1((countPV+count_s+1):size_Y,(countPV+count_s+1):size_Y);

J=[H N;M L];%拼接构成雅克比矩阵

ShortcircuitCalc.m文件的代码如下：
clear
%三相短路计算
format long;
%对YN进行修正，形成包括发电机内阻抗和负荷阻抗的节点导纳矩阵
re_admi();
Z=inv(rY);%计算节点阻抗矩阵
disp('4节点发生金属短路')
f=4;%短路点为4节点
%输出节点阻抗矩阵的短路点所在列
disp('节点阻抗矩阵的短路点所在列Z(:,f)=');
Z(:,f)

%短路电流If计算
load result U e_ang;
e_angf=e_ang(f);
Uf=U(f)*(cos(e_angf/180*pi)+1j*sin(e_angf/180*pi));
Zff=Z(f,f);
zf=0;
If=Uf/(Zff+zf);
i_angf=angle(If)*180/pi;
If=abs(If);
%短路电流计算结果显示
disp('短路电流幅值:')
If
disp('短路电流相角(单位为 °)')
i_angf

%短路时各节点电压计算
U_k=U;
for x1=1:size(Node_p,1)
U_k(x1)=U(x1)*(cos(e_ang(x1)*pi/180)+1j*sin(e_ang(x1)*pi/180))-Z(x1,f)*If*(cos(i_angf*pi/180)+1j*sin(i_angf*pi/180));
end;
%Uk为短路时各节点电压幅值
Uk=abs(U_k);
%uk_ang为短路时各节点电压相角（单位为 °）
uk_ang=angle(U_k)*180/pi;
%输出短路时各节点电压
disp('短路时1-9节点电压：')
disp('幅值：')
Uk
disp('相角：(单位为 °)')
uk_ang

%计算短路时各支路电流
%输出矩阵初始化
%第一列为首节点i
%第二列为末节点j
%第三列为Ii幅值
%第四列为Ii相角(单位为 °)
%第五列为Ij幅值
%第六列为Ij相角（单位为 °）
Ik=[1 4 1 0 1 0;
2 7 1 0 1 0;
3 9 1 0 1 0;
4 5 1 0 1 0;
4 6 1 0 1 0;
5 7 1 0 1 0;
6 9 1 0 1 0;
7 8 1 0 1 0;
8 9 1 0 1 0];
load data Branch_p k_t;%其中变压器变比为k_t
NI=size(Branch_p,1);
for m1=1:NI;
I=Branch_p(m1,1);
J=Branch_p(m1,2);
R=Branch_p(m1,3);
X=Branch_p(m1,4);
b=Branch_p(m1,5);
Style=Branch_p(m1,6);
if Style~=1 %判断为非变压器支路
Ik(m1,3)=(U_k(I)-U_k(J))/(R+1j*X)+U_k(I)*1j*b;
Ik(m1,4)=angle(Ik(m1,3))*180/pi;
Ik(m1,3)=abs(Ik(m1,3));
Ik(m1,5)=(U_k(J)-U_k(I))/(R+1j*X)+U_k(J)*1j*b;
Ik(m1,6)=angle(Ik(m1,5))*180/pi;
Ik(m1,5)=abs(Ik(m1,5));
else %否则为变压器支路
Ik(m1,3)=(U_k(I)-U_k(J))/(R+1j*X)/k_t;
Ik(m1,4)=angle(Ik(m1,3))*180/pi;
Ik(m1,3)=abs(Ik(m1,3));
Ik(m1,5)=(U_k(J)-U_k(I))/(R+1j*X)/k_t;
Ik(m1,6)=angle(Ik(m1,5))*180/pi;
Ik(m1,5)=abs(Ik(m1,5));
end
end
%输出短路时各节点电流
disp('第一列为首节点i')
disp('第二列为末节点j')
disp('第三列为Ii幅值')
disp('第四列为Ii相角(单位为 °)')
disp('第五列为Ij幅值')
disp('第六列为Ij相角（单位为 °）')
format short
Ik
save result_k

%形成包括发电机内阻抗、负荷阻抗的节点导纳矩阵
format long %规定数据显示精度
load data P Q Y Node_p ;
load result U;
rY=Y;
for n1=1:size(Node_p,1)
Xd=Node_p(n1,6);
Ed=Node_p(n1,7);
style=Node_p(n1,8);
if style~=0 %判断是否为发电机节点
rY(n1,n1)=rY(n1,n1)+1/(1j*Xd);%YN中与发电机节点对应的对角线元素增加发电机导纳
else %其他节点(包含负荷节点）
rY(n1,n1)=rY(n1,n1)+(-P(n1)+1j*Q(n1))/(U(n1)*U(n1));%与负荷节点对应的对角线元素增加负荷导纳
%对于非负荷节点 Y矩阵元素不做修改 但仍满足上式
end;
end;
rY;