含风电场的电力系统潮流程序matlab

function reactivepower
clear
baseMVA=100;                        %功率基值
%%读Data1中数据 
load Data1.txt
Bus=Data1（:1）;                        %节点号
Vtype=Data1（:2）;                    %节点类型
Pload=Data1（:3）;                    %负载有功
Qload=Data1（:4）;                    %负载无功
Vset=Data1（:8）;                      %母线电压幅值
Pgen=Data1（:14）;                  %发电机发出有功
Qgen=Data1（:15）;                  %发电机发出无功
Qsh=Data1（:16）;                  %并联电容电纳标幺值
%%读Data2中数据
load Data2.txt
II=Data2（:1）;                      
JJ=Data2（:2）;                        %支路端点号
R=Data2（:3）;                         %两点间电阻
X=Data2（:4）;                         %两点间电抗
B=Data2（:5）/2;                     %线路对地电纳
K=Data2（:9）;                         %变压器非标准变压比
Ltype=Data2（:14）;                %线路类型
%%求导纳矩阵Y
y1=zeros（14）;
y2=zeros（14）;
y3=zeros（14）;
lin=length（II）;                     %支路数
for x=1:lin
    switch Ltype（x）
        case 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*B（x）;
            y3（JJ（x）II（x））=i*B（x）;
        case 2
            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（14）;
for x=1:14
    Y（xx）=sum（y1（x:））+sum（y2（x:））+sum（y3（x:））+i*Qsh（x）;
end
clear x;
Y=Y-y1;
G=real（Y）;
B=imag（Y）;
%%设电压初值
U=Vset;
e=real（U）;
f=imag（U）;
%%
Ps=zeros（114）;
Qs=zeros（114）;
D=ones（261）;
for x=1:14
    Ps（x）=（Pgen（x）-Pload（x））/baseMVA;
    Qs（x）=（Qgen（x）-Qload（x））/baseMVA;
end
clear x;
N=0;
Jacbi=zeros（26）;
while max（abs（D））>0.000001
       for x=2:14                                            %节点功率及电压不平衡量
        switch Vtype（x）
            case 1                                           %PQ节点
                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                                           %PV节点
                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:14                                                   %求雅克比矩阵
      for J=2:14
          if I~=J                                              %非对角元素
              Jacbi（（2*I-3）（2*J-3））=B