牛拉法潮流计算程序(附3机9节点结果对比

摘要电力系统潮流计算是研究电力系统稳态运行的一种重要方法,它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态,包括各母线的电压、线路的功率分布以及功率损耗等等。潮流计算主要用于电网规划和静态安全分析，它可为扩建电力网络，以达到规划周期内所需要的输电能力提供依据；也可以对预想事故进行模拟和分析，校核预想事故下的电力系统安全性。本文简单介绍了牛顿-拉夫逊潮流计算的原理、模型与算法，然后用具体的实例，利用MATLAB对牛顿-拉夫逊法的算法进行了验证。关键词：电力系统潮流计算 牛顿-拉夫逊法 MATLAB一、牛拉法的数学模型对一个N节点的电力网路，列写节点电压方程，即I=Y V(1.1)式中，I为节点注入电流列相量，Y为节点导纳矩阵，V为节点电压列相量。由于异地测量的两个电流缺少时间同步信息，以注入功率替换注入电流作为已知量。即(1.2)其中,Yij=Gij+jBij，带入上式，得到有功功率和无功功率方程Pi=Vij=1nVjGijcosij+Bijsinij(1.3)Qi=Vij=1nVjGijsinij-Bijcosij(1.4)大部分情况下，已知PQ，求解V。考虑到电网的功率平衡，至少选择一台发电机来平衡全网有功功率，即至少有一个平衡节点，常选择调频或出线较多的发电机作为平衡节点。具有无功补偿的母线能保持电压幅值恒定，这类节点可作为PV节点。 潮流计算中节点分类总结如下：表 1.1 节点分类PQ节点负荷母线、变电站母线绝大部分是PQ节点已知PQ，求VPV节点有无功储备的发电机母线很小部分是PV节点已知PV（V=1），求Q平衡节点选择一个发电机母线已知V（V=1,=0）,求PQ电压相量作为参考轴=0已知电力系统有m个PQ节点，r个PV节点和1个平衡节点，则可以提取m+r个有功功率方程和m个无功功率方程，从而求解出m+r个和m个V，其余节点的有功和无功可通过式（1.3）、（1.4）求得，这样就完成了潮流计算。二、潮流计算流程用于潮流计算的有功和无功功率方程的阶数一般很高，常选择牛顿-拉夫逊迭代法进行求解高维的非线性方程组。结束maxF(Xk)？XkXk+1=Xk+XkmaxXk0.00001 P=zeros(Nbus,1);%存储所有节点的有功 Q=zeros(Nbus,1);%存储所有节点的无功 for m=1:Nbus %求pv、pq和平衡节点的有功 for n=1:Nbus P(m,1)=P(m,1)+b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3)+B(m,n)*sin(b(m,3)-b(n,3); end end for m=1:Nbus %求pq、pv和平衡节点的无功 for n=1:Nbus Q(m,1)=Q(m,1)+b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3)-B(m,n)*cos(b(m,3)-b(n,3); end end deltp=b(1:Nbus-1,4)-P(1:Nbus-1,1);%pq和pv节点有功差 deltq=b(1:pq,5)-Q(1:pq,1);%pq节点无功差 deltPQ=deltp;deltq; H=zeros(Nbus-1); for m=1:Nbus-1 %求H矩阵 for n=1:Nbus-1 if m=n H(m,n)=-b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3)-B(m,n)*cos(b(m,3)-b(n,3); else H(m,m)=b(m,2)*b(m,2)*B(m,m)+Q(m,1); end end end N=zeros(Nbus-1,pq); for m=1:Nbus-1 %求N矩阵 for n=1:pq if m=n N(m,n)=-b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3)+B(m,n)*sin(b(m,3)-b(n,3); else N(m,m)=-b(m,2)*b(m,2)*G(m,m)-P(m,1); end end end J=zeros(pq,Nbus-1); for m=1:pq %求J矩阵 for n=1:Nbus-1 if m=n J(m,n)=b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3)+B(m,n)*sin(b(m,3)-b(n,3); else J(m,m)=b(m,2)*b(m,2)*G(m,m)-P(m,1); end end end L=zeros(pq,pq); for m=1:pq %求L矩阵 for n=1:pq if m=n L(m,n)=-b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3)-B(m,n)*cos(b(m,3)-b(n,3); else L(m,m)=b(m,2)*b(m,2)*B(m,m)-Q(m,1); end end end Jacobi=H N;J L;%雅克比矩阵 Correction=-JacobideltPQ;%计算电压相角和幅值的修正量 for m=1:Nbus-1 b(m,3)=b(m,3)+Correction(m); end for m=1:pq b(m,2)=b(m,2)+Correction(Nbus-1+m)*b(m,2); endprecision=max(abs(deltPQ); t=t+1 bendb(Nbus,4)=P(Nbus,1);b(Nbus,5)=Q(Nbus,1);for m=pq+1:Nbus-1 b(m,5)=Q(m,1);endb运行结果：四、结果分析与结论运行