牛顿-拉夫逊迭代法电力网潮流计算方法与程序.docx
《牛顿-拉夫逊迭代法电力网潮流计算方法与程序.docx》由会员分享，可在线阅读，更多相关《牛顿-拉夫逊迭代法电力网潮流计算方法与程序.docx（6页珍藏版）》请在冰豆网上搜索。
6
牛顿-拉夫逊迭代法电力网潮流计算方法与程序
编写佘名寰
牛顿-拉夫逊迭代法在电力网潮流计算中因其收敛性较好获得广泛运用，该算法的难点是需反复计算功率方程中雅可比矩阵各个元素表达式。
本文简叙了牛顿-拉夫逊迭代法潮流计算的基本公式，通过例题介绍了用牛顿-拉夫逊法计算电力网潮流电压的MATLAB程序。
程序采用MATLAB语言的符号矩阵简化了雅可比矩阵系数的计算。
本文可供电力系统电气技术人员和大专院校电力类专业师生参考。
2．牛顿-拉夫逊迭代法潮流计算
2.1极坐标型式节点功率方程式
由电源注入或从负载流出节点的电流统称节点电流，由节点电压和节点电流所求得的节点注入功率为：
Sti=ViIi
Ii=YijVj（I,j=1,2,….n）
Pti=Vi1nVjGijcosθij+Bijsinθij
Qti=Vi1nVjGijsinθij-Bijcosθij（i,j=1,2…..n）（2-01）
式中Pti,QtI——节点注入的有功和无功功率
Vi，Vj——节点i,j电压幅值
θij=θi-θj节点i,j电压的相角差
Gij,Bij节点导纳矩阵的元素，YIJ=GIJ+jBIJ
节点功率平衡关系为：
Pgi-PlI=Pti
Qgi-Qli=Qti
ΔPi=Pgi-PlI-Pti=0
ΔQi=Qgi-Qli-Qti=0（2-02）
Pgi,Qgi ——节点i发电机输入有功和无功功率
Pli,Qli——节点i负荷有功和无功功率
ΔPi，ΔQi--节点i不平衡功率
不平衡功率的微分d（ΔPi）,d（ΔQi）
d∆pi=-（∂pti∂v1∆v1+∂pti∂v2∆v2…+∂pti∂ϑ1∆ϑ1+∂pti∂ϑ2∆ϑ2….）（i=1,2,…,n）
d∆qi=-（∂qti∂v1∆v1+∂q∂v2∆v2…+∂qti∂ϑ1∆ϑ1+∂qti∂ϑ2∆ϑ2….）（2-03）
对于n个节点系统可得如下矩阵形式修正方程式；
∆P1∆P2:
∆Pn∆Q1∆Q2:
∆Qn=∂∆P1∂ϑ1∂∆P2∂ϑ1:
∂∆Pn∂ϑ1∂∆Q1∂ϑ1∂∆Q2∂ϑ1:
∂∆Qn∂ϑ1∂∆P1∂ϑ2∂∆P2∂ϑ2:
∂∆Pn∂ϑ2∂∆Q1∂ϑ2∂∆Q2∂ϑ2:
∂∆Qn∂ϑ2⋯⋯⋯⋯⋯⋯⋯⋯∂∆P1∂ϑn∂∆P2∂ϑn:
∂∆Pn∂ϑn∂∆Q1∂ϑn∂∆Q2∂ϑn:
∂∆Qn∂ϑn∂∆P1∂v1∂∆P2∂v1:
∂∆Pn∂v1∂∆Q1∂v1∂∆Q2∂v1:
∂∆Qn∂v1∂∆P1∂v2∂∆P2∂v2:
∂∆Pn∂v2∂∆Q1∂v2∂∆Q2∂v2:
∂∆Qn∂v2⋯⋯⋯⋯⋯⋯⋯⋯∂∆P1∂vn∂∆P2∂vn:
∂∆Pn∂vn∂∆Q1∂vn∂∆Q2∂vn:
∂∆Qn∂vnΔϑ1∆ϑ2:
∆ϑn∆v1∆v2:
∆vn（2-04）
式中偏微分矩阵为雅可比矩阵。
若系统总共有n个节点，其中r个为PV节点，去掉平衡节点则电压角度θ变量有n-1-个，电压幅值v变量有n-1-r个，方程式共有2n-r-2个。
2.2牛顿-拉夫逊迭代法解算过程；
⑴指定各节点电压的初值；v1,v2,…vn;ϑ1,ϑ2,…ϑn.一般设V=1,θ=0；
⑵将所给定的电压值代入式（2-01），（2-02）计算各节点不平衡功率；
⑶将电压初值代入求雅可比矩阵系数；
⑷解方程式（2-04），求节点电压幅值和相角的修正值∆vi,∆ϑi;
⑸用所求得的电压修正值修正各节点电压：
v
（1）=v（0）+∆v（0）
ϑ
（1）=ϑ（0）+∆ϑ（0）
⑹用新的节点电压重复以上五个步骤计算直到计算精度合符要求。
牛顿-拉夫逊迭代法计算的难点是求雅可比矩阵系数，本文采用MATLAB语言符号矩阵运算jacobian（f,x）命令求解函数f（x）对x的雅可比矩阵。
2.3牛顿-拉夫逊迭代法计算程序
我们通过例题说明牛顿-拉夫逊迭代法计算潮流程序。
【例2.3.1】网络接线如图2-1所示。
支路阻抗分别为Z12=j0.1,Z23=j0.1,Z31=j0.1;三条支路两端的对地电纳皆是j0.01。
给定的注入功率分别为：
S1=-1.7192-j0.7346;S2=0.6661。
节点1是负荷节点即PQ节点；节点2给定的电压大小为V2=1.05,其属于PV节点；节点3是平衡节点，其电压保持定值，大小为V3=1.0。
Z13
Z12
Z23
YC3
YC1
YC1
V3=1.0
3
1
2
YC2
P1=1.7192
Q1=0.7346
P2=0.6661
V2=1.05
Z21=j0.1
Z23=j0.1
Z13=j0.1
YC1=J0.02
YC2=J0.02
YC3=J0.02
图2-1[例2.3.1]网络接线图
2.3.1.1程序参变量说明
Y=G+j*B节点导纳矩阵，G为电导，B为电纳；
P节点有功功率给定值，正为发电机注入，负为负荷功率；
q节点无功功率给定值，正为发电机注入，负为负荷功率；
u节点电压幅值的符号矩阵；
delt节点电压相角的符号矩阵；
pt节点有功功率符号表达式；
qt节点无功功率符号表达式；
pp节点有功功率不平衡值符号表达式；
qq节点无功功率不平衡值符号表达式；
uu节点电压幅值数值矩阵；
dd节点电压相角数值矩阵；
PP节点不平衡功率的数值矩阵
N1网络独立节点总数；
N2网络PV节点总数；
Sm节点功率矩阵；
Smn支路功率矩阵；
J1,J2,J节点不平衡功率雅可比符号矩阵.
JJ节点不平衡功率雅可比数值矩阵
2.3.1.2源程序
源程序如下：
%*****NUE78.mexample2-3-1fig2-1*******************************
%ThefollowingProgramforloadflowcalculationisbasedonMATLAB7.0
clear
u=sym（'[u1,u2,u3]'）;delt=sym（'[d1,d2,d3]'）;
G=zeros（3）;
B=[-19.98,10,10;10,-19.98,10;10,10,-19.98];
Y=G+j*B;
p
（1）=-1.7192;q
（1）=-0.7346;p
（2）=0.6661;
k=0;precision=1;
N1=2;%theN1istheamountofthePQandPVbus
N2=1;%theN2istheamountofthePVbus
form=1:
N1
forn=1:
N1+1
pt（n）=u（m）*u（n）*（G（m,n）*cos（delt（m）-delt（n））+B（m,n）*sin（delt（m）-delt（n）））;
end
pp（m）=p（m）-sum（pt）;
end
form=1:
N1-N2
forn=1:
N1+1
qt（n）=u（m）*u（n）*（G（m,n）*sin（delt（m）-delt（n））-B（m,n）*cos（delt（m）-delt（n）））;
end
qq（m）=q（m）-sum（qt）;
end
J1=jacobian（pp,['d1,d2,u1']）;
J2=jacobian（qq,['d1,d2,u1']）;
J=vertcat（J1,J2）;
uu=[1.0,1.05,1.0];dd=[0,0,0];
whileprecision>0.00001
u1=uu
（1）;u2=uu
（2）;u3=uu（3）;
d1=dd
（1）;d2=dd
（2）;d3=dd（3）;
form=1:
N1
PP（m）=eval（pp（m））;
end
form=1:
N1-N2
PP（N1+m）=eval（qq（m））;
end
JJ=eval（J）;
du=-inv（JJ）*PP';precision=max（abs（du））;
forn=1:
N1
dd（n）=dd（n）+du（n）;
end
forn=1:
N1-N2
uu（n）=uu（n）+du（N1+n）;
end
k=k+1;
end
k-1,dd,uu
%***********************************************
%thefollowingprogramisusedtocalculatetheSmandSmn
forn=1:
N1+1
U（n）=uu（n）*（cos（dd（n））+j*sin（dd（n）））;
end
Um=conj（U'）;
I=Y*Um;
Sm=diag（Um）*conj（I）
form=1:
N1+1
forn=1:
N1+1
Smn（m,n）=U（m）*（conj（U（m））-conj（U（n）））*conj（-Y（m,n））;
end
end
2.3.1.3程序说明
⑴节点导纳矩阵：
Y11=Y22=Y33=0+j（-10-10+0.02）=0-j19.98
Y12=Y23=Y13=0+j10.0
⑵节点电压符号矩阵：
一共三个节点，独立节点数N1=2，PV节点数N2=1
设节点电压幅值符号矩阵为u=sym（'[u1,u2,u3]'）;
节点电压相角符号矩阵为delt=sym（'[d1,d2,d3]'）;
⑶p
（1），q
（1），p
（2）为节点功率给定值，节点2为PV节点，q
（2）待求；
⑷pt（n）,pp（m）相关语句是按公式（2-01）,（2-02）建立ΔPi符号表达式；
ΔP1=pp
（1）=-2149/1250-10*u1*u2*sin（d1-d2）-10*u1*u3*sin（d1-d3）;
ΔP2=pp
（2）=6661/10000+10*u1*u2*sin（d1-d2）-10*u2*u3*sin（d2-d3）
⑸qt（n）,qq（m）相关语句是建立ΔQi符号表达式，只针对PQ节点；
ΔQ1=qq
（1）=-3673/5000-999/50*u1^2+10*u1*u2*cos（d1-d2）+10*u1*u3*cos（d1-d3）
⑹J1，J2，J语句是建立公式（2-03）中的雅可比符号矩阵；
J=
[-10*u1*u2*cos（d1-d2）-10*u1*u3*cos（d1-d3）,10*u1*u2*cos（d1-d2）,-10*u2*sin（d1-d2）-10*u3*sin（d1-d3）]
[10*u1*u2*cos（d1-d2）,-10*u1*u2*cos（d1-d2）-10*u2*u3*cos（d2-d3）,
升级会员 