电力系统稳态分析基于MatLab的潮流计算快速分解法.docx
《电力系统稳态分析基于MatLab的潮流计算快速分解法.docx》由会员分享，可在线阅读，更多相关《电力系统稳态分析基于MatLab的潮流计算快速分解法.docx（22页珍藏版）》请在冰豆网上搜索。
电力系统稳态分析基于MatLab的潮流计算快速分解法
研究生课程设计（论文）
电力系统稳态分析
基于MatLab的潮流计算快速分解法
教学单位自动化学院
姓名贾鹏辉
学号************
年级2013级
专业电力系统及其自动化
指导教师杨伟
职称副教授
2013年12月5日
摘要：
快速分解法是一种定雅克比法，形成系数矩阵、时忽略了支路电阻、对地导纳和理想变压器非标准变化，以及迭代过程中节点电压的不同取值，从而避免每次迭代重新形成雅克比矩阵及其因子表，计算效率大幅提高。
本文采用快速分解法中的XB型算法，并基于MATLAB软件仿真分析具体实例，发现具有较好的收敛性，特别适合在线计算。
关键词：
电力系统潮流计算快速分解法雅可比矩阵
1引言
用牛顿-拉夫逊法计算潮流时，每次迭代都需要重新形成雅克比矩阵并分解因子表。
为避免每次迭代重新形成雅克比矩阵及其因子表，20世纪70年代初Stott提出了快速分解法。
快速分解法采用了一些假设：
a、电力系统有功功率主要受电压相角影响，无功功率主要受电压幅值影响。
b、高压网线路的r<c、系数矩阵、忽略了支路电阻、对地导纳和理想变压器非标准变化。
d、迭代过程中的前的电压幅值用标幺值1代替。
对极坐标型定雅克比法的修正公式
（1）
处理得到简化的修方程式
（2）
在潮流计算中，
（2）式中两个修正方程式交替迭代，Stott把在此基础发展起来的潮流算法称为快速分解法。
快速分解法是目前电力系统进行潮流计算的主要方法。
由于该方法将节点电压的相位和幅值的迭代过程分开进行，使修正方程的系数矩阵维数降低很多，系数矩阵为常数且对称，因而快速分解法具有计算速度快,占用内存少等特点。
2理论基础
以定雅克比牛顿-拉夫逊迭代方程为出发点，具体过程如下：
通过高斯消去法，把牛顿-拉夫逊法的每一次迭代等价地细分为三步计算；对每一步运算作详细分析，证明了在连续的两次牛顿-拉夫逊迭代中，上一次迭代的第三步和下一次迭代的第一步可以合并，从而导出等效的两步式分解算法；论证了该两步式分解算法的系数矩阵与快速分解算法的系数矩阵是一致的。
推导过程并未因用任何解耦的假设。
将极坐标型定雅克比法的修正公式中和用和代替，用代替。
（3）
其中
2.1将原问题分解为,子问题
用高斯消去法消去子块
（4）
记
因，得到
于是定义
解为（5）
在给定的电压幅值和相角初值附近，保持电压相角不变，考虑只有电压幅值的变化时，有功功率的偏差量为
（6）
综合上述结果，如果当前的迭代点为，则第次迭代可以分为三步
（7）
（8）
（9）
2.2简化无功迭代步骤
按式（7）、式（8）和式（9）完成第次迭代后，下面考虑第+1次迭代，有
（10）
利用式（9），上式中的无功功率偏差为
（11）
代入式（10）中，经整理得
（12）
式（12）说明，如果将第次迭代计算出的和用于计算第迭代的无功功率偏差量，则所求得的第次迭代的电压修正量将自动包含第次迭代计算出的。
所以，的计算可以省略，相当于将式（9）与式（10）合并，因此，第次迭代可以两步完成
（13）
（14）
整理式（13）和式（14）得到简化的修正方程式
（15）
因和即为极坐标型定雅克比法的修正公式中和，所以式（15）与式
（2）格式相同。
与XB型快速分解法的修正公式相比，系数矩阵是导纳矩阵的虚部，与相同。
由的定义，有
（16）
对于一般的电网，可能有较复杂的结构。
为了对有直观的认识，假定网络中无接点，则式（16）中各矩阵的维数相等，并且接点导纳矩阵可用节点支路关联矩阵和支路导纳对角矩阵（分别用的和表示电纳和电导）表示。
下面将证明，对于树形电网或所有支路的比值都相同的环形网络，与相等。
如果网络是树状的，其关联矩阵是方阵且非奇异，此时对式（16）有
（17）
式中，为以为支路电纳组成的对角线矩阵；为以为支路电纳建立的节点电纳矩阵。
这说明对树形电网，就是XB型快速分解法中的阵。
对于环形网络，如果电网是均一网，即对任一条支路有，则得
并有
所以
，
故有
（18）
如果电网不是均一网，上述结论不再严格成立。
但和相比，在的零元素处，相应的元素近似等于零；在的非零元素处，相应的元素近似和的非零元素相等。
这可以用下面的例子来说明。
图1四节点电力系统
以图1所示的四节点系统为例，图中给出了支路阻抗。
该例中，和分别为
可见，比更接近于，而用代替即得到XB型快速分解法。
以上推导过程中，只在对有功功率偏差量和无功功率偏差量的计算处做了线性化近似，既没有r<3快速分解法在电力系统安全分析中的应用
某地区电力供应系统由2个发电厂，4个变电站组成，系统模型见图2，其中节点1、2、3为PQ节点，4为PV节点，5为平衡节点。
以100MVA、220KV为基准的节点初始标幺值和网络参数见表1、2。
利用快速分解法，当迭代功率偏差标准为=0.0001时，迭代结束。
图2电力系统模型
节点
1
2
3
4
5
P
-1.60
-2.00
-3.70
5.00
--
Q
-0.80
-1.00
-1.30
1.05
1.05
表1节点负荷值
支路
r
x
b
变比
1-2
0.040
0.250
0.250
-
1-3
0.10
0.350
0．000
-
2-3
0.080
0.300
0.300
-
2-4
0.000
0.015
-
1.050
3-5
0.000
0.030
-
1.050
表2线路参数
4．程序流程图
5．MatLab程序清单
本文采用两种程序来仿真电力系统的潮流，两个程序得到的仿真结果是相同的。
1.程序一
function[Y,U,P,Q,dSij,O,S,Sij,Sji,sumdS]=kuaisufenjie1776_1（）
%输出：
U——节点电压，P--节点有功，Q--节点无功，dSij--支路功率损耗，O为相角；Sij--从节点i流向节点j的功率,S--节点复功率,sumdS--网络总损耗；point为节点矩阵，branch为支路矩阵；
[x]=5;%节点数x
[y]=5;%支路数y
e=0.00001;%误差要求
point=[110-1.6-0.8;
110-2.0-1.0;
110-3.7-1.3;
31.0505.01.05;
21.05001.05];%节点矩阵
branch=[120.040.2500.2500;
130.100.35000;
230.0800.3000.3000;
240.00.01501.05;
350.00.03001.05;];%支路矩阵
TYPE=zeros（x,1）;%TYPE为节点类型矩阵
U=zeros（x,1）;%U为节点电压矩阵
a=zeros（x,1）;%a为节点电压相角矩阵
P=zeros（x,1）;%P为节点有功功率
Q=zeros（x,1）;%Q为节点无功功率
I=zeros（y,1）;%I为起始节点编号矩阵
J=zeros（y,1）;%J为终止节点编号矩阵
Rij=zeros（y,1）;%R为线路电阻
Xij=zeros（y,1）;%X为线路电抗
Zij=Rij+j*Xij;%Yij为线路阻抗
Y=zeros（x）;%Y为n阶节点导纳方阵
G=zeros（x）;%G为n阶节点电导方阵
B=zeros（x）;%B为n阶节点电纳方阵
B0=zeros（y,1）;%B0为n*1阶线路对地电纳值
K=zeros（y,1）;%K为ij支路y*1阶变压器变比,若K=0表示无变压器,K=1则为标准变比,K不等于1为非标准变比
%------------------------------矩阵赋初值：
TYPE=point（:
1）;%将point矩阵的第一列赋给TYPE,以下类似
U=point（:
2）;
O=point（:
3）;
P=point（:
4）;
Q=point（:
5）;
I=branch（:
1）;
J=branch（:
2）;
Rij=branch（:
3）;
Xij=branch（:
4）;
Zij=Rij+j*Xij;
B0=branch（:
5）;
KT=branch（:
6）;
%------------------------------求节点导纳矩阵Y
form=1:
y%求Y中非对角元的元素Yij
Y（I（m）,J（m））=-1/Zij（m）;
Y（J（m）,I（m））=-1/Zij（m）;
end;
form=1:
x%求Y中的Yii
forn=1:
y
if（I（n）==m|J（n）==m）
Y（m,m）=Y（m,m）-Y（I（n）,J（n））+j*B0（n）/2;%Yii为Yij加上线路对地电导的一半乘j
end;
end;
end;
Y
G=real（Y）;
%-----------------------求B'矩阵及其逆矩阵B1
B=imag（Y）;%求Y的虚部,节点电纳矩阵
B11=zeros（x,x）;
form=1:
y
B11（I（m）,J（m））=1/（Xij（m））;
B11（J（m）,I（m））=1/（Xij（m））;
end;
form=1:
x
forn=1:
y
if（I（n）==m|J（n）==m）
B11（m,m）=B11（m,m）-1/（Xij（n））;
end;
end;
end;
p=find（TYPE（:
1）==3）;%找出平衡节点编号
B11（:
p）=[];%平衡节点编号对应行置空
B11（p,:
）=[];%平衡节点编号对应列置空
B11
B1=B11;
B1=inv（B1）;%B1求逆后得B1矩阵
%-----------------------%求B''及其逆矩阵B2
npq=find（TYPE（:
1）>1）;%找出非PQ节点的编号
B22=B;%BB矩阵为中间变量
B22（:
npq）=[];%非PQ节点编号对应行置空
B22（npq,:
）=[];%非PQ节点编号对应行置空
B22
B2=B22;
B2=inv（B2）;%求得B2矩阵
%-------------计算各节点有功功率不平衡量dPi
k=0;%k为迭代次数
ep=0;%计算P不平衡量dPi的收敛标志（0表示不收敛,1表示收敛）
eq=0;%计算U不平衡量dQi的收敛标志（0表示不收敛,1表示收敛）
np=find（TYPE（:
1）<3）;%找出非平衡节点编号
dPi=zeros（x-1,1）;%dPi为x*1阶矩阵,x即为节点数
pq=find（TYPE（:
1）==1）;%找出PQ节点编号
pqn=size（B2）;
pqn=pqn
（1）;%求PQ节点的个数（因B1矩阵的维数等于PQ节点数）
dQi=zeros（pqn,1）;%dQi为pqn*1阶矩阵
while（（~eq|~ep）&（k<100））
k=k+1;
升级会员 