PQ分解法计算大电网潮流程序文件Word下载.docx
《PQ分解法计算大电网潮流程序文件Word下载.docx》由会员分享，可在线阅读，更多相关《PQ分解法计算大电网潮流程序文件Word下载.docx（13页珍藏版）》请在冰点文库上搜索。
%生成节点导纳矩阵
Y1=admittance（bus,line,2）;
%生成化简条件3的矩阵B1
Y2=admittance（bus,line,3）;
%生成化简条件3的矩阵B2
%-----------------------------------------------------
%%临时添加的测试数据
%nPQ=4;
nPV=0;
nSW=1;
nb=5;
%Y=[10.834-32.5i-1.667+5i-1.667+5i-2.5+7.5i-5+15i
%-1.667+5i12.917-38.75i-10+30i0-1.25+3.75i
%-1.667+5i-10+30i12.917-38.75i-1.25+3.75i0
%-2.5+7.5i0-1.25+3.75i3.75-11.25i0
%-5+15i-1.25+3.75i006.25-18.75i];
%
%Y1=[10.834-32.5i-1.667+5i-1.667+5i-2.5+7.5i-5+15i
%bus=[1100.20.21
%210-0.45-0.151
%310-0.4-0.051
%410-0.6-0.11
%51.060003];
%line=[521.25-3.75000
%2310-30000
%341.25-3.75000
%412.5-7.5000
%121.667-5000
%131.667-5000
%155-15000];
%-------------------------------------------------------
bus_PV0=bus（（nPQ+1）:
end,2）'
;
%1.05*ones（1,nPV+nSW）;
bus_U=[ones（1,nPQ）bus_PV0]'
%电压幅值
bus_e=zeros（nb,1）;
%电压角度
delta_P=zeros（nPQ+nPV,1）;
delta_Q=zeros（nPQ,1）;
%delta_e=zeros（nb-1,1）;
%delta_U=zeros（nPQ,1）;
c=0;
KP=1;
KQ=1;
%KPKQ用来判断有功、无功是否收敛
G=real（Y）;
B=imag（Y）;
B10=imag（Y1）;
B20=imag（Y2）;
%矩阵B0是进行化简三后的节点导纳矩阵虚部
%形成解耦潮流的系数矩阵B1和B2
B1=B10（1:
nb-1,1:
nb-1）;
B2=B20（1:
nPQ,1:
nPQ）;
whilec<
80
%求解PQ的不平衡量
forii=1:
nPQ+nPV
delta_P（ii）=bus（ii,4）;
forjj=1:
nb
delta_P（ii）=delta_P（ii）-bus_U（ii）*bus_U（jj）*（G（ii,jj）*cos（bus_e（ii）-bus_e（jj））+B（ii,jj）*sin（bus_e（ii）-bus_e（jj）））;
end
end
UP=diag（bus_U（1:
（nb-1）））;
%U矩阵利用各节点电压形成对角阵，来计算修正方程，对角线上的元素与bus_U列元素一一对应
error_P=UP\delta_P;
ifmax（abs（error_P））>
0.00001
delta_e=-（UP*B1）\error_P;
bus_e=bus_e+[delta_e;
0];
c=c+1;
KQ=1;
elseKP=0;
ifKQ~=0
else
break
nPQ
delta_Q（ii）=bus（ii,5）;
delta_Q（ii）=delta_Q（ii）-bus_U（ii）*bus_U（jj）*（G（ii,jj）*sin（bus_e（ii）-bus_e（jj））-B（ii,jj）*cos（bus_e（ii）-bus_e（jj）））;
UQ=diag（bus_U（1:
（nb-nPV-nSW）））;
error_Q=UQ\delta_Q;
ifmax（abs（error_Q））>
delta_U=-B2\error_Q;
bus_U=bus_U+[delta_U;
zeros（（nPV+nSW）,1）];
KP=1;
elseKQ=0;
ifKP~=0
%至此得到收敛的节点电压值
%----------------------------------------------
%--------------------------------------------
%对计算结果进行数据处理
%将节点结果用原节点编号表示
bus_Ue=zeros（nb,3）;
bus_Ue（:
[1,2,3]）=[nodenum（:
2）bus_Ubus_e/pi*180];
forjj=ii+1:
ifbus_Ue（ii,1）>
bus_Ue（jj,1）
t=bus_Ue（ii,:
bus_Ue（ii,:
）=bus_Ue（jj,:
bus_Ue（jj,:
）=t;
%r_U是收敛的电压表达成复数的形式
r_U=zeros（nb,1）;
fork=1:
r_U（k）=bus_U（k）*（cos（bus_e（k））+1i*sin（bus_e（k）））;
%计算平衡节点功率
SW_S=0;
SW_S=SW_S+r_U（nb）*conj（Y（nb,:
））*conj（r_U）;
%计算各支路功率Sij
line_S=zeros（nb,nb）;
line_S0=zeros（nb,nb）;
line_S（ii,jj）=r_U（ii）*（conj（r_U（ii））*conj（Y（ii,ii））+（conj（r_U（ii））-conj（r_U（jj）））*conj（Y（ii,jj）））;
%------------------------------------------
%把线路结果还原成原节点编号对应的结果
line_S0（nodenum（ii,2）,nodenum（jj,2））=line_S（ii,jj）;
line_P=real（line_S0）;
line_Q=imag（line_S0）;
%计算各支路损耗
delta_S=zeros（nl,1）;
nl
a=linenum（k,1）;
b=linenum（k,2）;
delta_S（k）=line_S0（a,b）+line_S（b,a）;
%计算网络总损耗
S0=sum（delta_S）;
%-------------------------------------------------
%将计算结果输入指定文件
fid=fopen（'
C:
\Users\lr\Desktop\matlab练习\训练题\大电网潮流计算\ieee14_out.txt'
'
wt'
fprintf（fid,'
节点号\t节点电压幅值\t节点电压角度\n'
fprintf（fid,'
%d\t%f\t%f\n'
k,bus_Ue（k,1）,bus_Ue（k,2））;
支路首端\t支路末端\t支路有功\t支路无功\t支路损耗\n'
%d\t\t%d\t\t%f\t%f\t%f\n'
linenum（k,1）,linenum（k,2）,line_P（linenum（k,1）,linenum（k,2））,line_Q（linenum（k,1）,linenum（k,2））,delta_S（k））;
平衡节点功率=%f\n'
SW_S）;
网络总损耗=%f\n'
S0）;
fclose（fid）;
function[bus,line,nPQ,nPV,nSW,nodenum]=change1_busline（bus,line）
%此函数用来对原始输入节点、线路数据进行重新编号
nPQ=0;
nSW=0;
%PQ=[];
PQ=zeros（nb,6）;
PV=zeros（nb,6）;
SW=zeros（nb,6）;
%PQPV平衡节点的个数
switchbus（k,6）
case1
nPQ=nPQ+1;
PQ（nPQ,:
）=bus（k,:
case2
nPV=nPV+1;
PV（nPV,:
case3
nSW=nSW+1;
SW（nSW,:
otherwise
disp（'
节点数据类型出错！
'
%生成重新编号后的节点数据矩阵
bus=[PQ;
PV;
SW];
nodenum=[（1:
%第一列为新的节点编号，第二列为对应的旧节点编号
1）=（1:
%至此实现了节点数据的重新编号
%------------------------------------------------------
%对线路数据重新编号
forii=1:
[r1,~]=find（nodenum（:
2）==line（ii,1））;
line（ii,1）=nodenum（r1,1）;
[r2,~]=find（nodenum（:
2）==line（ii,2））;
line（ii,2）=nodenum（r2,1）;
functionY=admittance（bus,line,c）
%此函数用来形成节点导纳矩阵
%line数组1.始端节点编号2.末端节点编号3.电阻4电抗5电导6电纳B/27.变比
%c是用来控制形成节点导纳矩阵的方式的，c=1形成一般的节点导纳矩阵，可以用来确定B2，c=2,形成化简条件3的节点导纳矩阵，确定B1
Y=zeros（nb,nb）;
zt=zeros（nl,1）;
yt=zeros（nl,1）;
ym=zeros（nl,1）;
I=zeros（nl,1）;
J=zeros（nl,1）;
K=zeros（nl,1）;
switchc
fork=1:
zt（k）=line（k,3）+1i*line（k,4）;
yt（k）=1/zt（k）;
ym（k）=line（k,5）+1i*line（k,6）;
I（k）=line（k,1）;
J（k）=line（k,2）;
K（k）=line（k,7）;
yt（k）=1/（line（k,3）+1i*line（k,4））;
ym=zeros（nl,1）;
ifline（k,7）~=1
K（k）=1;
elseK（k）=line（k,7）;
yt（k）=1/（1i*line（k,4））;
%针对不同线路设置节点导纳值
if（K（k）==1）&
&
（J（k）~=0）%普通线路
Y（I（k）,I（k））=Y（I（k）,I（k））+yt（k）+ym（k）;
Y（J（k）,J（k））=Y（J（k）,J（k））+yt（k）+ym（k）;
Y（I（k）,J（k））=Y（I（k）,J（k））-yt（k）;
Y（J（k）,I（k））=Y（I（k）,J（k））;
（J（k）==0）%对地支路
Y（I（k）,I（k））=Y（I（k）,I（k））+ym（k）;
if（K（k）>
1）&
（J（k）~=0）%非标准变比在j侧的变压器支路,ztym为折算到i侧的值
Y（J（k）,J（k））=Y（J（k）,J（k））+yt（k）/K（k）^2;
Y（I（k）,J（k））=Y（I（k）,J（k））-yt（k）/K（k）;
if（K（k）<
（J（k）~=0）
Y（J（k）,J（k））=Y（J（k）,J（k））+yt（k）*K（k）^2;
Y（I（k）,J（k））=Y（I（k）,J（k））-yt（k）*K（k）;
升级会员 