matlab实现电力系统潮流计算-PQ分解法

基于matlab开发，有详细备注，和案例图，其中案例为王锡凡主编的《现代电力系统分析》

%function [YUPQdeltaSijaSSijSjisumdeltaS]=PQchaoliu（）
clc 
close all
clear all
%function [Ua]=PQchaoliu（）
%电力系统潮流计算程序；
%输出：U——节点电压，P--节点有功，Q--节点无功，deltaSij--支路功率损耗，
%Sij--从节点i流向节点j的功率S--节点复功率sumdeltaS--网络总损耗
%输入参数：point为节点信息矩阵，zhilu为支路信息矩阵；
[x]=xlsread（‘shiyanshuju.xlsx‘‘A2:A2‘）;%从exel中读取节点数x
[y]=xlsread（‘shiyanshuju.xlsx‘‘B2:B2‘）;%从exel中读取支路数y
e=xlsread（‘shiyanshuju.xlsx‘‘B4:B4‘）;%误差要求
[point]=xlsread（‘shiyanshuju.xlsx‘‘D3:H100‘）;%从exel中读取节点信息矩阵
[zhilu]=xlsread（‘shiyanshuju.xlsx‘‘J3:Q100‘）;%从exel中读取支路信息矩阵
TYPE=zeros（x1）;%TYPE为节点类型矩阵
U=zeros（x1）;%U为节点电压矩阵
a=zeros（x1）;%a为节点电压相角矩阵******此处应当为弧度制
P=zeros（x1）;%P为节点有功功率
Q=zeros（x1）;%Q为节点无功功率
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，以下类似
U=point（:2）;
a=point（:3）;
P=point（:4）;
Q=point（:5）;
I=zhilu（:1）;
J=zhilu（:2）;
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（mm）-Y（I（n）J（n））+j*B0（n）/2;
         end
      elseif I（n）==m  %有变压器时，若支路起始节点为m则变压器等值模型的对地支路的（KT-1）/KT算到Ｉ（m）节点
        Y（mm）=Y（mm）-Y（I（n）J（n））+（KT（n）-1）/KT（n）*（1/ZT（n））; 
      elseif J（n）==m   %有变压器时，若支路终止节点为m则变压器等值模型的对地支路的（1-KT）/KT^2算到J（m）节点
         Y（mm）=Y（mm）-Y（I（n）J（n））+（1-KT（n））/（KT（n）^2）*（1/ZT（n））;
       else  Y（mm）=Y（mm）;
       end
   end
end
Y;
G=real（Y）;
%-----------------------求B‘矩阵及其逆矩阵B1
B=imag（Y）;%求Y的虚部节点电纳矩阵
 ph=find（TYPE（:1）==3）;%找出平衡节点编号
 BB=B;%BB矩阵为中间变量
 BB（:ph）=[];%平衡节点编号对应列置空并消去
 BB（ph:）=[];%平衡节点编号对应行置空
 B1=BB;
 B1=inv（B1）;%B1求逆后得B1矩阵
 %-----------------------%求B‘‘及其逆矩阵B2
phpv=find（TYPE（:1）>1）;%找出非PQ节点的编号
 BB=B;     %以下与求B‘相同
 BB（:phpv）=[];
 BB（phpv:）=[];
 B2=BB;    
 B2=inv（B2）;%求得B2矩阵
 
 %-------------计算各节点有功功率不平衡量deltaPi
 k=0;   %k为迭代次数
 kp=0;  %计算P不平衡量deltaPi的收敛标志（0表示不收敛1表示收敛）
 kq=0;  %计算U不平衡量deltaQi的收敛标志（0表示不收敛1表示收敛）
 notph=find（TYPE（:1）<3）;%找出非平衡节点编号
 deltaPi=zeros（x-11）;%deltaPi为（x-1）*1阶矩阵x即为节点数
 pq=find（TYPE（:1）==1）;%找出PQ节点编号
 pqnum=size（B2）;
 pqnum=pqnum（1）;%求PQ节点的个数（因B1矩阵的维数等于PQ节点数）
 deltaQi=zeros（pqnum1）;%deltaQi为pqnum*1阶矩

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件      329484  2019-03-03 19:37  pq分解法\2.jpg
     文件       10094  2019-03-04 20:38  pq分解法\PQchaoliu.m
     文件        2953  2019-03-04 20:39  pq分解法\outputpq.txt
     文件       16255  2019-03-03 11:19  pq分解法\shiyanshuju.xlsx
     目录           0  2019-03-04 20:38  pq分解法\