复杂网络NR法潮流分析与计算的设计.docx
《复杂网络NR法潮流分析与计算的设计.docx》由会员分享,可在线阅读,更多相关《复杂网络NR法潮流分析与计算的设计.docx(45页珍藏版)》请在冰豆网上搜索。
复杂网络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+