%function?[YUPQdeltaSijaSSijSjisumdeltaS]=PQchaoliu（）
clc?
close?all
clear?all
%function?[Ua]=PQchaoliu（）
%電力系統(tǒng)潮流計(jì)算程序；
%輸出：U——節(jié)點(diǎn)電壓，P--節(jié)點(diǎn)有功，Q--節(jié)點(diǎn)無(wú)功，deltaSij--支路功率損耗，
%Sij--從節(jié)點(diǎn)i流向節(jié)點(diǎn)j的功率S--節(jié)點(diǎn)復(fù)功率sumdeltaS--網(wǎng)絡(luò)總損耗
%輸入?yún)?shù)：point為節(jié)點(diǎn)信息矩陣，zhilu為支路信息矩陣；
[x]=xlsread（‘shiyanshuju.xlsx‘‘A2:A2‘）;%從exel中讀取節(jié)點(diǎn)數(shù)x
[y]=xlsread（‘shiyanshuju.xlsx‘‘B2:B2‘）;%從exel中讀取支路數(shù)y
e=xlsread（‘shiyanshuju.xlsx‘‘B4:B4‘）;%誤差要求
[point]=xlsread（‘shiyanshuju.xlsx‘‘D3:H100‘）;%從exel中讀取節(jié)點(diǎn)信息矩陣
[zhilu]=xlsread（‘shiyanshuju.xlsx‘‘J3:Q100‘）;%從exel中讀取支路信息矩陣
TYPE=zeros（x1）;%TYPE為節(jié)點(diǎn)類型矩陣
U=zeros（x1）;%U為節(jié)點(diǎn)電壓矩陣
a=zeros（x1）;%a為節(jié)點(diǎn)電壓相角矩陣******此處應(yīng)當(dāng)為弧度制
P=zeros（x1）;%P為節(jié)點(diǎn)有功功率
Q=zeros（x1）;%Q為節(jié)點(diǎn)無(wú)功功率
I=zeros（y1）;%I為起始節(jié)點(diǎn)編號(hào)矩陣
J=zeros（y1）;%J為終止節(jié)點(diǎn)編號(hào)矩陣
Rij=zeros（y1）;%R為線路電阻
Xij=zeros（y1）;%X為線路電抗
Zij=Rij+j*Xij;??%Yij為線路阻抗
Y=zeros（x）;??????%Y為n階節(jié)點(diǎn)導(dǎo)納方陣
G=zeros（x）;??????%G為n階節(jié)點(diǎn)電導(dǎo)方陣
B=zeros（x）;??????%B為n階節(jié)點(diǎn)電納方陣
B0=zeros（y1）;?%B0為n*1階線路對(duì)地電納值
RT=zeros（y1）;%RT為ij支路y（?矩陣zhilu的行數(shù)）*1階變壓器電阻
XT=zeros（y1）;%XT為ij支路y*1階變壓器電抗
ZT=RT+j*XT;%求變壓器阻抗
KT=zeros（y1）;?%K為ij支路y*1階變壓器變比，若k=0表示無(wú)變壓器K=1則為標(biāo)準(zhǔn)變比，k?不等于1為非標(biāo)準(zhǔn)變比
%------------------------------矩陣賦初值：
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）;
tic??%程序運(yùn)行時(shí)間開(kāi)始計(jì)時(shí)
%------------------------------求節(jié)點(diǎn)導(dǎo)納矩陣Y
for?m=1:y??????%求Y中非對(duì)角元元素Yij
????if?KT（m）==0%若無(wú)變壓器則Yij直接為線路阻抗分之一取負(fù)值．
????????Y（I（m）J（m））=-1/Zij（m）;
????????Y（J（m）I（m））=-1/Zij（m）;
????else?%有變壓器時(shí)，Ｙij為線路阻抗乘以ＫＴ后分之一再取負(fù)值
??????????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%無(wú)變壓器時(shí)Ｙii為Ｙij加上線路對(duì)地電導(dǎo)的一半乘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??%有變壓器時(shí)，若支路起始節(jié)點(diǎn)為m則變壓器等值模型的對(duì)地支路的（KT-1）/KT算到Ｉ（m）節(jié)點(diǎn)
????????Y（mm）=Y（mm）-Y（I（n）J（n））+（KT（n）-1）/KT（n）*（1/ZT（n））;?
??????elseif?J（n）==m???%有變壓器時(shí)，若支路終止節(jié)點(diǎn)為m則變壓器等值模型的對(duì)地支路的（1-KT）/KT^2算到J（m）節(jié)點(diǎn)
?????????Y（mm）=Y（mm）-Y（I（n）J（n））+（1-KT（n））/（KT（n）^2）*（1/ZT（n））;
???????else??Y（mm）=Y（mm）;
???????end
???end
end
Y;
G=real（Y）;
%-----------------------求B‘矩陣及其逆矩陣B1
B=imag（Y）;%求Y的虛部節(jié)點(diǎn)電納矩陣
?ph=find（TYPE（:1）==3）;%找出平衡節(jié)點(diǎn)編號(hào)
?BB=B;%BB矩陣為中間變量
?BB（:ph）=[];%平衡節(jié)點(diǎn)編號(hào)對(duì)應(yīng)列置空并消去
?BB（ph:）=[];%平衡節(jié)點(diǎn)編號(hào)對(duì)應(yīng)行置空
?B1=BB;
?B1=inv（B1）;%B1求逆后得B1矩陣
?%-----------------------%求B‘‘及其逆矩陣B2
phpv=find（TYPE（:1）>1）;%找出非PQ節(jié)點(diǎn)的編號(hào)
?BB=B;?????%以下與求B‘相同
?BB（:phpv）=[];
?BB（phpv:）=[];
?B2=BB;????
?B2=inv（B2）;%求得B2矩陣
?
?%-------------計(jì)算各節(jié)點(diǎn)有功功率不平衡量deltaPi
?k=0;???%k為迭代次數(shù)
?kp=0;??%計(jì)算P不平衡量deltaPi的收斂標(biāo)志（0表示不收斂1表示收斂）
?kq=0;??%計(jì)算U不平衡量deltaQi的收斂標(biāo)志（0表示不收斂1表示收斂）
?notph=find（TYPE（:1）<3）;%找出非平衡節(jié)點(diǎn)編號(hào)
?deltaPi=zeros（x-11）;%deltaPi為（x-1）*1階矩陣x即為節(jié)點(diǎn)數(shù)
?pq=find（TYPE（:1）==1）;%找出PQ節(jié)點(diǎn)編號(hào)
?pqnum=size（B2）;
?pqnum=pqnum（1）;%求PQ節(jié)點(diǎn)的個(gè)數(shù)（因B1矩陣的維數(shù)等于PQ節(jié)點(diǎn)數(shù)）
?deltaQi=zeros（pqnum1）;%deltaQi為pqnum*1階矩

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----
?????文件??????329484??2019-03-03?19:37??pq分解法\2.jpg
?????文件???????10094??2019-03-04?20:38??pq分解法\PQchaoliu.m
?????文件????????2953??2019-03-04?20:39??pq分解法\outputpq.txt
?????文件???????16255??2019-03-03?11:19??pq分解法\shiyanshuju.xlsx
?????目錄???????????0??2019-03-04?20:38??pq分解法\