%function?[YUPQdeltaSijaSSijSjisumdeltaS]=PQchaoliu（）
clc?
close?all
clear?all
%function?[Ua]=PQchaoliu（）
%電力系統潮流計算程序；
%輸出：U——節點電壓，P--節點有功，Q--節點無功，deltaSij--支路功率損耗，
%Sij--從節點i流向節點j的功率S--節點復功率sumdeltaS--網絡總損耗
%輸入參數：point為節點信息矩陣，zhilu為支路信息矩陣；
[x]=xlsread（‘shiyanshuju.xlsx‘‘A2:A2‘）;%從exel中讀取節點數x
[y]=xlsread（‘shiyanshuju.xlsx‘‘B2:B2‘）;%從exel中讀取支路數y
e=xlsread（‘shiyanshuju.xlsx‘‘B4:B4‘）;%誤差要求
[point]=xlsread（‘shiyanshuju.xlsx‘‘D3:H100‘）;%從exel中讀取節點信息矩陣
[zhilu]=xlsread（‘shiyanshuju.xlsx‘‘J3:Q100‘）;%從exel中讀取支路信息矩陣
TYPE=zeros（x1）;%TYPE為節點類型矩陣
U=zeros（x1）;%U為節點電壓矩陣
a=zeros（x1）;%a為節點電壓相角矩陣******此處應當為弧度制
P=zeros（x1）;%P為節點有功功率
Q=zeros（x1）;%Q為節點無功功率
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）;
tic??%程序運行時間開始計時
%------------------------------求節點導納矩陣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則變壓器等值模型的對地支路的（KT-1）/KT算到Ｉ（m）節點
????????Y（mm）=Y（mm）-Y（I（n）J（n））+（KT（n）-1）/KT（n）*（1/ZT（n））;?
??????elseif?J（n）==m???%有變壓器時，若支路終止節點為m則變壓器等值模型的對地支路的（1-KT）/KT^2算到J（m）節點
?????????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的虛部節點電納矩陣
?ph=find（TYPE（:1）==3）;%找出平衡節點編號
?BB=B;%BB矩陣為中間變量
?BB（:ph）=[];%平衡節點編號對應列置空并消去
?BB（ph:）=[];%平衡節點編號對應行置空
?B1=BB;
?B1=inv（B1）;%B1求逆后得B1矩陣
?%-----------------------%求B‘‘及其逆矩陣B2
phpv=find（TYPE（:1）>1）;%找出非PQ節點的編號
?BB=B;?????%以下與求B‘相同
?BB（:phpv）=[];
?BB（phpv:）=[];
?B2=BB;????
?B2=inv（B2）;%求得B2矩陣
?
?%-------------計算各節點有功功率不平衡量deltaPi
?k=0;???%k為迭代次數
?kp=0;??%計算P不平衡量deltaPi的收斂標志（0表示不收斂1表示收斂）
?kq=0;??%計算U不平衡量deltaQi的收斂標志（0表示不收斂1表示收斂）
?notph=find（TYPE（:1）<3）;%找出非平衡節點編號
?deltaPi=zeros（x-11）;%deltaPi為（x-1）*1階矩陣x即為節點數
?pq=find（TYPE（:1）==1）;%找出PQ節點編號
?pqnum=size（B2）;
?pqnum=pqnum（1）;%求PQ節點的個數（因B1矩陣的維數等于PQ節點數）
?deltaQi=zeros（pqnum1）;%deltaQi為pqnum*1階矩

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件??????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分解法\