clc
clear
N=41;
b=40;
UB = 10; % 电压基准 kV
SB = 10; % 功率基准 MVA
ZB = UB^2/SB; % 阻抗基准 omh
%% 数据初始化:branch、bus为matpower格式,A、S、U为本程序格式
branch = [
1 2 0.0207+0.0384i 0.0207+0.0384i 0.0207+0.0384i 0 0 0
2 3 0.0160+0.0220i 0.0160+0.0220i 0.0160+0.0220i 0 0 0
3 4 0.0160+0.0220i 0.0160+0.0220i 0.0160+0.0220i 0 0 0
4 5 0.0602+0.0191i 0.0602+0.0191i 0.0602+0.0191i 0 0 0
5 6 0.0602+0.0191i 0.0602+0.0191i 0.0602+0.0191i 0 0 0
6 7 0.0813+0.0256i 0.0813+0.0256i 0.0813+0.0256i 0 0 0
7 8 0.0813+0.0256i 0.0813+0.0256i 0.0813+0.0256i 0 0 0
8 9 0.0158+0.0255i 0.0158+0.0255i 0.0158+0.0255i 0 0 0
9 10 0.0203+0.0376i 0.0203+0.0376i 0.0203+0.0376i 0 0 0
10 11 0.0203+0.0376i 0.0203+0.0376i 0.0203+0.0376i 0 0 0
11 12 0.0247+0.0457i 0.0247+0.0457i 0.0247+0.0457i 0 0 0
12 13 0.0048+0.0078i 0.0048+0.0078i 0.0048+0.0078i 0 0 0
13 14 0.0063+0.0102i 0.0063+0.0102i 0.0063+0.0102i 0 0 0
14 15 0.0057+0.0092i 0.0057+0.0092i 0.0057+0.0092i 0 0 0
15 16 0.0046+0.0075i 0.0046+0.0075i 0.0046+0.0075i 0 0 0
16 17 0.0016+0.0015i 0.0016+0.0015i 0.0016+0.0015i 0 0 0
17 18 0.0022+0.0016i 0.0022+0.0016i 0.0022+0.0016i 0 0 0
18 19 0.0014+0.0009i 0.0014+0.0009i 0.0014+0.0009i 0 0 0
3 20 0.0116+0.0118i 0.0116+0.0118i 0.0116+0.0118i 0 0 0
20 21 0.0075+0.0103i 0.0075+0.0103i 0.0075+0.0103i 0 0 0
21 22 0.0066+0.0091i 0.0066+0.0091i 0.0066+0.0091i 0 0 0
22 23 0.0084+0.0115 0.0084+0.0115 0.0084+0.0115 0 0 0
23 24 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
24 25 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
25 26 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
26 27 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
20 28 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
28 29 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
28 30 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
21 31 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
31 32 0.0099+0.0135i 0.0099+0.0135i 0.0099+0.0135i 0 0 0
32 33 0.0013+0.0006i 0.0013+0.0006i 0.0013+0.0006i 0 0 0
22 34 0.0090+0.0123i 0.0090+0.0123i 0.0090+0.0123i 0 0 0
34 35 0.0090+0.0123i 0.0090+0.0123i 0.0090+0.0123i 0 0 0
35 36 0.0054+0.0023i 0.0054+0.0023i 0.0054+0.0023i 0 0 0
25 37 0.0097+0.0065i 0.0097+0.0065i 0.0097+0.0065i 0 0 0
37 38 0.0054+0.0023i 0.0054+0.0023i 0.0054+0.0023i 0 0 0
8 39 0.1313+0.0554i 0.1313+0.0554i 0.1313+0.0554i 0 0 0
13 40 0.0009+0.0011i 0.0009+0.0011i 0.0009+0.0011i 0 0 0
14 41 0.0006+0.0010i 0.0006+0.0010i 0.0006+0.0010i 0 0 0
];
bus = [
1 0 0 0 0 0 0 1 0 1 120 1 -120
2 0.032 0.019 0.033 0.02 0.035 0.021 1 0 1 120 1 -120
3 0.03 0.013 0.031 0.015 0.029 0.013 1 0 1 120 1 -120
4 0.045 0.03 0 0 0.035 0.024 1 0 1 120 1 -120
5 0.02 0.01 0.02 0.01 0.02 0.01 1 0 1 120 1 -120
6 0.02 0.006 0.02 0.007 0.02 0.007 1 0 1 120 1 -120
7 0.065 0.033 0.07 0.034 0.065 0.033 1 0 1 120 1 -120
8 0.07 0.034 0.065 0.033 0.065 0.033 1 0 1 120 1 -120
9 0.02 0.007 0.018 0.006 0.022 0.007 1 0 1 120 1 -120
10 0.021 0.007 0.02 0.007 0 0 1 0 1 120 1 -120
11 0.014 0.009 0.016 0.011 0.015 0.01 1 0 1 120 1 -120
12 0.02 0.011 0.02 0.012 0.02 0.012 1 0 1 120 1 -120
13 0.021 0.012 0.019 0.011 0.02 0.012 1 0 1 120 1 -120
14 0.04 0.028 0.038 0.027 0.042 0.025 1 0 1 120 1 -120
15 0 0 0.019 0.003 0.02 0.003 1 0 1 120 1 -120
16 0.019 0.006 0.02 0.007 0.021 0.007 1 0 1 120 1 -120
17 0.019 0.006 0.021 0.007 0.02 0.007 1 0 1 120 1 -120
18 0.03 0.014 0.03 0.013 0.03 0.013 1 0 1 120 1 -120
19 0.033 0.015 0.029 0.013 0.028 0.012 1 0 1 120 1 -120
20 0.029 0.013 0.028 0.012 0.033 0.015 1 0 1 120 1 -120
21 0.029 0.012 0.03 0.013 0.031 0.015 1 0 1 120 1 -120
22 0.028 0.012 0.033 0.015 0.029 0.013 1 0 1 120 1 -120
23 0.03 0.016 0.031 0.017 0.029 0.017 1 0 1 120 1 -120
24 0.13 0.06 0.14 0.07 0.15 0.07 1 0 1 120 1 -120
25 0.15 0.07 0.13 0.07 0.14 0.06 1 0 1 120 1 -120
26 0.02 0.008 0.02 0.008 0.02 0.009 1 0 1 120 1 -120
27 0.018 0.007 0.022 0.009 0.02 0.009 1 0 1 120 1 -120
28 0.019 0.006 0.022 0.008 0.019 0.006 1 0 1 120 1 -120
29 0.038 0.023 0.042 0.025 0.04 0.022 1 0 1 120 1 -120
30 0.06 0.18 0.07 0.21 0.07 0.21 1 0 1 120 1 -120
31 0.045 0.02 0.051 0.023 0.054 0.027 1 0 1 120 1 -120
32 0.07 0.033 0.072 0.035 0.068 0.032 1 0 1 120 1 -120
33 0.02 0.013 0.02 0.014 0.02 0.013 1 0 1 120 1 -120
34 0.03 0.015 0.18 0.06 0 0 1 0 1 120 1 -120
35 0.02 0.01 0.12 0.04 0.05 0.032 1 0 1 120 1 -120
36 0.15 0.08 0.09 0.03 0.06 0.026 1 0 1 120 1 -120
37 0.09 0.03 0.16 0.04 0.09 0.05 1 0 1 120 1 -120
38 0.1 0.05 0.1 0.03 0.08 0.038 1 0 1 120 1 -120
39 0.2 0.7 0.2 0.05 0.07 0.025 1 0 1 120 1 -120
40 0.05 0.02 0.08 0.04 0.06 0.025 1 0 1 120 1 -120
41 0.09 0.045 0.05 0.03 0 0 1 0 1 120 1 -120
];
A = [branch(:, 1:8);];
A(:, 3:end) = A(:, 3:end)./ZB; % 阻抗标幺化
num =1:size(A, 1);
A = [num', A]; % A 第一列支路号,第二列支路始端,第三列支路末端,第四到九列为相间电阻
S=[ bus(:, 1), -bus(:, 2)-i*bus(:, 3), -bus(:, 4)-i*bus(:, 5), -bus(:, 6)-i*bus(:, 7);]; %负号表示流出
S(:, 2:end) = S(:, 2:end)./SB; % 负荷标幺化
U = [bus(:, 1), bus(:, 8).*exp(-i*pi*bus(:, 9)/180), ...
bus(:, 10).*exp(-i*pi*bus(:, 11)/180), bus(:, 12).*exp(-i*pi*bus(:, 13)/180);]; % 电压赋初值
U = [U ones(size(U, 1), 1)];
%% 生成节点的关联矩阵
A0 = zeros(N);
for m = 1:b
A0(A(m,2),A(m,3)) = 1;
end % 形成 A0
%% 赋初值
DeltaU=1;
IterMax = 500;
k=0;
IN=zeros(size(S, 1),4); %IN=zeros(N,4);
IN(:,1)=S(:,1);
IL = zeros(b, 6);
IL(:,1:3) = A(:,1:3);
%% 开始迭代
while max(max(DeltaU))>1e-5
if (k > IterMax)
break;
end
k = k + 1;
IL(N, 4:6) = conj(S(N, 2:4) ./ U(N, 2:4));
PQ(N,4:6)=abs(IL(N,4:6)).^2.*real(A(N-1,4:6))
IN(:, 2:4) = conj(S(:, 2:4) ./ U(:,2:4)); % 计算节点的负荷注入电流,为负
for m=1:N-1
IL(N-m,4:6) = IN(N-m,2:4) + A0(N-m,N-m+1:end) * IL(N-m+1:end, 4:6) ; % 前推计算支路电流
end
TempU = U(:, 2:4); % 保存上一次电压结果
for m=2:N % 回代计算节点电压
ZL=[A(m-1,4),A(m-1,7),A(m-1,8);
A(m-1,7),A(m-1,5),A(m-1,9);
U(A(m-1, 3),2:4) = U(A(m-1, 2),2:4)+(ZL*(IL(m,4:6).')).'; %更新回代电压数据,不必排除PV节点
end
DeltaU=abs(abs(TempU) - abs(U(:, 2:4))); %收敛判断
end
%% 结果输出
Uabs = [U(:, 1) abs(U(:,2:4))]
Uang = [U(:, 1) angle(U(:,2:4))*180/pi]
disp(['Message:']);
if k <= IterMax
disp([' Successfully solved in ' num2str(k) ' iterations.']);
else
disp([' Maximum ' num2str(IterMax) ' iterations reached.']);
end