牛顿拉夫逊法概要

一、牛顿－拉夫逊法概要首先对一般的牛顿－拉夫逊法作一简单说明。已知一个变量X的函数（4－6）解此方程式时，由适当的近似值X（0）出发，根据（4－7）反复进行计算，当X（n）满足适当的收敛判定条件时就是（4－6）式的根。这样的方法就是所谓的牛顿－拉夫逊法。式（4－7）就是取第n次近似解X（n）在曲线上的点处的切线与X轴的交点作下一次X（n+1）值的方法。参考图4－2（a）。在这一方法中为了能收敛于真解，初值X（0）的选取及函数f（X）必须满足适当的条件，如图4－2（b）所示的那种情况就不能收敛或收敛到别的根上去。这一方法还可以做下面的解释，设第n次迭代得到的解与真值之差，即的误差为时，则（4-8）把在附近对用泰勒级数展开（4-9）上式略去以下的项（4-10）的误差可近似由上式计算出来图4-2（4-11）比较式（4-7）和（4-11），可以看出牛顿-拉夫逊法的修正量和的误差的一次项相等。用同样的方法考虑，给出对n个变量的n个方程式（4－12）对其近似解的修正量，可以解下面的方程式来确定（4－13）式（4－13）的右边的矩阵的等都是对于的值。这一矩阵称为雅可比（Jacobi）矩阵。按上述得到的修正量后，得到如下关系：这比进一步接近于真值。这一步骤在收敛到希望的值以前重复进行。一般要反复计算到满足时为止。ε为预先规定的小正数，此处是第n次迭代Xi的近似值。一、牛顿－拉夫逊法潮流计算把牛顿法用于潮流计算，要求将潮流方程改写成形如方程式（4－12）所示的形式。为此，首先应将潮流方程（4－5）的变形式的右端展开，并且分开实部和虚部。采用直角坐标时，节点电压可表示为：节点导纳矩阵元素则表示为：将上述表示式代入的右端，展开并分出实部和虚部，便得：（4－14）按照上节的分类，PQ节点的有功功率和无功功率给定的，第I个节点的给这功率设为Pis和Qis。假定系统中的第1，2，………m号节点为PQ节点，对其中每一个节点可列（i=1，2，…………，m）（4－15）PV节点的有功功率和节点电压幅值是给定的。假定系统中的第m+1，m+2，………n-1号节点为PV节点，则对其中每一节点可以列写方程：（4－16）第n号节点为平衡节点，其电压是给定的，故不参加迭代。式（4－15）和（4－16）总共包含了2（n-1）个方程，待求的变量有也是2（n-1）个。我们还可以看到，方程式（4－15）和（4－16）已经具备方程组（4－12）的形式：（4－16）'式中上述方程中雅可比矩阵的各元素，可以对（4－15）和（4－16）式求偏导数获得。当时，对角元素是（4-17）当时，矩阵中非对角元素是（4－18）由以上表达式不难看出，雅可比矩阵有以下特点：（1）雅可比矩阵中的诸元素都是节点电压的函数，因此在迭代过程中，它们将随着各节点电压的变化而不断地改变；（2）矩阵是不对称的；（3）由式（4－18）可以看出，当导纳矩阵中的非对角元素Yij为零时，雅可比矩阵中相对应的元素也是零，即矩阵是非常稀疏的。因此，修正方程的求解同样可以应用稀疏矩阵的求解技巧。正是由于这一点才使N－R法获得广泛的应用。三、牛顿法的框图及求解过程程序框图如下：用牛顿法计算潮流时，有以下的步骤：（1）给这各节点电压初始值（2）将以上电压初始值代入式（4－15）和（4－16），求出修正方程式的常数项向量；（3）将电压初始值再代入式（4－17）和（4－18），求出修正方程式中系数矩阵（雅可比矩阵）的各元素；（4）解修正方程式（4－16），求出修正量；（5）修正各节点电压（6）将再代入（4－15），（4－16）式求；（7）校验是否收敛，即（8）如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果。如果不收敛，转回第（2）步进行下一次迭代计算，直到收敛为止。四、实例及程序清单例4－1试用牛顿－拉夫逊法计算图4－4所示电力系统的潮流分布。解㈠需要输入的数据⑴NN－节点数、NL－线路数、ISB－平衡母线节点号（此例为节点1）、PR－误差精度,如果要输入则输入eps即可.⑵请输入由支路参数形成的矩阵B1矩阵B1的每行是由下列参数构成的：○1某支路的首端号P。○2末端号Q;且PQ。○3支路的阻抗（R+jX）。○4支路的对地容抗。○5支路的变比K。○6折算到哪一侧的标志（如果支路的首端P处于高压侧则请输入1，否则请输入0）。⑶请输入各节点参数形成的矩阵B2矩阵B2的每行是由下列参数构成的：○1节点所接发电机的功率SG。○2节点负荷的功率SL。○3节点电压的初始值。○4PV节点电压V的给定值。○5节点所接的无功补偿设备的容量。○6节点分类标号igl。1－－平衡节点igl=2－－PQ节点3－－PV节点（二）先形成节点导纳矩阵，同例2－3（三）根据式（4－15）和（4－16）求出修正方程式的常数项向量。（四）据式（4－17）和（4－18）求出雅可比矩阵各元素值，即可得到第一次迭代时的修正方程式。（五）方程式，求。（六）各节点电压，即得出第一次迭代后各节点的电压值。（七）计算步骤迭代下去，当收敛精度取（即PR=0.0001）时，需要进行四次迭代。求出了各节点电压后，即可求各支路的潮流分布。见图4－5程序清单及打印结果：%本程序的功能是用牛顿－拉夫逊法进行潮流计算n=input('请输入节点数:n=');nl=input('请输入支路数:nl=');isb=input('请输入平衡母线节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入由支路参数形成的矩阵:B1=');B2=input('请输入各节点参数形成的矩阵:B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n);S1=zeros(nl);fori=1:nlifB1(i,6)==0p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endY(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)+B1(i,4)./2;Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;end%求导纳矩阵disp('导纳矩阵Y=');disp(Y);G=real(Y);B=imag(Y);fori=1:ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfori=1:nS(i)=B2(i,1)-B2(i,2);B(i,i)=B(i,i)+B2(i,5);endP=real(S);Q=imag(S);ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;whileIT2~=0IT2=0;a=a+1;fori=1:nifi~=isbC(i)=0;D(i)=0;forj1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);endP1=C(i)*e(i)+f(i)*D(i);Q1=f(i)*C(i)-D(i)*e(i);%求'P,Q'V2=e(i)^2+f(i)^2;ifB2(i,6)~=3DP=P(i)-P1;DQ=Q(i)-Q1;forj1=1:nifj1~=isb&j1~=iX1=-G(i,j1)*e(i)-B(i,j1)*f(i);X2=B(i,j1)*e(i)-G(i,j1)*f(i);X3=X2;X4=-X1;p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;elseifj1==i&j1~=isbX1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;endendelseDP=P(i)-P1;DV=V(i)^2-V2;forj1=1:nifj1~=isb&j1~=iX1=-G(i,j1)*e(i)-B(i,j1)*f(i);X2=B(i,j1)*e(i)-G(i,j1)*f(i);X5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;elseifj1==i&j1~=isbX1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);X5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;endendendendend%求雅可比矩阵fork=3:N0k1=k+1;N1=N;fork2=k1:N1J(k,k2)=J(k,k2)./J(k,k);endJ(k,k)=1;ifk~=3k4=k-1;fork3=3:k4fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endifk==N0,break;endfork3=k1:N0fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endelsefork3=k1:N0fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endendendfork=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N);k1=k+1;f(L)=f(L)-J(k1,N);endfork=3:N0DET=abs(J(k,N));ifDET=prIT2=IT2+1;endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消云法解w=-J*Vdisp('迭代次数');disp(ICT1);disp('没有达到精度要求的个数');disp(ICT2);fork=1:nV(k)=sqrt(e(k)^2+f(k)^2);O(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;enddisp('各节点的实部电压标么值E为(节点号从小到大排):');disp(E);disp('各节点的电压大小V为(节点号从小到大排):');disp(V);disp('各节点的电压相角O为(节点号从小到大排):');disp(O);forp=1:nC(p)=0;forq=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp('各节点的功率S为(节点号从小到大排):');disp(S);disp('各条支路的首端功率Si为(顺序同您输入B1时一样):');fori=1:nlifB1(i,6)==0p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));disp(Si(p,q)