极坐标牛顿法潮流计算的matlab通用程序

极坐标牛顿法潮流计算的matlab通用程序，精度较高 可以通用 收敛

clear
clc
%极坐标形式的牛顿法潮流计算程序
%B1矩阵为支路参数矩阵，其中1、2：分别为支路首端、末端节点编号；3：支路阻抗；4：接地导纳；5：支路变比；6：支路首端处于K侧为11侧为0；7：单侧对地导纳
%B2矩阵为节点参数矩阵，其中1：该节点的发电机功率；2：该节点负荷功率；3：节点电压初始幅值；4：节点电压初始相位；5：PV节点的电压给定值；6：节点所接的无功补偿设备容量；
%7：节点分类编号：1为平衡节点，2为PQ节点，3为PV节点
n=input（‘请输入节点数:n=‘）;
nl=input（‘请输入支路数:nl=‘）;
isb=input（‘请输入平衡节点编号:isb=‘）;
pr=input（‘请输入误差精度:pr=‘）;
B1=input（‘请输入支路参数矩阵:B1=‘）;
B2=input（‘请输入节点参数矩阵:B2=‘）;
%创建矩阵
Y=zeros（n）;
V=zeros（n1）;
for i=1:nl
    if B1（i6）==0 %支路不含变压器
        p=B1（i1）;q=B1（i2）;
        Y（pq）=Y（pq）-1./B1（i3）;%导纳矩阵的非对角元
        Y（qp）=Y（pq）;
        Y（pp）=Y（pp）+1./B1（i3）+0.5*B1（i4）+B1（i7）;%导纳矩阵的对角元
        Y（qq）=Y（qq）+1./B1（i3）+0.5*B1（i4）+B1（i7）;
    else 
        p=B1（i2）;q=B1（i1）;%支路含有变压器
        Y（pq）=Y（pq）-1./（B1（i3）*B1（i5））;
        Y（qp）=Y（pq）;
        Y（pp）=Y（pp）+1./B1（i3）+B1（i7）; %pi型等效电路
        Y（qq）=Y（qq）+1./（（B1（i5）^2*B1（i3）））+B1（i7）;
    end
end
%分解导纳的实部虚部
G=real（Y）;B=imag（Y）;
for i=1:n
    theta（i1）=B2（i4）;
    V（i1）=B2（i3）;
end
%求出两节点的相位差
Times=0;
%求出PQ节点的个数
m=0;
for i=1:n
    if B2（i7）==2
        m=m+1;
    end
end
M=ones（n-1+m1）;
while max（abs（M））> pr
%形成潮流方程
P=zeros（n1）;
Q=zeros（n1）;
for i=1:n
    for j=1:n
        P（i1）=V（i）*sum（（G（i:）.*cos（theta（i）-theta（:））‘+B（i:）.*sin（theta（i）-theta（:））‘）*V（:））;
        Q（i1）=V（i）*sum（（G（i:）.*sin（theta（i）-theta（:））‘-B（i:）.*cos（theta（i）-theta（:））‘）*V（:））;
    end
end
%求出功率的不平衡量
dP=zeros（n-11）;
dQ=zeros（m1）;
for i=1:n-1
    dP（i1）=real（B2（i1）-B2（i2））-P（i）;
end
for i=1:m
    dQ（i1）=imag（B2（i1）-B2（i2））-Q（i）;
end
H=zeros（n-1n-1）;
N=zeros（n-1m）;
K=zeros（mn-1）;
L=zeros（mm）;
%求出雅克比矩阵元素的表达式
for i=1:n-1
    for j=1:n-1
        if i~=j
            H（ij）=-V（i）*V（j）*（G（ij）*sin（theta（i）-theta（j））-B（ij）*cos（theta（i）-theta（j）））; 
        else
            H（ii）=V（i）^2*B（ii）+Q（i）;
        end
    end
end
for i=1:n-1
   for j=1:m
       if i~=j
           N（ij）=-V（i）*V（j）*（G（ij）*cos（theta（i）-theta（j））+B（ij）*sin（theta（i）-theta（j）））;
       else
           N（ii）=-V（i）^2*G（ii）-P（i）; 
       end
   end
end
for i=1:m
    for j=1:n-1
        if i~=j
            K（ij）=V（i）*V（j）*（G（ij）*cos（theta（i）-theta（j））+B（ij）*sin（theta（i）-theta（j）））;
        else 
            K（ii）=V（i）^2*G（ii）-P（i）;
        end
    end
end
for i=1:m
    for j=1:m
        if i~=j
            L（ij）=-V（i）*V（j）*（G（ij）*sin（theta（i）-theta（j））-B（ij）*cos（theta（i）-theta（j）））;
        else
            L（ii）=V（i）^2*B（ii）-Q（i）;
        end
    end
end
Jacbi=-[HN;KL];
Vd=diag（V）;
A=[dP;dQ];
M=Jacbi\A;
for i=1:n-1
    dtheta（i1）=M（i）;
end
for i=1:n-1
    theta（i1）=theta（i）+dtheta（i）;
end
for i=1:m
    dV（i1）=V（i）*M（i+n-1）;
end
for i=1:m
    V（i）=V（i）+dV（i）;
end
for i=1:n-1
    w（i1）=theta（i）*180/pi;
end
Times=Times+1;
end
Y
dP
dQ
H
N
K
L
Jacbi
Vd
dtheta
dV
Times
w
V

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

     文件       3239  2015-04-18 19:29  PolarCoordinates - 副本.m

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

                 3239                    1