% 配电网前推回代潮流计算程序
% 使用IEEE 33节点配电系统作为算例，可实现弱环网情况下的潮流计算
clear
clc
countnum=0;
BranchData = [
1 2 0.0922 0.0470;
2 3 0.4930 0.2511;
3 4 0.3660 0.1864;
4 5 0.3811 0.1941;
5 6 0.8190 0.7070;
6 7 0.1872 0.6188;
7 8 0.7114 0.2351;
8 9 1.0300 0.7400;
9 10 1.0440 0.7400;
10 11 0.1966 0.0650;
11 12 0.3744 0.1238;
12 13 1.4680 1.1550;
13 14 0.5416 0.7129;
14 15 0.5910 0.5260;
15 16 0.7463 0.5450;
16 17 1.2890 1.7210;
17 18 0.7320 0.5740;
2 19 0.1640 0.1565;
19 20 1.5042 1.3554;
20 21 0.4095 0.4784;
21 22 0.7089 0.9373;
3 23 0.4512 0.3083;
23 24 0.8980 0.7091;
24 25 0.8960 0.7011;
6 26 0.2030 0.1034;
26 27 0.2842 0.1447;
27 28 1.0590 0.9337;
28 29 0.8042 0.7006;
29 30 0.5075 0.2585;
30 31 0.9744 0.9630;
31 32 0.3105 0.3619;
32 33 0.3410 0.5302;
]; % 支路，阻抗
NodeData = [
2 100.00 60.00;
3 90.00 40.00;
4 120.00 80.00;
5 60.00 30.00;
6 60.00 20.00;
7 200.00 100.00;
8 200.00 100.00;
9 60.00 20.00;
10 60.00 20.00;
11 45.00 30.00;
12 60.00 35.00;
13 60.00 35.00;
14 120.00 80.00;
15 60.00 10.00;
16 60.00 20.00;
17 60.00 20.00;
18 90.00 40.00;
19 90.00 40.00;
20 90.00 40.00;
21 90.00 40.00;
22 90.00 40.00;
23 90.00 50.00;
24 420.00 200.00;
25 420.00 200.00;
26 60.00 25.00;
27 60.00 25.00;
28 60.00 20.00;
29 120.00 70.00;
30 200.00 600.00;
31 150.00 70.00;
32 210.00 100.00;
33 60.00 40.00;
]; % 节点，负荷
UB = 12.66; % 电压基准 kV
SB = 10; % 功率基准 MVA
ZB = UB^2/SB; % 阻抗基准 ohm
BranchData(:,[3,4]) = BranchData(:,[3,4]) / ZB; % 阻抗标幺化
NodeData(:,[2,3]) = NodeData(:,[2,3]) / SB / 1000;% 功率标幺化
NN = 33; % 节点数
A0 = zeros(NN);
for n = 1:NN-1
A0(BranchData(n,1),BranchData(n,2)) = 1;
end % 形成 A0
% AssociatedMatrix=0;
%
% for n=2:NN-1
% AssociatedMatrix(n,n)=1;
% temp=BranchData(n-1,1);
% AssociatedMatrix(n,1:n-1)=AssociatedMatrix(temp,1:n-1);
% end
A0T = A0'; % 形成 A0 的转置
S = [0;-NodeData(:,2) - i*NodeData(:,3)]; % 形成 S
ZL = [0;BranchData(:,3) + i*BranchData(:,4)]; % 形成 ZL
V = ones(NN,1);
V(1) = 1; % 各个节点电压赋初值
IL(NN,1) = -conj(S(NN) / V(NN)); % 最末支路电流赋初值
Delta = 1; % 收敛判据赋初值
TempV = V; % 赋初值，用于记忆上次迭代结果
while Delta > 1e-8
countnum=countnum+1;
IN = conj(S ./ V); % 节点注入电流
for n = 1:NN-1
IL(NN-n) = A0(NN-n,NN-n+1:end) * IL(NN-n+1:end) - IN(NN-n);
end % 电流回代过程
for n = 2:NN
V(n) = A0T(n,1:n-1) * V(1:n-1) - ZL(n) * IL(n);
end % 电压前推过程
Delta = max(abs(V-TempV)); % 更新收敛判据
TempV = V; % 记忆迭代结果
end
Vangle(:,1)=abs(V);
Vangle(:,2)=angle(V)/3.1415*180;
for i=1:NN-1
st=BranchData(i,1);
en=BranchData(i,2);
Sij(i,1)=V(st)*conj((V(st)-V(en))/ZL(i+1));
Sji(i,1)=V(en)*conj((V(en)-V(st))/ZL(i+1));
end
