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

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

使用方法：

按后文中给出的代码建立.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))

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 %判断为变压器元件