matlab潮流程序牛顿拉夫逊+PQ分解法有算例

matlab潮流程序（牛顿拉夫逊+PQ分解法）有算例，有注解，亲测可用，小规模节点下的结果与BPA结果一致，有兴趣的朋友可以对程序进行改进（比如PV节点转PQ节点处理等）

%function [YUPQdeltaSijaSSijSjisumdeltaS]=NRchaoliu（）
function [Ua]=NRchaoliu（）
%电力系统潮流计算程序；
%输出：U——节点电压，P--节点有功，Q--节点无功，deltaSij--支路功率损耗，
%Sij--从节点i流向节点j的功率S--节点复功率sumdeltaS--网络总损耗
%输入参数：point为节点信息矩阵，zhilu为支路信息矩阵；
[x]=xlsread（‘NRinput.xls‘‘A2:A2‘）;%从exel中读取节点数x
[y]=xlsread（‘NRinput.xls‘‘B2:B2‘）;%从exel中读取支路数y
e=xlsread（‘NRinput.xls‘‘B4:B4‘）;%误差要求
[point]=xlsread（‘NRinput.xls‘‘C3:H50‘）;%从exel中读取节点信息矩阵 zys
[zhilu]=xlsread（‘NRinput.xls‘‘J3:Q50‘）;%从exel中读取支路信息矩阵
sizeHL=size（zhilu）;
zhilu0=zeros（sizeHL（1）2）;%增加2列存储首末新节点编号zys
zhilu=[zhilu zhilu0];%zys
%-------节点优化编号定义start  zys-------------------
nSW = 0;                   % nSW为平衡节点个数 
nPV = 0;                   % nPV为PV节点个数 
nPQ = 0;                   % nPQ为PQ节点个数 
for lx=1:x
    type0=point（lx2）;
    if type0 == 3      
        nSW = nSW + 1;     
        SW（nSW:）=point（lx:）;% 计算并储存平衡节点   
    elseif type0 == 2       
        nPV = nPV +1;      
        PV（nPV:）=point（lx:）;% 计算并储存PV节点   
    else       
        nPQ = nPQ + 1;     
        PQ（nPQ:）=point（lx:）;% 计算并储存PQ节点   
    end 
end
point=[PQ;PV;SW];            % 对bus矩阵按PQ、PV、平衡节点的顺序重新排序
nodenum=[[1:x]‘ point（:1）];% 生成新旧节点对照表
for hh=1:y     
    for ll=1:2         
        for k=1:x             
            if zhilu（hhll）==nodenum（k2）
                zhilu（hh8+ll）=nodenum（k1）;                
                break             
            end         
        end     
    end 
end    
%-------节点优化编号定义end zys-------------------


TYPE=zeros（x1）;%TYPE为节点类型矩阵
U=zeros（x1）;%U为节点电压矩阵
a=zeros（x1）;%a为节点电压相角矩阵
P=zeros（x1）;%P为节点有功功率
Q=zeros（x1）;%Q为节点无功功率
S=ones（x1）+j*ones（x1）;%S为节点复功率
I=zeros（y1）;%I为起始节点编号矩阵
J=zeros（y1）;%J为终止节点编号矩阵
Rij=zeros（y1）;%R为线路电阻
Xij=zeros（y1）;%X为线路电抗
Zij=Rij+j*Xij;  %Yij为线路阻抗
Y=zeros（x）;      %Y为n阶节点导纳方阵
G=zeros（x）;      %G为n阶节点电导方阵
B=zeros（x）;      %B为n阶节点电纳方阵
B0=zeros（y1）; %B0为n*1阶线路对地电纳值
RT=zeros（y1）;%RT为ij支路y（ 矩阵zhilu的行数）*1阶变压器电阻
XT=zeros（y1）;%XT为ij支路y*1阶变压器电抗
ZT=RT+j*XT;%求变压器阻抗
KT=zeros（y1）; %K为ij支路y*1阶变压器变比，若k=0表示无变压器K=1则为标准变比，k 不等于1为非标准变比
%------------------------------矩阵赋初值：

%TYPE=point（:1）;%将point矩阵的第一列赋给TYPE，以下类似
TYPE=point（:2）;%将point矩阵的第一列赋给TYPE，以下类似
%U=point（:2）;a=point（:3）;P=point（:4）;Q=point（:5）;
U=point（:3）;
a=point（:4）;
P=point（:5）;
Q=point（:6）;
%I=zhilu（:1）;J=zhilu（:2）;
I=zhilu（:9）;
J=zhilu（:10）;
Rij=zhilu（:3）;
Xij=zhilu（:4）;
Zij=Rij+j*Xij;
B0=zhilu（:5）;
RT=zhilu（:6）;
XT=zhilu（:7）;
ZT=RT+j*XT;
KT=zhilu（:8）;
tic  %程序运行时间开始计时
%------------------------------求节点导纳矩阵Y
for m=1:y      %求Y中非对角元元素Yij
    if KT（m）==0%若无变压器则Yij直接为线路阻抗分之一取负值．
        Y（I（m）J（m））=-1/Zij（m）;
        Y（J（m）I（m））=-1/Zij（m）;
    else %有变压器时，Ｙij为线路阻抗乘以ＫＴ后分之一再取负值
          Y（I（m）J（m））=-1/（KT（m）*ZT（m））;
          Y（J（m）I（m））=-1/（KT（m）*ZT（m））;
   end
end
for m=1:x    %求Y中的Ｙii
    for n=1:y
       if KT（n）==0%无变压器时Ｙii为Ｙij加上线路对地电导乘j
         if（I（n）==m|J（n）==m）
            Y（mm）=Y（

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

     文件      21671  2017-09-19 15:00  outputNR.txt

     文件      32256  2017-09-19 14:34  NRinput.xls

     文件      15553  2017-09-19 14:22  NRchaoliu.m

     文件       9174  2017-09-19 12:38  outputpq.txt

     文件      11214  2017-09-18 16:10  PQchaoliu.m

     文件      27648  2017-09-15 01:05  PQinput.xls

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

               117516                    6