matlab潮流计算PQ解耦法14节点计算程序

本程序通过matlab采用PQ解耦法对电力系统潮流进行计算，可以实现14节点的计算

%      fir	fin        r          x           b     k;
branch=...	
       [1	3	0.082     0.427	0.028	1;
        1	5	0.163   0.754	0.014	1;
        3	4	0.145   0.581	0.021	1;
        4	5	0.104   0.518	0.018	1;
        2	5	0.031   0.248	0  	0.95];

%     no type    p        q       v        a;
bus=...
    [1    0  1            1   1.05   0;
     2    1  0.35      0.35   1.0    0;
     3    2  0.22     0.14    1.0   0; 
     4    2  0.18     0.09    1.0   0;
     5    2  0.27     0.13    1.0   0];
 
Yc=[000.04*i00];

baseMVA=100;               %功率基值%%读Data1中数据 
NO=bus（:1）;
Type=bus（:2）;             %节点类型
P=bus（:3）;                %有功
Q=bus（:4）;                %无功
U=bus（:5）;                %电压设定点
a=bus（:6）;                %并联电容电纳标幺值
II=branch（:1）;                     
JJ=branch（:2）;            %支路端点号
R=branch（:3）;             %两点间电阻
X=branch（:4）;             %两点间电抗
BD=branch（:5）;           %线路对地电纳
K=branch（:6）;             %变压器非标准变压比
y1=zeros（5）;
y2=zeros（5）;
y3=zeros（5）;
lin=length（II）;                     %支路数
for x=1:1:lin
    if（K（x）==1）
            y1（II（x）JJ（x））=1/（R（x）+i*X（x））;
            y1（JJ（x）II（x））=y1（II（x）JJ（x））;
            y3（II（x）JJ（x））=i*BD（x）;
            y3（JJ（x）II（x））=i*BD（x）;
    else
            y1（II（x）JJ（x））=1/（（R（x）+i*X（x））*K（x））;
            y1（JJ（x）II（x））=y1（II（x）JJ（x））;
            y2（II（x）JJ（x））=（1-K（x））/（（R（x）+i*X（x））*K（x）^2）;
            y2（JJ（x）II（x））=（K（x）-1）/（（R（x）+i*X（x））*K（x））;
    end  
end
clear x
Y=zeros（5）;
for x=1:1:5
    Y（xx）=sum（y1（x:））+sum（y2（x:））+sum（y3（x:））+Yc（x）;
end
clear x;
Y=Y-y1;
B=imag（Y）;
G=real（Y）;
ph=find（Type（:1）==0）; 
BB=B;  
BB（:ph）=[]; 
BB（ph:）=[]; 
B1=BB;
phpv=find（Type（:1）<2）;         
BB=B;   
BB（:phpv）=[]; 
BB（phpv:）=[];        
B2=BB;
count=0;
E=0.00001;
DPI=zeros（41）;
DQI=zeros（31）;
nph=find（Type（:1）>0）;
pq=find（Type（:1）==2）;
PD1=1;
PD2=1;
while（（PD1>E）||（PD2>E））
    for m=1:1:4
        sum1=0;
        for n=1:1:5
        sum1=sum1+U（nph（m）1）*U（n1）*（G（nph（m）n）*cos（a（nph（m）1）-a（n1））+...
		B（nph（m）n）*sin（a（nph（m）1）-a（n1）））;
        end
        DPI（m1）=P（nph（m）1）-sum1;
    end
    DAA=（-inv（B1）*（DPI（1:4）./U（2:5）））./U（2:5）;
    a（2:5）=a（2:5）+DAA（1:4）;
    PD1=max（abs（DPI（1:4）./U（2:5）））;
    for m=1:1:3
        sum2=0;
        for n=1:1:5
        sum2=sum2+U（pq（m）1）*U（n1）*（G（pq（m）n）*sin（abs（a（pq（m）1）-a（n1）））-...
		B（pq（m）n）*cos（abs（a（pq（m）1）-a（n1））））;
        end
        DQI（m1）=Q（pq（m）1）-sum2;
    end
    Uq=U;
    Uq（phpv:）=[];
    Upq=Uq;
    DUI=（-inv（B2）*（DQI（1:3）./Upq（1:3）））;
    for m=1:1:3
        U（pq（m）1）=U（pq（m）1）+DUI（m1）;
    end
    PD2=max（abs（DQI（1:3）./Upq（1:3）））;
    count=count+1;
end
S0=U（1）*（conj（Y（1:））*conj（U））;                            %平衡节点功率
S1=zeros（51）;                                            %始端功率 
S2=zeros（51）;                             

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3917  2014-12-17 13:18  Power system\bandan3_command.asv

     文件       3924  2014-12-17 21:52  Power system\bandan3_command.m

     文件        784  2014-12-17 11:39  Power system\branchpq.m

     文件        691  2014-12-17 11:39  Power system\buspq.m

     文件     315261  2014-12-14 23:49  Power system\稳态课设题目相关.rar

     目录          0  2015-12-17 16:29  Power system

----------- ---------  ---------- -----  ----

               324577                    6