1.初衷
在做课程设计的时候，老师要求手算潮流，而非采用计算机算法求解（emmm，老师的心思真难猜啊），也就是采用电力系统分析课本上的前推回代求解潮流电压分布，嫌麻烦，试着用matalab编个程序实现。
2.程序的功能
这个程序可以用来针对辐射性网络求解潮流电压分布，分为一个主程序和两个子程序。如果有环网，那么就要配合此程序手动解环后，将环网分解为两个辐射性网络再来求解。程序中采用的都是标幺值。
3.思路
在编写代码时，首先需要考虑：
如何识别电网的拓扑结构，也就是说，如何知道一个节点是否还有分支。我的想法是：利用Branch矩阵中的第一列和第二列，规定功率流向必须是从第一列流向第二列，以此来确定电网的辐射情况。
4.算例
原件参数图如下所示，此参数图为有铭值，程序计算时使用的标幺值，所以需要在计算时化为标幺值，SB = 100MW，UB = 110KV
各节点功率情况
| 节点编号 | 有功功率 | 功率因数 |
|---|---|---|
| 14 | 35 | 0.85 |
| 17（发出） | 150 *0.9 | 115*0.9 |
| 2 | 15 | 0.8 |
| 3 | 25 | 0.8 |
| 4 | 20 | 0.8 |
| 5 | 20 | 0.8 |
| 6 | 21 | 0.80 |
| 7 | 15 | 0.82 |
| 16 | 30 | 0.82 |
| 18 （发出） | 平衡节点 | 平衡节点 |
| 9 | 26 | 0.85 |
| 10 | 24 | 0.83 |
| 11 | 16 | 0.80 |
5.代码
分为1个主程序和3个子程序：
主程序 main.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 前推回代手算潮流
clear;
% 母线Bus参数
%如果节点下面没有变压器，那么Bus参数的意思就是：
%bus_i（编号） type（节点类型） Pd Qd Gs Bs 变比（节点下面没有变压器的支路为0）
%如果这个节点后面还跟着变压器，那么Bus参数的意思就是：
%bus_i（编号） type（节点类型） Pd Qd Y Zt 变比
%节点编号规则为：PQ节点-PV节点-平衡节点
% Pd,Qd指的是消耗的功率，负荷消耗功率是负值
Bus = [
1 1 0 0 1.3674*10^(-5)*110^2/100 -6.844^(-5)*110^2/100 0;
2 1 -0.15 -0.15*tan(acos(0.8)) (13.24-1i*19.042)*10^-5*110^2/100 (1.915+1i*39.703)*100/110^2 1.1;
3 1 -0.25 -0.25*tan(acos(0.8)) (4.528-1i*27.68)*10^-5*110^2/100 (1.123+1i*39.703)*100/110^2 1.1;
4 1 -0.20 -0.20*tan(acos(0.8)) (3.834-1i*23.8)*10^-5*110^2/100 (1.498+1i*31.763)*100/110^2 1.1;
5 1 -0.20 -0.20*tan(acos(0.8)) (3.834-1i*23.8)*10^-5*110^2/100 (1.498+1i*31.763)*100/110^2 1.1;
6 1 -0.21 -0.21*tan(acos(0.8)) (13.24-1i*19.042)*10^-5*110^2/100 (1.915+1i*39.703)*100/110^2 1.1;
7 1 -0.15 -0.15*tan(acos(0.82)) (2.71-1i*16.116)*10^-5*110^2/100 (2.556+1i*50.82)*100/110^2 1.1;
8 1 0 0 1.448*10^(-5)*110^2/100 -8.032^(-5)*110^2/100 0;
9 1 -0.26 -0.26*tan(acos(0.85)) (3.834-1i*23.8)*10^-5*110^2/100 (1.498+1i*31.763)*100/110^2 1.1;
10 1 -0.24 -0.24*tan(acos(0.83)) (3.834-1i*23.8)*10^-5*110^2/100 (1.498+1i*31.763)*100/110^2 1.1;
11 1 -0.16 -0.16*tan(acos(0.8)) (13.24-1i*19.042)*10^-5*110^2/100 (1.915+1i*39.703)*100/110^2 1.1;
12 1 0 0 0 0 0;
13 1 0 0 0 0 0;
14 1 -0.35 -0.35*tan(acos(0.85)) 0 0 0;
15 1 0 0 0 0 0;
16 1 -0.3 -0.3*tan(acos(0.82)) 0 0 0;
17 1 1.50*0.9 1.15*0.9 0 0 0;
18 3 0 0 0 0 0;
];
% 发电机参数：母线节点编号 节点类型 P Q U θ
% 此处P，Q为发电机发出功率（正值）
Gen = [
17 1 1.50*0.9 1.15*0.9 1.05 0;
18 3 0 0 1.05 0;
];
% 交流分支线Branch参数：I侧母线 J侧母线 r x 1/2接地导纳 变比（没有变压器的支路为0）
% 变压器阻抗归算到I侧
Branch = [
1 2 8.10*100/110^2 12.27*100/110^2 0 0;
1 3 9.1*100/110^2 27.16*100/110^2 0 0;
10 1 6.50*100/110^2 19.4*100/110^2 0 0;
12 1 5.10*100/110^2 11.85*100/110^2 0 0;
1 13 0.175*100/110^2 -0.387*100/110^2 0 0;
3 4 5.78*100/110^2 13.43*100/110^2 0 0;
3 5 6.51*100/110^2 12.493*100/110^2 0 0;
8 6 5.5*100/110^2 19.15*100/110^2 0 0;
6 11 7.56*100/110^2 14.508*100/110^2 0 0;
11 7 13.5*100/110^2 20.45*100/110^2 0 0;
8 10 4.03*100/110^2 12.028*100/110^2 0 0;
15 8 0.11*100/110^2 -0.61*100/110^2 0 0;
12 9 3.40*100/110^2 7.90*100/110^2 0 0;
11 12 8.61*100/110^2 16.523*100/110^2 0 0;%此处11和12没有反映出他们的向下辐射情况。需要手算
13 14 0.175*100/110^2 5.423*100/110^2 0 1;
13 17 0.175*100/110^2 8.521*100/110^2 0 0.9091;
15 16 0.22*100/110^2 5.491*100/110^2 0 1;
15 18 0.11*100/110^2 8.541*100/110^2 0 0.9091;
];
%Bus参数后面跟着的变压器和Branch参数中的变压器的区别:
%为什么Bus参数会跟着变压器，而不都放到Branch里面呢？因为如果放到Branch里面，需要变压器支路的两端节点，
%那么就要再多标一些节点编号才行，就导致节点过于多，而这些节点编号就是为了求变压器损耗，所以就放到Bus里面，省下了节点编号的工作。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 求节点数n和PQ节点数m
%PQ节点包含联络节点，其P和Q都为0
n=size(Bus,1); %返回数组Bus的行数,即总节点数
%求PQ节点数m
m=n-1; % 先把平衡节点减掉
%然后把所有的PV节点都减掉，就求出PQ的节点数
for i=1:size(Bus,1)%返回Bus数组行数
if Bus(i,2)==2
m=m-1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 节点导纳矩阵的计算
[Y,B,G,Trans_pi]= MakeYbus( n,Branch,Bus );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 把Bus数据中的变压器参数变为Π型等效参数，更新Bus数据
for i=1:size(Bus,1)
if Bus(i,7) ~=0
%Bus(i,5)是Y，Bus(i,6)是Zt,Bus(i,7)是变比k
Bus(i,:)=[Bus(i,1:4) Bus(i,6)*Bus(i,7) (Bus(i,7)-1)/(Bus(i,7).*Bus(i,6)) + Bus(i,5) (1-Bus(i,7))/(Bus(i,7).^2*Bus(i,6))];
%如果有变压器，Trans更新以后参数变为%bus_i type Pd Qd Z_pi Y_pi1 Y_pi2
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 计算考虑变压器损耗的节点真实负荷(最大负荷，发电机出力为0.9)
Sj =[Bus(:,1) -(Bus(:,3)+1j*Bus(:,4))];%此处将所有的负荷功率变为正数
Si_true = zeros(18,2);%Si_true表示的意义是所有流出该节点的功率之和
for i = 1:size(Si_true,1)
%如果节点后面跟着变压器，考虑变压器的损耗
if Bus(i,7) ~=0
%Bus(:,5)是 Z_pi，Bus(:,6)是Y_pi1，Bus(:,7)是 Y_pi2
Sj_pie= Sj(i,2) + 1^2*conj(Bus(i,7));
Si_true(i,:)=[i Sj_pie+(real(Sj_pie)^2+imag(Sj_pie)^2)/1^2*Bus(i,5)+1^2*conj(Bus(i,6))];
else
%如果节点后面没有变压器，那么Si_true就是原始的负荷功率
Si_true(i,:) = Sj(i,:);
end
end
%把节点的自导纳也算进去，也当是功率经自导纳流出该节点。此题中就8和15节点有自导纳，原因是他们所接的三绕组变压器
%的等效电路中，有一个导纳，这个导纳不像线路的导纳参数一样，两边都有，所以作为节点的自导纳处理
for i = 1:size(Si_true,1)
if (Bus(i, 7) == 0)&&(Bus(i, 5) ~=0 || Bus(i, 6) ~=0)
Si_true(i,2) = Si_true(i,2) + 1^2*conj(Bus(i,5) +1i*Bus(i,6));
end
end
S1 = Si_true(1,2);%后面1需要重置时使用
S11 = Si_true(11,2);%后面11需要重置时使用
S12 = Si_true(12,2);%后面12需要重置时使用
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 前推
%%求辐射性网络
%求3变电站的向下辐射的功率分布
i=3;
[s3,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%求13变电站的向下辐射的功率分布
i = 13;
[s13,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%求1变电站的向下辐射的功率分布
i = 1;
[s1,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%此时求出的s1即为s右
%求11变电站的向下辐射的功率分布,此值不能使用子程序，因为Branch中不能反映真实的辐射情况
a = real(Si_true(Branch(10,2),2));
b = imag(Si_true(Branch(10,2),2));
c =Si_true(Branch(10,2),2) +(a^2 + b^2)*(Branch(10,3) + 1i*Branch(10,4));
s11(1,:) = [Branch(10,1) Branch(10,2) c];
Si_true(11,2) = Si_true(11,2) + c;
%求12变电站的向下辐射的功率分布
a = real(Si_true(Branch(13,2),2));
b = imag(Si_true(Branch(13,2),2));
c =Si_true(Branch(13,2),2) +(a^2 + b^2)*(Branch(13,3) + 1i*Branch(13,4));
s12(1,:) = [Branch(13,1) Branch(13,2) c];
Si_true(12,2) = Si_true(12,2) + c;
%%求环网
%先解环
z = -[1/Y(8,10) 1/Y(10,1) 1/Y(1,12) 1/Y(12,11) 1/Y(11,6) 1/Y(6,8);];
z_sum = sum(z(1,:));
s8_10 =(Si_true(6,2)*z(1,6) + Si_true(11,2)*sum(z(1,5:6))+Si_true(12,2)*sum(z(1,4:6))+Si_true(1,2)*sum(z(1,3:6))+Si_true(10,2)*sum(z(1,2:6)))/z_sum;
s8_6 =(Si_true(6,2)*sum(z(1,1:5)) + Si_true(11,2)*sum(z(1,1:4))+Si_true(12,2)*sum(z(1,1:3))+Si_true(1,2)*sum(z(1,1:2))+Si_true(10,2)*z(1,1))/z_sum;
s10_1 = s8_10 - Si_true(10,2);
s1_12 = s10_1 - Si_true(1,2);
s6_11 = s8_6 - Si_true(6,2);
s11_12 = s6_11 - Si_true(11,2);
disp('s8-10');disp(s8_10);
disp('s8_6');disp(s8_6);
disp('s10-1');disp(s10_1);
disp('s1-12');disp(s1_12);
disp('s6-11');disp(s6_11);
disp('s11-12');disp(s11_12);
%确定功率分点为：12,人工修正Branch矩阵
%使用时要在此处设置断点，从而确定功率分点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 人工修正Branch矩阵
Branch = [
1 2 8.10*100/110^2 12.27*100/110^2 0 0;
1 3 9.1*100/110^2 27.16*100/110^2 0 0;
10 1 6.50*100/110^2 19.4*100/110^2 0 0;
1 12 5.10*100/110^2 11.85*100/110^2 0 0;
1 13 0.175*100/110^2 -0.387*100/110^2 0 0;
3 4 5.78*100/110^2 13.43*100/110^2 0 0;
3 5 6.51*100/110^2 12.493*100/110^2 0 0;
8 6 5.5*100/110^2 19.15*100/110^2 0 0;
6 11 7.56*100/110^2 14.508*100/110^2 0 0;
11 7 13.5*100/110^2 20.45*100/110^2 0 0;
8 10 4.03*100/110^2 12.028*100/110^2 0 0;
15 8 0.11*100/110^2 -0.61*100/110^2 0 0;
12 9 3.40*100/110^2 7.90*100/110^2 0 0;
11 12 8.61*100/110^2 16.523*100/110^2 0 0;
13 14 0.175*100/110^2 5.423*100/110^2 0 1;
13 17 0.175*100/110^2 8.521*100/110^2 0 0.9091;
15 16 0.22*100/110^2 5.491*100/110^2 0 1;
15 18 0.11*100/110^2 8.541*100/110^2 0 0.9091;
];
%成功在12节点处解环，环网转化为两个辐射性网络，再用子程序求解
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 求环网解环后的第一个辐射性网络
x = 12;
Si_true(x,2) = s1_12;
% 求1变电站的向下辐射的功率分布
% i = 1;
% Si_true(1,2) = S1;%因为在计算右边辐射性网络时，此节点已经算了一次，所以要
%重置为该节点最初的负荷情况，否则，1-3，1-2，1-13支路就算了两次。
% [s1,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%
%求1变电站的向下辐射的功率分布
i =1;
Si_true(1,2) = S1;%因为在计算右边辐射性网络时，此节点已经算了一次，所以要
% 重置为该节点最初的负荷情况，否则，12-9支路就算了两次。
[s1,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%求10变电站的向下辐射的功率分布
i = 10;
[s10,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 求求环网解环后的第二个辐射性网络
x = 12;
Si_true(x,2) = s11_12;
%求11变电站的向下辐射的功率分布
i = 11;
S11 = Si_true(11,2);%重置
[s11,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%求6变电站的向下辐射的功率分布
i = 6;
[s6,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%求8变电站的向下辐射的功率分布
i = 8;
[s8,Si_true]=qiantui( Si_true,Branch,i,Trans_pi);
%求15变电站的向下辐射的功率分布,Branch参数体现不出辐射情况,不能调用子程序
a = real(Si_true(Branch(12,2),2));
b = imag(Si_true(Branch(12,2),2));
c =Si_true(Branch(12,2),2) +(a^2 + b^2)*(Branch(12,3) + 1i*Branch(12,4));
s15(1,:) = [Branch(12,1) Branch(12,2) c];
Si_true(15,2) = Si_true(15,2) + c;
d = Si_true(Branch(17,2),2) + 1^2*conj(Trans_pi(3,5));
c = d + (real(d)^2 + imag(d)^2)/1^2*(Trans_pi(3,3)) + 1^2*conj(Trans_pi(3,4));
s15(2,:) = [Branch(17,1) Branch(17,2) c];
Si_true(15,2) = Si_true(15,2) + c;
%求18变电站的向下辐射的功率分布,Branch参数体现不出辐射情况,不能调用子程序
d = Si_true(Branch(18,1),2) + 1^2*conj(Trans_pi(4,4));
c = d + (real(d)^2 + imag(d)^2)/1^2*(Trans_pi(4,3)) + 1^2*conj(Trans_pi(4,5));
s18(1,:) = [Branch(18,2) Branch(18,1) c];
Si_true(18,2) = Si_true(18,2) + c;
s_huizong =[s1;s3;s6;s8;s10;s11;s12;s13;s15;s18];
disp(vpa(Si_true,3));
disp(vpa(s_huizong,3));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 求解各节点电压
%初始化U，8节点为1.05+1i*0，其他暂时未知
U = zeros(18,2);
U(8,:) =[8 1.05+1i*0];
%求i是8节点的时候，j的电压,确定6，10的电压。
i = 8;
[u8,U]= huidai( s_huizong,Y,U,i);
%求i是6节点的时候，j的电压，确定11的电压。
i = 6;
[u6,U]= huidai( s_huizong,Y,U,i);
%求i是11节点的时候，j的电压，确定7，12的电压。
i = 11;
[u11,U]= huidai( s_huizong,Y,U,i);
%求i是12节点的时候，j的电压，确定1，9的电压。
i = 12;
[u12,U]= huidai( s_huizong,Y,U,i);
%求i是10节点的时候，j的电压，确定1的电压。
i = 10;
[u10,U]= huidai( s_huizong,Y,U,i);
%求i是1节点的时候，j的电压，确定2，3，13的电压。
i = 1;
[u1,U]= huidai( s_huizong,Y,U,i);
%求i是13节点的时候，j的电压，确定14，17的电压。
i = 13;
[u13,U]= huidai( s_huizong,Y,U,i);
%求i是3节点的时候，j的电压，确定4，5的电压。
i = 3;
[u3,U]= huidai( s_huizong,Y,U,i);
disp(U);
子程序 MakeYbus.m
function [Y,B,G,Trans_pi ]= MakeYbus( n,Branch,Bus )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 变压器π型等效阻抗参数
%Branch(:,4)*1i是阻抗，Branch(:,6)是变比
% 求出变压器个数count
count = 0;
for i=1:size(Branch,1)
if Branch(i,6) ~=0
count = count + 1;
end
end
% 求变压器π型等效阻抗参数
Zt=zeros(count,3);
Trans_pi=zeros(count,5);
count = 0;
for i=1:size(Branch,1)
if Branch(i,6) ~=0
count= count+1;
Zt(count,1)=(Branch(i,3)+Branch(i,4)*1i)*Branch(i,6);
Zt(count,2)=(Branch(i,6)-1)/(Branch(i,6).*(Branch(i,3)+Branch(i,4)*1i));
Zt(count,3)=(1-Branch(i,6))/(Branch(i,6).^2*(Branch(i,3)+Branch(i,4)*1i));
%变压器π型参数：I侧母线 J侧母线 Z_pi Y_pi1 Y_pi2
Trans_pi(count,:)=[Branch(i,1:2) Zt(count,1) Zt(count,2) Zt(count,3)];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 节点导纳矩阵的计算
Y=zeros(n); %新建节点导纳矩阵
y=zeros(n); %网络中的真实导纳标幺值yij
%计算y(i,j)
for i=1:size(Branch,1) %与交流线联结的真实导纳
ii=Branch(i,1); jj=Branch(i,2);
y(ii,jj)=vpa(1/(Branch(i,3)+Branch(i,4)*1i),7);%保留七位有效数字
y(jj,ii)=vpa(y(ii,jj),7);
end
% Trans_pi=[Trans(:,1:2) Zt(:,1) 1./Zt(:,2) 1./Zt(:,3)];
for i=1:size(Trans_pi,1) %与变压器联结的真实导纳
ii=Trans_pi(i,1); jj=Trans_pi(i,2);
y(ii,jj)=1/Trans_pi(i,3);
y(jj,ii)=y(ii,jj);
end
%计算y(i,i)
for i=1:size(Branch,1) %与交流线联结的对地导纳
ii=Branch(i,1); jj=Branch(i,2);
y(ii,ii)=y(ii,ii) + Branch(i,5)*1i;
y(jj,jj)=y(jj,jj) + Branch(i,5)*1i;
end
% Trans_pi=[Trans(:,1:2) Zt(:,1) Zt(:,2) Zt(:,3)];
for i=1:size(Trans_pi,1) %与变压器联结的对地导纳
ii=Trans_pi(i,1); jj=Trans_pi(i,2);
y(ii,ii)=y(ii,ii)+Trans_pi(i,4); %朝I侧归算。输入数据时要把归算测当i
y(jj,jj)=y(jj,jj)+Trans_pi(i,5);
end
%算上补偿电容
for i=1:size(Bus,1)
if (Bus(i, 7) == 0)&&(Bus(i, 5) ~=0 || Bus(i, 6) ~=0)
ii=Bus(i,1);
y(ii,ii)=y(ii,ii)+Bus(i,5)+Bus(i, 6)*1i;
end
end
%由y计算Y
ysum=sum(y,1); %每一行求和
for i=1:n
for j=1:n
if i==j
Y(i,j)=ysum(i);
else
Y(i,j)=-y(i,j);
end
end
end
G=real(Y); %电导矩阵
B=imag(Y); %电纳矩阵
子程序 qiantui.m
%此子函数用于求某节点向下辐射（流出该节点）的功率分布，存放于s矩阵。
%同时可以将流出该节点的功率相加，存放于Si_true中
function [s,Si_true]=qiantui( Si_true,Branch,i,Trans_pi)
%i是节点，Si_true表示的是流出节点的功率和
count = 0;
%求出i节点向下辐射（功率流出该节点）的支路数
for x = 1:size(Branch(:,1))
if Branch(x,1)== i
count = count + 1;
end
end
%建立s矩阵，存放某节点向下辐射（流出该节点）的功率分布
s = zeros(count,3);
count = 0;
n = 0;
%识别向下辐射（功率流出节点）的信息，是靠Branch矩阵的第一列和第二列实现，
%功率由第一列向第二列辐射，在未知功率分布情况下，先不管，因为也没求这些未知
%流向的功率分布。等到知道了功率流向以后，需要手动更新Branch矩阵的第一、二列。
for x = 1:size(Branch(:,1))
if Branch(x,1)== i%识别功率流向
count = count + 1;
%分为两种情况，支路有没有变压器
%没有变压器
if Branch(x,6) == 0
a = real(Si_true(Branch(x,2),2));
b = imag(Si_true(Branch(x,2),2));
c = Si_true(Branch(x,2),2) +(a^2 + b^2)*(Branch(x,3) + 1i*Branch(x,4));
else
%有变压器
n = n + 1;
d = Si_true(Branch(x,2),2) + 1^2*conj(Trans_pi(n,5));
c = d + (real(d)^2 + imag(d)^2)/1^2*(Trans_pi(n,3)) + 1^2*conj(Trans_pi(n,4));
end
%s矩阵存放该节点向下辐射（流出该节点）的功率分布，格式为
%节点i 节点j 从i流向j的功率
s(count,:) = [Branch(x,1) Branch(x,2) c];
%更新Si_true中流出该节点的功率和
Si_true(i,2) = Si_true(i,2) + s(count,3);
end
end
子程序 huidai.m
% 回代求解各节点电压。
function [u,U]=huidai( s_huizong,Y,U,i)
%i是节点，s_huizong 表示的是从i节点到j节点的功率分布,U矩阵存放j的的电压分布
count = 0;
%求出i节点向下辐射的支路数
for x = 1:size(s_huizong(:,1))
if s_huizong(x,1)== i
count = count + 1;
end
end
u =zeros(count,2);
count = 0;
for x = 1:size(s_huizong(:,1))
if s_huizong(x,1)== i%识别功率流向
j = s_huizong(x,2);
count =count + 1;
u(count,:) =[j U(i,2) - conj(s_huizong(x,3))/abs( U(i,2))*(-1/Y(i,j))];
U(j,:) =[j u(count,2)];
end
end
6.结果
先设置一个断点，打印出环网的粗略估计的功率分布。
s8-10
0.4993 + 0.9347i
s8_6
0.5231 + 0.7508i
s10-1
0.2522 + 0.6918i
s1-12
0.3141 + 0.0735i
s6-11
0.2951 + 0.5278i
s11-12
-0.0432 + 0.1805i
确定功率分点为12，修正Branch矩阵，继续执行程序
结果如下
流出某节点的功率和：
[ 1.0, 0.257 + 0.702i]
[ 2.0, 0.168 + 0.174i]
[ 3.0, 0.68 + 0.745i]
[ 4.0, 0.207 + 0.228i]
[ 5.0, 0.207 + 0.228i]
[ 6.0, 0.729 + 1.0i]
[ 7.0, 0.155 + 0.157i]
[ 8.0, 1.38 + 2.42i]
[ 9.0, 0.267 + 0.245i]
[ 10.0, 0.534 + 1.03i]
[ 11.0, 0.458 + 0.697i]
[ 12.0, - 0.0432 + 0.181i]
[ 13.0, - 0.994 - 0.54i]
[ 14.0, 0.35 + 0.217i]
[ 15.0, 1.69 + 2.6i]
[ 16.0, 0.3 + 0.209i]
[ 17.0, - 1.35 - 1.04i]
[ 18.0, 1.69 + 3.01i]
功率分布：
[ 1.0, 2.0, 0.172 + 0.18i]
[ 1.0, 3.0, 0.757 + 0.974i]
[ 1.0, 12.0, 0.319 + 0.0837i]
[ 1.0, 13.0, - 0.992 - 0.544i]
[ 3.0, 4.0, 0.211 + 0.238i]
[ 3.0, 5.0, 0.212 + 0.238i]
[ 6.0, 11.0, 0.501 + 0.781i]
[ 8.0, 6.0, 0.799 + 1.25i]
[ 8.0, 10.0, 0.579 + 1.17i]
[ 10.0, 1.0, 0.287 + 0.792i]
[ 11.0, 7.0, 0.16 + 0.165i]
[ 11.0, 12.0, - 0.0407 + 0.185i]
[ 12.0, 9.0, 0.271 + 0.254i]
[ 13.0, 14.0, 0.35 + 0.225i]
[ 13.0, 17.0, - 1.34 - 0.764i]
[ 15.0, 8.0, 1.39 + 2.38i]
[ 15.0, 16.0, 0.3 + 0.215i]
[ 18.0, 15.0, 1.69 + 3.01i]
电压分布：
1.0000 + 0.0000i 0.7665 - 0.0214i
2.0000 + 0.0000i 0.7276 - 0.0284i
3.0000 + 0.0000i 0.4072 - 0.1475i
4.0000 + 0.0000i 0.3228 - 0.1754i
5.0000 + 0.0000i 0.3242 - 0.1686i
6.0000 + 0.0000i 0.8273 - 0.0665i
7.0000 + 0.0000i 0.6096 - 0.0929i
8.0000 + 0.0000i 1.0500 + 0.0000i
9.0000 + 0.0000i 0.6065 - 0.0690i
10.0000 + 0.0000i 0.9210 - 0.0177i
11.0000 + 0.0000i 0.6768 - 0.0801i
12.0000 + 0.0000i 0.7383 - 0.0575i
13.0000 + 0.0000i 0.7661 - 0.0266i
14.0000 + 0.0000i 0.7523 - 0.0467i
0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i
17.0000 + 0.0000i 0.8322 + 0.0844i
0.0000 + 0.0000i 0.0000 + 0.0000i
1090
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
