牛顿法计算潮流的matlab程序

function [bus_resS_res] = PowerFlow_NR_2          % 牛顿-拉夫逊法解潮流方程的主程序
                         
bus = [   1    1  0       0       0     3;
          2    1  0       -0.4    -0.3   1;
          3    1  0       -0.5    -0.3   1;
          4    1  0      -0.147  -0.27  1;
          5    1  0      -0.344  -0.17  1;
          6    1  0      -0.8    -0.36  1;      
          7    1  0      -0.3    -0.12  1;
          8    1  0      -0.6    -0.3   1;
          9    1  0      -0.18   -0.4   1;
         10    1  0       1    0     2  ];
line=[  1       2	0.245	0.399	0	0	25;
        2	3	0.233	0.379	0	0	1;
        3	4	0.196	0.319	0	0	1;
        4	5	0.196	0.319	0	0	1;
        5	6	0.221	0.359   0	0	1;
        5	9	0.453	0.738	0	0	1;
        4	6	0.392	0.638	0	0	1;
        6	3	0.576	0.938	0	0	1;
        6	7	0.343	0.559	0	0	1;
        7	8	0.294	0.479	0	0	1;
        9	8	0.282	0.459	0	0	1;
        7	1	0.882	1.436	0	0	0.04;
        9	1	0.380	0.619	0	0	0.04];

[nbmb]=size（bus）;
[nlml]=size（line）;                                % 计算bus和line矩阵的行数和列数
[buslinenPQnPVnodenum] = Num_（busline）;       % 对节点重新排序的子程序
Y = y_（busline）;                                  % 计算节点导纳矩阵的子程序
myf = fopen（‘Result.m‘‘w‘）;
fprintf（myf‘-- by tongchao ----\n\n‘）;
fclose（myf）;                                       % 在当前目录下生成“Result.m”文件，写入节点导纳矩阵
format long
EPS = 1.0e-5;                                     % 设定误差精度
for t = 1:100                                      % 开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出
 
    [dPdQ] = dPQ_（YbusnPQnPV）;                 % 计算功率偏差dP和dQ的子程序
    J = Jac_（busYnPQ）;                           % 计算雅克比矩阵的子程序
    UD = zeros（nPQnPQ）;
 
    for i = 1:nPQ
        UD（ii） = bus（i2）;                        % 生成电压对角矩阵
    end
    dAngU = J \ [dP;dQ];
    dAng = dAngU（1:nb-11）;                        % 计算相角修正量
    dU = UD*（dAngU（nb:nb+nPQ-11））;                % 计算电压修正量
 
    bus（1:nPQ2） = bus（1:nPQ2） - dU;              % 修正电压
    bus（1:nb-13） = bus（1:nb-13） - dAng;          % 修正相角
  
    if （max（abs（dU））        break
    end                                            % 判断是否满足精度误差，如满足则跳出，否则返回继续迭代
 
end
bus = PQ_（busYnPQnPV）;                          % 计算每个节点的有功和无功注入的子程序
 
[busline] = ReNum_（buslinenodenum）;             % 对节点恢复编号的子程序
YtYm = YtYm_（line）;                                % 计算线路的等效Yt和Ym的子程序，以计算线路潮流
bus_res = bus_res_（bus）;                           % 计算节点数据结果的子程序
S_res = S_res_（buslineYtYm）;                     % 计算线路潮流的子程序
myf = fopen（‘Result.m‘‘a‘）;
fprintf（myf‘----牛顿－拉夫逊法潮流计算结果------\n\n 节点计算结果：\n节点     节点电压     节点相角（角度）     节点注入功率\n‘）;
for i = 1:nb
    fprintf（myf‘%f  ‘bus_res（i1））;
    fprintf（myf‘%f ‘bus_res（i2））;
    fprintf（myf‘%f ‘bus_res（i3））;
    fprintf（myf‘%f + j %f\n‘real（bus_res（i4））imag（bus_res（i4）））;
end
    fprintf（myf‘\n线路计算结果：\n节点I    节点J          线路功率S（IJ）               线路功率S（JI）               线路损耗dS（IJ）\n‘）