% m发电机节点,n负荷节点,ij为节点
% IEEE39节点算例中:m=10,n=21,46 条支路,N=46
% chaoliujieshu=min(p_m,p_n)*p_mn/p_ij_of_mn;% p_m为m节点发出的功率,p_n为n节点需要的功率,p_ij_of_mn为mn节点队中ij占的份额
% p_ij_from_m=A_ij_nixu'*p_ij*p_m/p_i;%p_ij,p_i可从matpower中获得,p_m输入m文件里
% p_ij_of_mn=p_ij_from_m*p_ij_to_n/p_ij;
% p_mn=p_n*A_ij_nixu*p_m/p_ij;
% result.branch(:,1:2)%用来计算zhilu_ij矩阵;
% result.branch(:,14/16)%用来计算p_ij,p_i,p_j;
%%拟返回函数
% function x=cljs(a)
% [paixu,pos]=sort(zhiluchaoliujieshu);
% x=pos(N);
%%统计支路
result=runpf('case39');
s1=size(result.branch);
s2=size(result.bus);
N=s1(1,1);
M=s2(1,1);
zhilu(M,M)=0;
for i=1:N
zhilu(result.branch(i,1),result.branch(i,2))=1;
end
for i=1:M
for j=1:M
if zhilu(j,i)~=0
zhilu(i,j)=zhilu(j,i);
else
end
end
end
%%算支路潮流和节点潮流
zhiluchaoliu(M,M)=0;
for i=1:N
zhiluchaoliu(result.branch(i,1),result.branch(i,2))=(abs(result.branch(i,14))+abs(result.branch(i,16)))/2;
end
for i=1:M
for j=1:M
if zhiluchaoliu(j,i)~=0
zhiluchaoliu(i,j)=zhiluchaoliu(j,i);
else
end
end
end
jiedianchaoliu(M)=0;
jiedianchaoliuru(M)=0;
jiedianchaoliuchu(M)=0;
for i=1:M
for j=1:N
if result.branch(j,1)==i
if result.branch(j,14)>0
jiedianchaoliuchu(i)=jiedianchaoliuchu(i)+result.branch(j,14);
else
jiedianchaoliuru(i)=jiedianchaoliuru(i)+result.branch(j,14);
end
elseif result.branch(j,2)==i
if result.branch(j,16)>0
jiedianchaoliuchu(i)=jiedianchaoliuchu(i)+result.branch(j,16);
else
jiedianchaoliuru(i)=jiedianchaoliuru(i)+result.branch(j,16);
end
else
end
end
jiedianchaoliu(i)=max(abs(jiedianchaoliuchu(i)),abs(jiedianchaoliuru(i)));
end
%%顺序矩阵,出线集
shunxu(M,M)=0;
for i=1:M
for j=1:M
if i==j;
shunxu(i,j)=1;
else
shunxu(i,j)=0;
end
end
end
%%逆序矩阵,入线集
nixu(M,M)=0;
for i=1:M
for j=1:M
if i==j;
nixu(i,j)=1;
else
nixu(i,j)=0;
end
end
end
for i=1:N
if result.branch(i,14)>0
shunxu(result.branch(i,1),result.branch(i,2))=-1;
nixu(result.branch(i,2),result.branch(i,1))=-1;
else
shunxu(result.branch(i,2),result.branch(i,1))=-1;
nixu(result.branch(i,1),result.branch(i,2))=-1;
end
end
for i=1:M
for j=1:M
if shunxu(i,j)==-1
shunxu(i,j)=(-1)*zhiluchaoliu(i,j)/jiedianchaoliu(j);
else
end
end
end
for i=1:M
for j=1:M
if nixu(i,j)==-1
nixu(i,j)=(-1)*zhiluchaoliu(i,j)/jiedianchaoliu(j);
else
end
end
end
%%计算潮流介数
for i=1:size(result.gen)
G(i,1)=result.gen(i,1);
G(i,2)=result.gen(i,2);
end
for i=1:size(result.bus)
if result.bus(i,3)~=0
L(i,1)=result.bus(i,1);
L(i,2)=result.bus(i,3);
end
end
L(all(L==0,2),:)=[];
s3=size(G);
s4=size(L);
chaoliujieshu(M,M)=0;
nixuni=inv(nixu);
shunxuni=inv(shunxu);
for m=1:s3(1,1)
for n=1:s4(1,1)
for i=1:M
for j=1:M
if zhilu(i,j)==1&&nixuni(L(n,1),G(m,1))~=0
chaoliujieshu(i,j)=chaoliujieshu(i,j)+min(G(m,2),L(n,2))*zhiluchaoliu(i,j)*jiedianchaoliu(L(n,1))*nixuni(i,G(m,1))*shunxuni(i,L(n,1))/(jiedianchaoliu(i)*jiedianchaoliu(i)*nixuni(L(n,1),G(m,1)));
else
end
end
end
end
end
zhiluchaoliujieshu(N)=0;
for i=1:N
zhiluchaoliujieshu(i)=chaoliujieshu(result.branch(i,1),result.branch(i,2));
end
%%结果显示
daipaixu(N,2)=0;
for i=1:N
daipaixu(i,1)=i;
daipaixu(i,2)=zhiluchaoliujieshu(i);
end
for i=1:N
if result.branch(i,9)~=0
daipaixu(i,:)=[];
else
end
end
A=max(daipaixu);
B=min(daipaixu);
s5=size(daipaixu);
for i=1:s5(1,1)
guiyihua(i)=(daipaixu(i,2)-B)/(A-B);
end
plot(guiyihua)
评论9