function?[YUPQdeltaSijaSSijSjisumdeltaS]=chaoliujisuan（）
%電力系統潮流計算程序；
%輸出：U——節點電壓，P--節點有功，Q--節點無功，deltaSij--支路功率損耗，
%Sij--從節點i流向節點j的功率S--節點復功率sumdeltaS--網絡總損耗
%輸入參數：point為節點信息矩陣，zhilu為支路信息矩陣；
[x]=xlsread（‘NRinput.xls‘‘A2:A2‘）;%從exel中讀取節點數x
[y]=xlsread（‘NRinput.xls‘‘B2:B2‘）;%從exel中讀取支路數y
e=xlsread（‘NRinput.xls‘‘B4:B4‘）;%誤差要求
[point]=xlsread（‘NRinput.xls‘‘D3:H50‘）;%從exel中讀取節點信息矩陣
[zhilu]=xlsread（‘NRinput.xls‘‘J3:Q50‘）;%從exel中讀取支路信息矩陣
TYPE=zeros（x1）;%TYPE為節點類型矩陣
U=zeros（x1）;%U為節點電壓矩陣
a=zeros（x1）;%a為節點電壓相角矩陣
P=zeros（x1）;%P為節點有功功率
Q=zeros（x1）;%Q為節點無功功率
S=ones（x1）+j*ones（x1）;%S為節點復功率
I=zeros（y1）;%I為起始節點編號矩陣
J=zeros（y1）;%J為終止節點編號矩陣
Rij=zeros（y1）;%R為線路電阻
Xij=zeros（y1）;%X為線路電抗
Zij=Rij+j*Xij;??%Yij為線路阻抗
Y=zeros（x）;??????%Y為n階節點導納方陣
G=zeros（x）;??????%G為n階節點電導方陣
B=zeros（x）;??????%B為n階節點電納方陣
B0=zeros（y1）;?%B0為n*1階線路對地電納值
RT=zeros（y1）;%RT為ij支路y（?矩陣zhilu的行數）*1階變壓器電阻
XT=zeros（y1）;%XT為ij支路y*1階變壓器電抗
ZT=RT+j*XT;%求變壓器阻抗
KT=zeros（y1）;?%K為ij支路y*1階變壓器變比，若k=0表示無變壓器K=1則為標準變比，k?不等于1為非標準變比
%------------------------------矩陣賦初值：
TYPE=point（:1）;%將point矩陣的第一列賦給TYPE，以下類似
U=point（:2）;
a=point（:3）;
P=point（:4）;
Q=point（:5）;
I=zhilu（:1）;
J=zhilu（:2）;
Rij=zhilu（:3）;
Xij=zhilu（:4）;
Zij=Rij+j*Xij;
B0=zhilu（:5）;
RT=zhilu（:6）;
XT=zhilu（:7）;
ZT=RT+j*XT;
KT=zhilu（:8）;
%------------------------------求節點導納矩陣Y
for?m=1:y??????%求Y中非對角元元素Yij
????if?KT（m）==0%若無變壓器則Yij直接為線路阻抗分之一取負值．
????????Y（I（m）J（m））=-1/Zij（m）;
????????Y（J（m）I（m））=-1/Zij（m）;
????else?%有變壓器時，Ｙij為線路阻抗乘以ＫＴ后分之一再取負值
??????????Y（I（m）J（m））=-1/（KT（m）*ZT（m））;
??????????Y（J（m）I（m））=-1/（KT（m）*ZT（m））;
???end
end
for?m=1:x????%求Y中的Ｙii
????for?n=1:y
???????if?KT（n）==0%無變壓器時Ｙii為Ｙij加上線路對地電導乘j
?????????if（I（n）==m|J（n）==m）
????????????Y（mm）=Y（mm）-Y（I（n）J（n））+j*B0（n）/2;
?????????end
??????elseif?I（n）==m??%有變壓器時，若支路起始節點為m則變壓器等值模型的對地支路的（1-KT）/KT^2算到Ｉ（m）節點
???????????Y（mm）=Y（mm）-Y（I（n）J（n））+（1-KT（n））/（KT（n）^2）*（1/ZT（n））;
??????elseif?J（n）==m???%有變壓器時，若支路終止節點為m則變壓器等值模型的對地支路的（KT-1）/KT算到J（m）節點
???????????Y（mm）=Y（mm）-Y（I（n）J（n））+（KT（n）-1）/KT（n）*（1/ZT（n））;??
???????else??Y（mm）=Y（mm）;
???????end
???end
end
Y
G=real（Y）;%求Y的實部，節點電導矩陣
B=imag（Y）;%求Y的虛部節點電納矩陣
%-------------計算各節點有功功率不平衡量deltaPi
?k=0;???%k為迭代次數
?kp=0;??%計算P不平衡量deltaPi的收斂標志（0表示不收斂1表示收斂）
?kq=0;??%計算U不平衡量deltaQi的收斂標志（0表示不收斂1表示收斂）
?ph=find（TYPE（:1）==3）;%找出平衡節點編號
?notph=find（TYPE（:1）<3）;%找出非平衡節點編號
?deltaPi=zeros（x-11）;%deltaPi為（x-1）*1階矩陣x即為節點數
?pq=find（TYPE（:1）==1）;%找出PQ節點編號
?pqnum=size（pq）;
?pqnum=pqnum（1）;%求PQ節點個數
?deltaQi=zeros（pqnum1）;%deltaQi為pqnum*1階矩陣
while（（~kq|~kp）&（k<50））?
????k=k+1;
????%------書上式4-45（a）求deltaPi?
????for?m=1:（x-1）
?????????sum1=0;
?????????for?n=1:x
?????????????sum1=sum1+U（notph（m））*U（n）*（G（notph（m）n）*cos（a（notph（m））-a（n））+B（notph（m）n）*sin（a（notph（m））-a（n）））;
?????????end
?????????deltaPi（m）=P（notph（m））-sum1;
????end

???%------書上式4-45（b）求deltaQi
?????for?m=1:pqnum??
?????????sum2=0;
?????????for?n=1:x

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件??????12731??2010-05-19?11:06??PQandNRchaoliulast\NRchaoliu.m

?????文件??????17408??2010-05-19?11:32??PQandNRchaoliulast\NRinput.xls

?????文件??????10307??2010-05-19?11:39??PQandNRchaoliulast\PQchaoliu.m

?????文件??????17408??2010-05-19?11:32??PQandNRchaoliulast\PQinput.xls

?????目錄??????????0??2010-05-19?11:33??PQandNRchaoliulast

-----------?---------??----------?-----??----

????????????????57854????????????????????5