分层前推回代潮流计算

function [CL]=powerflow（）
b=32;
n=33;
LL=5;           %联络开关数  
Sb=10;          %MW
Vb=12.66;       %KV
Zb=Vb^2/Sb;     %阻抗
check=1;
checkhl=1;
checkgd=1;
jy=0;
BranchM=[1 1 2    0.0922+i*0.047
    2 2 3    0.4930+i*0.2511
    3 3 4    0.3660+i*0.1864
    4 4 5    0.3811+i*0.1941
    5 5 6    0.8190+i*0.7070
    6 6 7    0.1872+i*0.6188
    7 7 8    0.7114+i*0.2351
    8 8 9    1.0300+i*0.7400
    9 9 10   1.0440+i*0.7400
    10 10 11 0.1966+i*0.0650
    11 11 12 0.3744+i*0.1238
    12 12 13 1.4680+i*1.1550
    13 13 14 0.5416+i*0.7129
    14 14 15 0.5910+i*0.5260
    15 15 16 0.7463+i*0.5450
    16 16 17 1.2890+i*1.7210
    17 17 18 0.3720+i*0.5740
    18 2 19 0.1640+i*0.1565
    19 19 20 1.5042+i*1.3554
    20 20 21 0.4095+i*0.4784
    21 21 22 0.7089+i*0.9373
    22 3 23 0.4512+i*0.3083
    23 23 24 0.8980+i*0.7091
    24 24 25 0.8960+i*0.7011
    25 6 26 0.2030+i*0.1034
    26 26 27 0.2842+i*0.1447
    27 27 28 1.0590+i*0.9337
    28 28 29 0.8042+i*0.7006
    29 29 30 0.5075+i*0.2585
    30 30 31 0.9744+i*0.9630
    31 31 32 0.3105+i*0.3619
    32 32 33 0.3410+i*0.5362
    33 8  21 2.0+i*2.0
    34 9  15 2.0+i*2.0
    35 12 22 2.0+i*2.0
    36 18 33 0.5+i*0.5
    37 25 29 0.5+i*0.5   ];%支路参数矩阵，支路号、首末节点、支路参数
NodeM=[1 0
    2 0.1000+i*0.0600
    3 0.0900+i*0.0400
    4 0.1200+i*0.0800
    5 0.0600+i*0.0300
    6 0.0600+i*0.0200
    7 0.2000+i*0.1000
    8 0.2000+i*0.1000
    9 0.0600+i*0.0200
    10 0.0600+i*0.0200
    11 0.0450+i*0.0300
    12 0.0600+i*0.0350
    13 0.0600+i*0.0350
    14 0.1200+i*0.0800
    15 0.0600+i*0.0100
    16 0.0600+i*0.0200
    17 0.0600+i*0.0200
    18 0.0900+i*0.0400
    19 0.0900+i*0.0400
    20 0.0900+i*0.0400
    21 0.0900+i*0.0400
    22 0.0900+i*0.0400
    23 0.0900+i*0.0500
    24 0.4200+i*0.2000
    25 0.4200+i*0.2000
    26 0.0600+i*0.0250
    27 0.0600+i*0.0250
    28 0.0600+i*0.0200
    29 0.1200+i*0.0700
    30 0.2000+i*0.6000
    31 0.1500+i*0.0700
    32 0.2100+i*0.1000
    33 0.0600+i*0.0400 ]; %节点参数矩阵（电源节点负荷为0），节点号，功率

layerM=[1];      %节点分层矩阵，电源节点号记“1”
NU=zeros（1n）;   %上层节点矩阵（有33列的行矩阵）
BranchM（:4）=BranchM（:4）/Zb; %阻抗标幺化
NodeM（:2）=NodeM（:2）/Sb; %功率标幺化
A=zeros（n）;   %节点-节点关联矩阵A
for i1=1:b
    A（BranchM（i12）BranchM（i13））=1;
    A（BranchM（i13）BranchM（i12））=1;
end

%以下用循环求取矩阵layerM和NU
h=1;
while（min（NU（2:33）~=0）==0） %NU矩阵的2-最后都有上层节点了，表示循环结束了
    m=max（find（layerM（:h）））;  %m为矩阵layerM第h列非零元素的个数
    k=1;
    for i1=1:m
        g=layerM（i1h）;           %layerM的第i1行第h列元素
        ss=f