基于牛拉法的电力系统潮流计算程序设计

clear
%手动导入节点数据后保存为.mat文件
%所有数据均为标幺值
load（‘matlab.mat‘）
%Bus：节点号
%Vtype：节点类型
%Pload：负荷有功
%Qload：负荷无功
%Pgen：发电机有功
%Qgen：发电机无功
%Vset：节点电压
%Qsh：电纳标幺值
%II、JJ：支路端点号
%R：支路电阻
%X：支路电抗
%B：对地电纳/2
%K：变压器非标准变压比
%Ltype：是否有变比支路
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求导纳矩阵Y
y1=zeros（30）;
y2=zeros（30）;
y3=zeros（30）;
lin=length（II）;%求支路数
for x=1:lin
    switch Ltype（x）
        case 1
            y1（II（x）JJ（x））=1/（R（x）+1i*X（x））;
            y1（JJ（x）II（x））=y1（II（x）JJ（x））;
            y3（II（x）JJ（x））=1i*B（x）;
            y3（JJ（x）II（x））=1i*B（x）;
        case 2
            y1（II（x）JJ（x））=1/（（R（x）+1i*X（x））*K（x））;
            y1（JJ（x）II（x））=y1（II（x）JJ（x））;
            y2（II（x）JJ（x））=（1-K（x））/（（R（x）+1i*X（x））*K（x）^2）;%IIII=y1+y2           
            y2（JJ（x）II（x））=（K（x）-1）/（（R（x）+1i*X（x））*K（x））;%JJJJ=y1+y2
    end
end
clear x;
Y=zeros（30）;%定义导纳矩阵
for x=1:30%求对角线元素
    Y（xx）=sum（y1（x:））+sum（y2（x:））+sum（y3（x:））+1i*Qsh（x）;
end
clear x;
Y=Y-y1;%非对角线元素
G=real（Y）;
B=imag（Y）;
%%%%%%%%%%%%%%%%%%
%设电压初值
U=Vset;
e=real（U）;%实部
f=imag（U）;
%
Ps=zeros（301）;%注入有功
Qs=zeros（301）;%注入无功
D=ones（301）;
for x=1:30
    Ps（x）=Pgen（x）-Pload（x）;
    Qs（x）=Qgen（x）-Qload（x）;
end
clear x;
N=0;
%24*2+5+2=55
Jacbi=zeros（55）;
while max（abs（D））>0.0000000001
for x=2:30%节点功率及电压不平衡量
    switch Vtype（x）
        case 1
            D（2*x-3）=Ps（x）-e（x）*（G（x:）*e-B（x:）*f）-f（x）*（G（x:）*f+B（x:）*e）;
            D（2*x-2）=Qs（x）-f（x）*（G（x:）*e-B（x:）*f）+e（x）*（G（x:）*f+B（x:）*e）;
        case 2
            D（2*x-3）=Ps（x）-e（x）*（G（x:）*e-B（x:）*f）-f（x）*（G（x:）*f+B（x:）*e）;
            D（2*x-2）=Vset（x）.*Vset（x）-（e（x）.^2+f（x）.^2）;
    end
end
clear m
for I=2:30
    for J=2:30
        if I~=J%非对角元素
            Jacbi（（2*I-3）（2*J-3））=B（IJ）*e（I）-G（IJ）*f（I）;
            Jacbi（（2*I-3）（2*J-2））=-（G（IJ）*e（I）+B（IJ）*f（I））;
          s