电力系统中matlab纯交流潮流计算

最近刚好要做一个含VSC的交直流潮流计算，先用matlab做一个纯交流潮流计算练练手。照着陈珩第四版的《电力系统分析》编的代码，已经验证过了，结果正确。代码中注释比较全，阅读起来问题不大。分享给大家参考，给一些网友排雷吧！

%%function flow_calculation
clear;clc;
tic
%1节点编号 2电压幅值 3电压相角 4输出有功 5输出无功 6输入有功 7输入无功  8并联电纳 9节点类型
Bus.con=[...
    1  1.030000    0.000000  7.000000  1.850000  0.000000  0.000000  0.000000  2
    2  1.010000    0.000000  7.000000  2.350000  0.000000  0.000000  0.000000  2
    3  1.030000   -6.800000  7.190000  1.760000  0.000000  0.000000  0.000000  1
    4  1.010000    0.000000  7.000000  2.020000  0.000000  0.000000  0.000000  2
    5  1.000000    0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  3
    6  1.000000    0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  3
    7  1.000000    0.000000 -6.000000  0.000000  3.670000  1.000000  2.000000  3
    8  1.000000    0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  3
    9  1.000000    0.000000 -12.00000  0.000000  5.670000  1.000000  3.500000  3
   10  1.000000    0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  3
   11  1.000000    0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  3
   ];

%1首节点 2末节点 3线路电阻 4线路电抗  5线路电纳  6变比
Line.con=[...
   1      5     0.000000  0.016700   0.000000  1
   2      6     0.000000  0.016700   0.000000  1
   3     11     0.000000  0.016700   0.000000  1
   4     10     0.000000  0.016700   0.000000  1  
   5      6     0.002500  0.025000   0.043750  1    
   6      7     0.001000  0.010000   0.017500  1   
   7      8     0.011000  0.110000   0.192500  1    
   7      8     0.011000  0.110000   0.192500  1    
   8      9     0.011000  0.110000   0.192500  1    
   8      9     0.011000  0.110000   0.192500  1 
   9     10     0.001000  0.010000   0.017500  1
  10     11     0.002500  0.025000   0.043750  1
  ];

T=0.00001;

kk=0;%迭代次数

line=size(Line.con,1);%支路数目

n=max(Bus.con(:,1));%找出第1列的最大值，即节点数
TT=ones(n,1)*T;%收敛误差设置

Y=zeros(n);%创建导纳矩阵Y
U=ones(n,1);
dU=zeros(n,1);
U(Bus.con(:,1))=Bus.con(:,2);
th(Bus.con(:,1),1)=(Bus.con(:,3)/180)*pi;
dth=zeros(n,1);
Ps1=Bus.con(:,4)-Bus.con(:,6);
Qs1=Bus.con(:,5)-Bus.con(:,7);

%% 计算交流导纳矩阵
for m=1:line
    a=Line.con(m,1);%首节点a
    b=Line.con(m,2);%末节点b
    z=Line.con(m,3)+1i*Line.