基于牛顿-拉夫逊法进行MATLAB电力系统潮流分析

1.程序包括三机九节点数据，潮流分析程序，程序可直接运行。
2.求解方法包括极坐标和直角坐标系下的牛顿-拉夫逊求解方法，以及考虑雅克比矩阵的对称稀疏性时的牛顿-拉夫逊方法。
3.只需要更换输入数即可用来求解其他电力系统潮流

% IEEE-9标准试验系统潮流求解，直角坐标系求解
% 1号节点为平衡节点
% 基准功率100MW，基准电压230KV

clear

% time=zeros（1001）;
% for t=1:100
% tic

%% 读Data1中数据
load ieee9_node.txt
Bus=ieee9_node（:1）;                    %节点号
Vtype=ieee9_node（:3）;                  %节点类型
Pload=ieee9_node（:4）;                  %负载有功功率
Qload=ieee9_node（:5）;                  %负载无功功率
Pgen=ieee9_node（:6）;                   %发电机发出的有用功率
Qgen=ieee9_node（:7）;                   %发电机发出的无功功率
Vset=ieee9_node（:9）;                   %电压设定值
Qsh=ieee9_node（:12）;                   %并联电容器的电纳（标幺值）

swing=find（Vtype==3）;                   %平衡节点编号
PV=find（Vtype==2）;                      %PV节点编号
PQ=find（Vtype==1）;                      %PQ节点编号
PVPQ=find（Vtype<3）;                     %PV/PQ节点编号
%% 读取Data2中的数据
load ieee9_branch.txt
II=ieee9_branch（:1）;                      %支路首端节点
JJ=ieee9_branch（:2）;                      %支路末端节点
Ltype=ieee9_branch（:3）;                   %支路类型
R=ieee9_branch（:4）;                       %支路电阻
X=ieee9_branch（:5）;                       %支路电抗
B_=ieee9_branch（:6）;                      %1/2线路电纳
K=ieee9_branch（:7）;                       %变压器非标准变比
%% 节点导纳矩阵
n=length（Bus）;                      %节点数量
y1=zeros（n）;                        %支路/变压器导纳
y3=zeros（n）;                        %支路对地导纳

y1（（JJ-1）*n+II）=1./（（R+1i*X）.*K）;
y1（（II-1）*n+JJ）=y1（（JJ-1）*n+II）;
y3（（JJ-1）*n+II）=（1-K）./（（R+1i*X）.*K.^2）+1i*B_;
y3（（II-1）*n+JJ）=（K-1）./（（R+1i*X）.*K）+1i*B_;
Y=diag（sum（y12）+sum（y32）+1i*Qsh）;       %节点导纳矩阵
Y=Y-y1;
G=real（Y）;                                          %节点电导
B=imag（Y）;                                          %节点电纳
%% 设电压初始值
U=Vset;
e=real（U）;
f=imag（U）;
%% 牛顿-拉夫逊法迭代求解
Ps=Pgen-Pload;                                      %节点注入有功功率
Qs=Qgen-Qload;                                      %节点注入无功功率
D=ones（2*（n-1）1）;                                  %节点功率不平衡量
N=0;                                                %迭代次数
J=zeros（2*（n-1））;                                   %雅克比矩阵
tolerance=1e-5;                                     %误差精度
while max（abs（D））>tolerance
    %计算节点功率和节点电压平方的不平衡量（不包括平衡节点）
    
    D（2*PVPQ-3）=Ps（PVPQ）-e（PVPQ）.*（G（PVPQ:）*e-B（PVPQ:）*f）-f（PVPQ）.*（G（PVPQ:）*f+B（PVPQ:）*e）;     %有功不平衡
    D（2*PQ-2）=Qs（PQ）-f（PQ）.*（G（PQ:）*e-B（PQ:）*f）+e（PQ）.*（G（PQ:）*f+B（PQ:）*e）;                     %PQ节点无功不平衡
    D（2*PV-2）= Vset（PV）.^2-（e（PV）.^2+f（PV）.^2）;                                                     %PV节点电压不平衡
    %计算雅克比矩阵
    for a=1:n-1
        j=PVPQ（a）;
        for b=1:n-1
            k=PVPQ（b）;
            if j~=k
                J（2*j-32*k-3）=-B（jk）*e（j）+G（jk）*f（j）;
                J（2*j-32*k-2）=G（jk）*e（j）+B（jk）*f（j）;
                if Vtype（j）==1                     %PQ节点
                    J（2*j-22*k-3）=-J（2*j-32*k-2）;
                    J（2*j-22*k-2）=J（2*j-32*k-3）;
                end
            else
                J（2*j-32*k-3）=-B（jk）*e（j）+G（jk）*f（j）+G（j:）*f+B（j:）*e;
                

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        4982  2018-11-15 11:00  ieee9_1.m
     文件        5883  2018-11-15 10:50  ieee9_12.m
     文件        7040  2018-11-15 10:50  ieee9_2.m
     文件        9298  2018-11-15 10:50  ieee9_22.m
     文件        1932  2018-11-11 21:36  ieee9_branch.txt
     文件        2080  2018-11-14 21:25  ieee9_node.txt