基于MATLAB的潮流计算程序

%本程序的功能是用牛顿――拉夫逊法进行潮流计算
% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳
%         5、支路的变比；6、支路首端处于K侧为1，1侧为0
% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值
%         4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量
%         6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；
%            3为PV节点；
clear;
n=10;%input('请输入节点数：n=');
nl=10;%input('请输入支路数：nl=');
isb=1;%input('请输入平衡母线节点号：isb=');
pr=0.00001;%input('请输入误差精度：pr=');
B1=[1   2   0.03512+0.08306i   0.13455i    1          0;
    2   3   0.0068+0.18375i    0           1.02381    1;
    1   4   0.05620+0.13289i   0.05382i    1          0;
    4   5   0.00811+0.24549i   0           1.02381    1;
    1   6   0.05620+0.13289i   0.05382i    1          0;
    4   6   0.04215+0.09967i   0.04037i    1          0;
    6   7   0.0068+0.18375i    0           1.02381    1;
    6   8   0.02810+0.06645i   0.10764i    1          0;
    8   10  0.00811+0.24549i   0           1          1;
    8   9   0.03512+0.08306i   0.13455i    1          0]
B2=[0     0                1.1      1.1      0         1;
    0     0                1        0        0         2;
    0     0.343+0.21256i   1        0        0         2;
    0     0                1        0        0         2;
    0     0.204+0.12638i   1        0        0         2;
    0     0                1        0        0         2;
    0     0.306+0.18962i   1        0        0         2;
    0     0                1        0        0         2;

    0.5   0                1.1      1.1      0         3;
    0     0.343+0.21256i   1        0        0         2]


;%input('请输入各节点参数形成的矩阵： B2=');
Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);
% % %---------------------------------------------------
for i=1:nl %支路数
    if B1(i,6)==0 %左节点处于1侧
       p=B1(i,1);q=B1(i,2);
    else                                %左节点处于K侧
       p=B1(i,2);q=B1(i,1);
    end
    Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元
    Y(q,p)=Y(p,q);                      %非对角元
    Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元K侧
    Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧 
end
%求导纳矩阵
disp('导纳矩阵 Y=');
disp(Y)
%----------------------------------------------------------
G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部  
for i=1:n %给定各节点初始电压的实部和虚部    
    e(i)=real(B2(i,3));
    f(i)=imag(B2(i,3));
    V(i)=B2(i,4); %PV节点电压给定模值 
end
for i=1:n %给定各节点注入功率    
    S(i)=B2(i,1)-B2(i,2);  %i节点注入功率SG-SL