IEEE算例潮流程序

%% 运用牛顿-拉夫逊法解IEEE潮流  
%%%%%%把原始文件分割成两个独立的纯数据文件分别保存为Bus.DAT和Line.DAT
 
clear;
clc;
tic;                        %开始计时
%% 读取Bus.DAT中的数据
load（‘C:\Users\xyphhunj\Desktop\普通潮流\030bus\Bus.DAT‘）;
Node=Bus（:1） ;             %设定节点号
Voltbase=Bus（:3） ;         %电压标准值
Type1=Bus（:6）;             %节点类型
V=Bus（:7）;                 %节点电压表幺值
Angle=Bus（:8）;             %电压相角
Pload=Bus（:9）;             %负荷有功
Qload=Bus（:10）;            %负荷载无功
Pgen=Bus（:11）;             %发电有功
Qgen=Bus（:12）;             %发电无功
Vset=Bus（:14）;             %电压初始值
Qgmax=Bus（:15）;            %Q的最大值
Qgmin=Bus（:16）;            %Q的最小值
Qsh=Bus（:18）;              %并联电容电纳为正，并联电抗电纳为负
base=100;                   %功率标准值
%%  读取Line.DAT中的数据
load（‘C:\Users\xyphhunj\Desktop\普通潮流\030bus\Line.DAT‘）;
II=Line（:1）;                  
JJ=Line（:2）;                  %两端点
Type2=Line（:5）;               %线路类型
R=Line（:6）;                   %两点间电阻
X=Line（:7）;                   %两点间电抗
B0=Line（:8）/2;                %线路两端电纳
K=Line（:14）;                  % 变压器变比
Kmax=Line（:16）;
Kmin=Line（:17）; 
Kset=Line（:18）;               %变比要求
%% 开始求导纳 Y
num=size（Bus1）;
y0=zeros（num） ;                  %两点间的支路导纳
yt=zeros（num） ;                  %变压器的对地导纳
for x=1:size（Line1） 
    switch Type2（x）             %比较线路类型  
        case 2                  %含变压器的
            y0（II（x）JJ（x））=1./（K（x）.*（R（x）+1i*X（x）））+y0（II（x）JJ（x）） ;                   
            y0（JJ（x）II（x））=y0（II（x）JJ（x））    ;               
            yt（II（x）JJ（x））=（1-K（x））./（（R（x）+1i*X（x））.*K（x）^2）+yt（II（x）JJ（x））; 
            yt（JJ（x）II（x））=（K（x）-1）./（（R（x）+1i*X（x））.*K（x））+yt（JJ（x）II（x））;    
        case 1                 %不含变压器的
            y0（II（x）JJ（x））=1/（R（x）+1i*X（x））+y0（II（x）JJ（x））;                 
            y0（JJ（x）II（x））=y0（II（x）JJ（x））;                  
    end
end
Y=zeros（num）;                  %导纳矩阵
clear x
for x=1:num
    Y（xx）=sum（y0（x:）+yt（x:））+1i*Qsh（x） ;     %求出对角导纳（不含B/2）
end
Y=Y-y0;
clear x
for x=1:size（Line1）
    Y（II（x）II（x））=Y（II（x）II（x））+1i*B0（x） ;
    Y（JJ（x）JJ（x））=Y（JJ（x）JJ（x））+1i*B0（x） ;          %加B
end
clear x                           % Y导纳矩阵结束
G=real（Y）;
B=imag（Y）;
%% 
D=ones（2*num1）;
%给定节点电压值，a为电压相角，u为电压幅值
for m=1:num
    a（m）=0;          %初始设定角度
    u（m）=Bus（m14）;  %初始迭代电压值
    e（m）=u（m）*cos（a（m））;
    f（m）=u（m）*sin（a（m））;  %电压实部和虚部
end
t=0;                    %设定迭代次数变量t
while （max（abs（D））>1e-8）

   for m=1:num
      Ps（m1）=（Pgen（m1）-Pload（m1））/base;
      Qs（m1）=（Qgen（m1）-Qload（m1））/base;      %节点注入有功和注入无功                                                                            
   end
%计算功率不平衡量有功不平衡量为dP，无功不平衡量为dQ电压平方的不平衡量为dUU
   sum1=zeros（1num）;
   sum2=zeros（1num）;
   sum3=zeros（1num）;
   for m=1:num
       for n=1:num
          sum1（m）=sum1（m）+G（mn）*e（n）-B（mn）*f（n）;
          sum2（m）=sum2（m）+G（mn）*f（n）+B（mn）*e（n）;     
       end
       sum3（m）=e（m）^2+f（m）^2;
   end
   for m=1:nu