电力系统稳态潮流计算程序PQ和NR法-PQandNRchaoliulast.rar

电力系统稳态潮流计算程序PQ和NR法-PQandNRchaoliulast.rar
RT。程序根据中国电力出版社的《电力系统稳态分析》新版编的，里面的归算方法和变压器模型也是这本书上讲的，可以针对不同电压等级的系统进行潮流计算，提供两个方法，PQ法和NR法，NR法相对精确一点，但速度稍慢，数据采用表格输入， PQandNRchaoliulast.rar

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