复杂网络NR法潮流分析与计算的设计.docx
《复杂网络NR法潮流分析与计算的设计.docx》由会员分享，可在线阅读，更多相关《复杂网络NR法潮流分析与计算的设计.docx（37页珍藏版）》请在冰豆网上搜索。
复杂网络NR法潮流分析与计算的设计
电气工程及其自动化专业课程设计
复杂网络N-R法潮流分析与计算的设计
课程设计报告撰写内容
一、设计要求（宋体，小四号字，加黑）
用matlab编程，N_R法计算潮流分布
具体要求为：
（1）给出程序，并给出注释
（2）输出迭代次数，各节点电压，各支路电流
（3）在图中标明功率流向
节点数据如下表所示（标幺值）
1
2
3
4
5
6
P
3
1.8
0.6
3.5
5
Q
1
0.5
0.8
1.3
V
1
1.05
支路及变压器数据
线路
T1
T2
L2
L3
L4
L5
阻抗
j0.04
j0.02
0.06+j0.025
0.01+j0.2
0.06+j0.5
0.05+j0.3
导纳/2
j0.25
j0.25
j0.25
j0.25
变比
1.05:
1
1.05:
1
精度要求：
0.0001
二、设计方案（要求给出详细的设计思路及其必要的论证）
（1.）潮流计算的方法
（1）高斯雅克比迭代法
（2）高斯-塞得尔法（对初值要求底，迭代次数多）
（3）牛顿-拉夫逊法（使用广泛）
（4）PQ快速分解法（提升运算速度）
目前广泛应用的潮流计算方法都是基于节点电压法的，以节点导纳矩阵Y
作为电力网络的数学模型。
节点电压Ui和节点注入电流Ii由节点电压方程
YV=I
（1）
根据S=VI﹡（I﹡为I的共轭）可得非线性的节点方程
YV=I=（S/V）﹡
（2）
在实际的电力系统中，已知的运行条件不是节点的注入电流，而是负
荷和发电机的功率，而且这些功率一般不随节点电压的变化而变化。
由于
各节点注入功率与注入电流的关系为Si＝Pi＋jQi=ViIi﹡，因此可将式
（2）
改写为
Ii=Si/Vi=Pi+jQi/Vi（i=1,2,3⋯n）（3）
式中，Pi和Qi分别为节点i向网络注入的有功功率和无功功率，当
i为发电机节点时Pi﹥0；当i为负荷节点时Pi﹤0；当i为无源节点Pi＝
0，Qi＝0；Vi和Ii分别为节点电压相量Vi和节点注入电流相量Ii的共
轭。
式（3）亦即潮流计算的基本方程式,它可以在直角坐标也可以在极坐
标上建立2n个实数形式功率方程式。
发电机Pi、Qi为正，负荷Pi、Qi为
负。
展开YV=I为
Ii=ΣYijVj=YiiVi+ΣYijVi（i=123⋯n）（4）
将式（4）代入式（3），得n维的非线性复数的电压方程组
潮流计算的基本方程为
（Pi-jQi）/Vi=YiiVi+ΣYijVi（i=1,2,⋯n）（5）
（2.）变量的分类
假设系统中有n个节点，构成n个复数方程，2n个实数方程，变量总数为6n个。
a）不可控变量（2n个）：
负荷消耗的有功功率LiP和无功功率LiQ.由于该类变
量无法控制，取决于用户，而且出现事先没有预计的变动，使系统偏离原始运行
状态，因此又称为不可控变量或扰动变量。
b）控制变量（2n个）:
发电机发出的有功功率GiP和无功功率GiQ，因为该类
变量可控。
也称独立变量。
c）状态变量（2n个）：
母线电压或节点电压的幅值大小iV与相角大小iδ，又
称依从变量或因变量。
并且iV受GiP控制，iδ受GiQ控制。
其中2n个扰动变量是给定的，2n个控制变量和2n个状态变量中给定两个，求
另外两个。
（3.）变量的约束条件
a）扰动变量没有约束条件。
b）控制变量约束条件：
为满足发电机的技术经济特性指标。
c）状态变量的iV的约束条件：
保证良好的电能质量。
d）状态变量的iδ的约束条件：
保证系统的稳定运行。
（4.）系统节点的分类，根据给定的控制变量和状态变量进行分类如下：
（1）PQ节点（即负荷节点）：
给定Pi、Qi，求Vi和iδ（iie,f）。
通常变电所都是这一类型的节点，由于
没有发电设备，因而发电功率为零电力系统中的绝大多数节点属于这一节点。
其
包含变电站节点（即联络节点或浮游节点）。
（2）PV节点（即调节节点、电压控制节点）：
给定Pi和Vi，求Qi和iδ（iie,f）。
这类节点必须有足够的可调无功容量，
用以维持给定的电压幅值。
一般时选择有一顶武功储备的发电厂和具有可调无功
电源设备的变电所作为PV节点。
在电力系统中，这类节点数很少。
（3）平衡节点（即松弛节点、参考节点、基准节点）：
给定Vi和iδ（iδ=0），求Pi和Qi。
（只有一个）有功功率不能给定，这个节
点承担了系统的有功功率平衡。
同时其电压幅值也是给定的，相位为零。
（5.）P-Q分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思
想是：
把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功功率
误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，
把有功功率和无功功率迭代分开来进行。
牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系
统时，修正方程式展开为：
ΔP=HΔΘ+NΔV/V
ΔQ=JΔΘ+LΔV/V
以上方程式是从数学上推倒出来的，并没有考虑电力系统这个具体对象的特
点。
电力系统中有功功率主要与各节点电压向量的角度有关，无功功率则主要受
各节点电压幅值的影响。
大量运算经验也告诉我们，矩阵N及J中各元素的数
值相对是很小的，因此对牛顿法的第一步简化就是把有功功率和无功功率分开
来进行迭代，即将式（4）化简为：
ΔP=HΔΘ
ΔQ=LΔV/V（5）
这样，由于我们把2n阶的线性方程组变成了二个n阶的线性方程组，对牛
顿法的第二个化简，也是比较关键的一个化简，即把式（5）中的系数矩阵简化为
在迭代过程中不变的对称矩阵。
众所周知，一般线路两端电压的相角差是不大
的（通常不超过10~20度），因此可以认为：
（6）
此外，与系统各节点无功功率相应的导纳LiB必定远远小于该节点自导纳的
虚部，即：
因此，
（7）
考虑到以上关系后，式（5）中系数矩阵中的元素表达式可以化简为：
（8）
这样，式（5）中系数矩阵可以表示为：
（9）
进一步可以把它们表示为以下矩阵的乘积：
（10）
将它代入（5）中，并利用乘法结合率，我们可以把修正方程式变为：
将以上两式的左右两侧用以下矩阵左乘
就可得到
以上两式就是P-Q分解法达到修正方程式，其中系数矩阵只不过是系统导纳
矩阵的虚部，因而是对称矩阵，而且在迭代过程中维持不变。
它们与功率误差
方程式
构成了P-Q分解法迭代过程中基本计算公式，其迭代步骤大致是：
（1）给定各节点电压向量的电压初值Vi（0）,θi（0）;
（2）根据（12）计算各节点有功功率误差ΔPi，并求出;ΔPi/Vi
（3）解修正方程式（11），并进而计算各节点电压向量角度的修正量iΔθ
（4）修正各节点电压向量角度θi；
（5）根据式（16）计算各节点无功功率误差iΔQ，并求出/;iiΔQV
（6）解修正方程式（11），求出各节点电压幅值的修正量iΔV
（7）修正各节点电压幅值iV（i）（i1）（i1）
iiiV=V−−ΔV−（18）
（8）返回
（2）进行迭代，直到各节点功率误差及电压误差都满足收敛条件。
P-Q分解法与牛顿法潮流程序的主要差别表现在它们的修正方程式上。
P-Q
分解法通过对电力系统具体特点的分析，对牛顿法修正方程式的雅可比矩阵进
行了有效的简化和改进，有以下三个特点：
（1）在提高计算速度和减少内存方面的作用是明显的，不再叙述。
（2）使我们得到以下好处。
首先，因为修正方程式的系数矩阵就是导
纳矩阵的虚部，因此在迭代过程中不必象牛顿法那样进行形成雅可比矩阵的
计算，这样不仅是仅减少了运算量，而且也大大简化了程序。
其次，由于系数
矩阵在迭代过程中维持不变，因此在求解修正方程式时，可以迅速求得修正量，
从而显著提高了迭代速度。
（3）可以使我们减少形成因子表时的运算量，而且由于对称矩阵三角分解
后，其上三角矩阵和下三角矩阵有非常简单的关系，所以在计算机中可以只储
存上三角矩阵或下三角矩阵，从而也进一步节约了内存。
三、设计内容
%本程序的功能是用牛顿——拉夫逊法进行潮流计算
%B1矩阵：
1、支路首端号；2、末端号；3、支路阻抗；4、线路对地电纳（或变压器导纳）；
%5、支路的变比；6、支路首端处于K侧为1，1侧为0；
%7、线路/变压器标识（0/1）变压器参数当支路首端处于K侧标识为1时归算至末端侧，0归算至首端侧
%B2矩阵：
1、该节点发电机功率；2、该节点负荷功率；
%3、PQ节点电压初始值，或平衡节点及PV节点电压的给定值
%4、节点所接无功补偿并联电容（感）的电纳
%5、节点分类标号：
1为平衡节点（应为1号节点）；2为PQ节点；3为PV节点；
clear;
isb=1;%input（'请输入平衡母线节点号：
isb='）;
pr=0.0001;%input（'请输入误差精度：
pr='）;
%---------------------------------------------------
n=6;%input（'请输入节点数：
n='）;
nl=6;%input（'请输入支路数：
nl='）;
B1=[120+0.04i01.0511;
230.06+0.025i0+0.5i100;
250.01+0.2i0+0.5i100;
340.06+0.50i0+0.5i100;
450.05+0.3i0100;
650+0.02i01.0511]
B2=[00101;
03+1i1.0002;
01.8+0.50i1.0002;
00.6+0.8i1.0002;
03.5+1.3i1.0002;
0-5+0i1.0503]
%input（'请输入各节点参数形成的矩阵：
B2='）;
%X=[10;20;30;40;50;60]
%-------------------------------------------------------------
%n=4;%input（'请输入节点数：
n='）;nl=4;%input（'请输入支路数：
nl='）;
%B1=[124+16i0100;134+16i0100;232+8i0100;241.49+48.02i011/11001]%input（'请输入由支路参数形成的矩阵：
B1='）;
%B2=[0011501;0011002;020+4i11002;010+6i1002]%input（'请输入各节点参数形成的矩阵：
B2='）;
%-------------------------------------------------------------
Y=zeros（n）;e=zeros（1,n）;f=zeros（1,n）;V=zeros（1,n）;sida=zeros（1,n）;S1=zeros（nl）;
%%%-----------求导纳矩阵------------------------
%fori=1:
n
%ifX（i,2）~=0;
%p=X（i,1）;
%Y（p,p）=1/X（i,2）;
%end
%end
fori=1:
nl%从1到n1（总支路数）
ifB1（i,7）==1%-----------如果是变压器支路--------
ifB1（i,6）==0%左节点（首端）处于1侧
p=B1（i,1）;q=B1（i,2）;
else%左节点（首端）处于K侧
p=B1（i,2）;q=B1（i,1）;
end
Y（p,q）=Y（p,q）-1./（B1（i,3）*B1（i,5））;%非对角元
Y（q,p）=Y（p,q）;%非对角元
Y（q,q）=Y（q,q）+1./（B1（i,3）*B1（i,5）^2）;%对角元K侧
Y（p,p）=Y（p,p）+1./B1（i,3）+B1（i,4）;%对角元1侧+励磁导纳
else%------------否则为线路支路--------------------
p=B1（i,1）;q=B1（i,2）;
Y（p,q）=Y（p,q）-1./B1（i,3）;%非对角元
Y（q,p）=Y（p,q）;%非对角元
Y（q,q）=Y（q,q）+1./B1（i,3）+B1（i,4）./2.0000;%对角元j侧+线路电纳的一半
Y（p,p）=Y（p,p）+1./B1（i,3）+B1（i,4）./2.0000;%对角元i侧+线路电纳的一半
end
end
disp（'导纳矩阵Y='）;disp（Y）;
%-----------给定各节点初始电压及给定各节点注入功率--------------------------
G=real（Y）;B=imag（Y）;%分解出导纳阵的实部和虚部
fori=1:
n%给定各节点初始电压的实部和虚部
e（i）=real（B2（i,3））;
f（i）=imag（B2（i,3））;
V（i）=abs（B2（i,3））;%PV、平衡节点及PQ节点电压模值
end
fori=1:
n%给定各节点注入功率
S（i）=B2（i,1）-B2（i,2）;%i节点注入功率SG-SL
B（i,i）=B（i,i）+B2（i,4）;%i节点无功补偿量（电纳值）
end
%==================用牛顿-拉夫逊法迭代求解非线性代数方程（功率方程）=======================
P=real（S）;Q=imag（S）;%分解出各节点注入的有功和无功功率
ICT1=0;IT2=1;N0=2*n;N1=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2
whileIT2~=0%N0=2*n雅可比矩阵的阶数；N1=N0+1扩展列
IT2=0;a=a+1;
JZ=['Jacobi矩阵第（',num2str（a）,'）次消去运算'];JZ1=['Jacobi矩阵第（',num2str（a）,'）次回代运算'];JZ0=['功率方程第（',num2str（a）,'）次差值：
'];
%----------------求取各个节点的功率及功率偏差及PV节点的电压偏差--------------------
fori=1:
n%n个节点2n行（每节点两个方程P和Q或U）
p=2*i-1;m=p+1;C（i）=0;D（i）=0;
forj1=1:
n%第i行共n列（n个节点间互导纳及节点电压相乘即电流）
C（i）=C（i）+G（i,j1）*e（j1）-B（i,j1）*f（j1）;%Σ（Gij*ej-Bij*fj）
D（i）=D（i）+G（i,j1）*f（j1）+B（i,j1）*e（j1）;%Σ（Gij*fj+Bij*ej）
end
%求i节点有功和无功功率P',Q'的计算值
P1=C（i）*e（i）+f（i）*D（i）;%节点功率P计算eiΣ（Gij*ej-Bij*fj）+fiΣ（Gij*fj+Bij*ej）
Q1=C（i）*f（i）-e（i）*D（i）;%节点功率Q计算fiΣ（Gij*ej-Bij*fj）-eiΣ（Gij*fj+Bij*ej）
V2=e（i）^2+f（i）^2;%电压模平方
%===求取功率差及PV节点电压模平方差=========
ifi~=isb%非平衡节点（PQ或PV节点）
ifB2（i,5）~=3%非PV节点（只能是PQ节点）
J（m,N1）=P（i）-P1;%PQ节点有功功率差J（m,N1）扩展列△P
J（p,N1）=Q（i）-Q1;%PQ节点无功功率差J（p,N1）扩展列△Q
else%PV节点==================
J（m,N1）=P（i）-P1;%PV节点有功功率差J（m,N1）扩展列△P
J（p,N1）=V（i）^2-V2;%PV节点电压模平方差J（p,N1）扩展列△U
end
end%（ifi~=isb）非平衡节点（PQ或PV节点）
end%（fori=1:
n）n个节点2n行（每节点两个方程P和Q或U）
form=1:
N0
JJN1（m）=J（m,N1）;
end
disp（JZ0）;disp（JJN1）;
%-------------判断功率偏差量及PV节点的电压偏差量是否满足要求-----------------
fork=3:
N0%除去平衡节点1、2号以外的所有节点
DET=abs（J（k,N1））;
ifDET>=pr;%PQ节点的功率偏差量及PV节点的电压偏差量是否满足要求
IT2=IT2+1;%不满足要求的节点数加1
end
end
ICT2（a）=IT2;%不满足要求的节点数；a为迭代次数
ICT1=ICT1+1;%迭代次数
ifICT2（a）==0;%当前不满足要求的节点数为零
break%退出迭代运算
end
%--------------------以上为求取各个节点的功率及功率偏差及PV节点的电压偏差-------------
%=================求取Jacobi矩阵形成修正方程===================
fori=2:
n%n个节点2n行（每节点两个方程P和Q或U）
ifi~=isb%非平衡节点（PQ或PV节点）
ifB2（i,5）~=3%下面是针对PQ节点来求取Jacobi矩阵的元素===========
C（i）=0;D（i）=0;
forj1=1:
n%第i行共n列（n个节点间互导纳及节点电压相乘即电流）
C（i）=C（i）+G（i,j1）*e（j1）-B（i,j1）*f（j1）;%Σ（Gij*ej-Bij*fj）
D（i）=D（i）+G（i,j1）*f（j1）+B（i,j1）*e（j1）;%Σ（Gij*fj+Bij*ej）
end
forj1=2:
n%第i行共n列（2n个Jacobi矩阵元素dP/de及dP/df或dQ/de及dQ/df）
ifj1~=isb&j1~=i%非平衡节点&非对角元
X1=-G（i,j1）*e（i）-B（i,j1）*f（i）;%X1=dP/de=-dQ/df=-X4
X2=B（i,j1）*e（i）-G（i,j1）*f（i）;%X2=dP/df=dQ/de=X3
X3=X2;%X2=dp/dfX3=dQ/de
X4=-X1;%X1=dP/deX4=dQ/df
p=2*i-1;q=2*j1-1;
J（p,q）=X3;m=p+1;%X3=dQ/deJ（p,N）=DQ节点无功功率差J（p,N）=DQ;
J（m,q）=X1;q=q+1;%X1=dP/deJ（m,N）=DP节点有功功率差J（m,N）=DP;
J（p,q）=X4;J（m,q）=X2;%X4=dQ/dfX2=dp/df
elseifj1==i&j1~=isb%非平衡节点&对角元
X1=-C（i）-G（i,i）*e（i）-B（i,i）*f（i）;%dP/de
X2=-D（i）+B（i,i）*e（i）-G（i,i）*f（i）;%dP/df
X3=D（i）+B（i,i）*e（i）-G（i,i）*f（i）;%dQ/de
X4=-C（i）+G（i,i）*e（i）+B（i,i）*f（i）;%dQ/df
p=2*i-1;q=2*j1-1;J（p,q）=X3;%扩展列△QJ（p,N）=DQ;
m=p+1;
J（m,q）=X1;q=q+1;J（p,q）=X4;%扩展列△PJ（m,N）=DP;
J（m,q）=X2;
end
end
else%ifB2（i,5）~=3%否则（即为PV节点）
%===============下面是针对PV节点来求取Jacobi矩阵的元素===========
forj1=1:
n
ifj1~=isb&j1~=i%非平衡节点&非对角元
X1=-G（i,j1）*e（i）-B（i,j1）*f（i）;%dP/de
X2=B（i,j1）*e（i）-G（i,j1）*f（i）;%dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J（p,q）=X5;%PV节点电压误差J（p,N）=DV;
m=p+1;
J（m,q）=X1;q=q+1;J（p,q）=X6;%PV节点有功误差J（m,N）=DP;
J（m,q）=X2;
elseifj1==i&j1~=isb%非平衡节点&对角元
X1=-C（i）-G（i,i）*e（i）-B（i,i）*f（i）;%dP/de
X2=-D（i）+B（i,i）*e（i）-G（i,i）*f（i）;%dP/df
X5=-2*e（i）;
X6=-2*f（i）;
p=2*i-1;q=2*j1-1;J（p,q）=X5;%PV节点电压误差J（p,N）=DV;
m=p+1;
J（m,q）=X1;q=q+1;J（p,q）=X6;%PV节点有功误差J（m,N）=DP;
J（m,q）=X2;
end
end
end%（ifB2（i,5）~=3else）
end%（ifi~=isb）
end%（fori=1:
n）n个节点2n行（每节点两个方程P和Q或U）
JZ0=['形成的第（',num2str（a）,'）次Jacobi矩阵：
'];
disp（JZ0）;disp（J）;
%===============================以上为形成完整的Jacobi矩阵================================
%====下面用高斯消去法对由Jacobi矩阵形成的修正方程进行求解（按列消去、回代）==========
fork=3:
N0%N0=2*n（从第三行开始，第一、二行是平衡节点）
fork1=k+1:
N1%从k+1列的Jacobi元素到扩展列的△P、△Q或△U
J（k,k1）=J（k,k1）./J（k,k）;%用K行
升级会员 