Matlab牛拉法潮流计算程序
潮流计算(matlab)实例计算
潮流例题：根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择Matlab进行设计）。
2.在给定的电力网络上画出等值电路图。
3.运用计算机进行潮流计算。
4.编写设计说明书。
一、设计原理1．牛顿-拉夫逊原理牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。
牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。
电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。
为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。
牛顿—拉夫逊迭代法的一般步骤：（1）形成各节点导纳矩阵Y。
（2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。
（3）计算各个节点的功率不平衡量。
（4）根据收敛条件判断是否满足，若不满足则向下进行。
（5）计算雅可比矩阵中的各元素。
（6）修正方程式个节点电压（7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。
（8）计算平衡节点输出功率和各线路功率2．网络节点的优化１）静态地按最少出线支路数编号这种方法由称为静态优化法。
matlab3节点牛顿拉夫逊潮流计算
1. Matlab3节点牛顿拉夫逊潮流计算简介Matlab是一种高度灵活的编程语言和数学工具，被广泛应用于科学计算和工程领域。
3节点牛顿拉夫逊潮流计算是一种电力系统分析方法，用于计算电力系统中各节点的电压和相角。
在本文中，我们将介绍如何使用Matlab进行3节点牛顿拉夫逊潮流计算，以及该方法的原理和应用。
2. 3节点牛顿拉夫逊潮流计算原理3节点牛顿拉夫逊潮流计算是一种基于潮流方程的迭代算法。
它通过不断迭代求解节点电压和相角，以达到系统在给定负荷下的稳态。
其核心原理是利用牛顿拉夫逊法迭代求解潮流方程，即功率平衡方程和节点电压方程，直至收敛得到结果。
3. Matlab在3节点牛顿拉夫逊潮流计算中的应用Matlab提供了丰富的数学工具和函数库，使其成为进行电力系统分析的理想工具。
在3节点牛顿拉夫逊潮流计算中，我们可以利用Matlab编写相应的算法和程序，对实际电力系统进行分析和计算。
通过Matlab的矩阵运算和迭代算法，可以高效地求解潮流方程，得到系统各节点的电压和相角。
4. 3节点牛顿拉夫逊潮流计算的应用3节点牛顿拉夫逊潮流计算在电力系统规划、运行和故障分析中具有重要的应用价值。
通过计算系统各节点的电压和相角，可以评估系统的电压稳定性和潮流分布，指导电力系统的规划和调度。
在系统发生故障时，可以利用3节点牛顿拉夫逊潮流计算分析系统的稳定性和可靠性，为故障处理提供依据。
5. 结语3节点牛顿拉夫逊潮流计算是一种重要的电力系统分析方法，Matlab作为一种强大的数学工具，为其提供了理想的支持和实现。
通过Matlab进行3节点牛顿拉夫逊潮流计算，可以高效地进行电力系统分析和计算，为电力系统的规划和运行提供科学依据。
希望本文可以帮助读者更加深入地了解3节点牛顿拉夫逊潮流计算及其在Matlab 中的应用。
6. Matlab3节点牛顿拉夫逊潮流计算的优势利用Matlab进行3节点牛顿拉夫逊潮流计算具有许多优势。
Matlab提供了丰富的数学函数和工具，能够快速高效地进行矩阵运算、迭代求解等操作，极大地简化了算法的实现。
电力系统潮流计算牛顿拉夫迅法与PQ分解法通用MATLAB计算程序
此程序经40余同学使用检验，无误。
这是一个电气狗熬两个礼拜图书馆的成果，根据华中科技大学《电力系统分析》中原理编写，可用牛顿-拉夫逊和PQ分解法计算给定标幺值条件的潮流。
本人水平有限，仅供参考，欢迎一起找Bug。
2018/07/06 说明：由于本人变压器建模与PSASP不同，本人使用模型如下图，参数输入时请按该模型计算。
2018/06/18 主程序更新：增加补偿电容参数主程序% file name:chaoliu_lj.m% auther: 山东科技大学罗江% function:使用牛顿-拉夫逊法、PQ分解法计算潮流% updata:2018/6/18 13:22 增加补偿电容参数%节点类型标号%PQ节点 1%PV节点 2%slack节点 3%能计算给定标幺值网络，有且仅有一个平衡节点的潮流%注意：母线标号顺序要求：PQ节点-PV节点-平衡节点%若某元件不存在，其导纳为0，阻抗为infclear %清除工作空间变量clc %清屏%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%数据输入（标幺值）SB=100; %基准容量，单位MVA%母线基准电压Bus=[115 10.5 115 115];%交流线参数：I侧母线J侧母线阻抗1/2接地导纳Line=[4 1 0.06125+0.09527i 0;4 3 0.08469+0.12703i 0;1 3 0.13989+0.15501i 0];%变压器参数：I侧母线J侧母线阻抗变比%变压器阻抗归算到I侧Trans=[2 3 0.0137+0.2881i 0.9504];%加接地电容器补偿: 母线导纳Cap=[2 0.5i];%发电机参数：母线节点类型P V/U θGen=[4 3 1 0];%负荷参数：母线节点类型P Q%按参考方向，发电机发出功率（正值），负荷消耗功率（负值）Load=[1 1 -0.18 -0.06;2 1 -0.32 -0.12];mode=1; %1-极坐标下牛拉法，2-PQ分解法Tmax=50; %最大迭代次数limit=1.0e-4; %要求精度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%变压器π型等效阻抗参数Zt=zeros(size(Trans,1),3);Zt(:,1)=Trans(:,3)./Trans(:,4);Zt(:,2)=Trans(:,3)./(1-Trans(:,4));Zt(:,3)=Trans(:,3)./(Trans(:,4).^2-Trans(:,4));Trans_pi=[Trans(:,1:2) Zt(:,1) 1./Zt(:,2) 1./Zt(:,3)];n=numel(Bus); %总节点数m=n-1; %PQ节点数for i=1:size(Gen,1)%数组行数if Gen(i,2)==2 %除去PV节点就是PQ节点m=m-1;endendfor i=1:size(Load,1)if Load(i,2)==2m=m-1;endend%PQ节点包含浮游节点，其PQ=0%提取P,Q,U向量P=zeros(1,n); %P,Q为原始数据，Pi,Qi为计算结果Q=zeros(1,n);U=ones(1,n); %电压初始值由此确定cita=zeros(1,n); %此处未知节点皆设为1.0∠0 %注意：此处角度单位为度，提取后再转换成弧度，后面计算使用弧度for i=1:size(Gen,1)if Gen(i,2)==1 %PQ节点P(Gen(i,1))=Gen(i,3);Q(Gen(i,1))=Gen(i,4);endif Gen(i,2)==2 %PV节点P(Gen(i,1))=Gen(i,3);U(Gen(i,1))=Gen(i,4);endif Gen(i,2)==3 %slack节点U(Gen(i,1))=Gen(i,3);cita(Gen(i,1))=Gen(i,4);endendfor i=1:size(Load,1)if Load(i,2)==1 %PQ节点P(Load(i,1))=Load(i,3);Q(Load(i,1))=Load(i,4);endif Load(i,2)==2 %PV节点P(Load(i,1))=Load(i,3);U(Load(i,1))=Load(i,4);endif Load(i,2)==3 %slack节点U(Load(i,1))=Load(i,3);cita(Load(i,1))=Load(i,4);endenddisp('初始条件：')disp('各节点有功：')disp(P);disp('各节点无功：')disp(Q);disp('各节点电压幅值：')disp(U);cita=(deg2rad(cita)); %角度转换成弧度disp('各节点电压相角(度)：')disp(rad2deg(cita)); %显示依然使用角度%节点导纳矩阵的计算Y=zeros(n); %新建节点导纳矩阵y=zeros(n); %网络中的真实导纳%计算y(i,j)for i=1:size(Line,1) %与交流线联结的真实导纳ii=Line(i,1); jj=Line(i,2);y(ii,jj)=1/Line(i,3);y(jj,ii)=y(ii,jj);endfor i=1:size(Trans_pi,1) %与变压器联结的真实导纳ii=Trans_pi(i,1); jj=Trans_pi(i,2);y(ii,jj)=1/Trans_pi(i,3);y(jj,ii)=y(ii,jj);end%计算y(i,i)for i=1:size(Line,1) %与交流线联结的对地导纳ii=Line(i,1); jj=Line(i,2);y(ii,ii)=y(ii,ii)+Line(i,4);y(jj,jj)=y(jj,jj)+Line(i,4);endfor i=1:size(Trans_pi,1) %与变压器联结的对地导纳ii=Trans_pi(i,1); jj=Trans_pi(i,2);y(ii,ii)=y(ii,ii)+Trans_pi(i,4);y(jj,jj)=y(jj,jj)+Trans_pi(i,5);end%算上补偿电容for i=1:size(Cap,1)ii=Cap(i,1);y(ii,ii)=y(ii,ii)+Cap(i,2);end%由y计算Yysum=sum(y,1); %每一行求和for i=1:nfor j=1:nif i==jY(i,j)=ysum(i);elseY(i,j)=-y(i,j);endendenddisp('节点导纳矩阵：');disp(Y);G=real(Y); %电导矩阵B=imag(Y); %电纳矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%以上为基础数据整理%接下来是牛拉法的大循环%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%计算功率不平衡量[dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita );disp('有功不平衡量：');disp(dP);disp('无功不平衡量：');disp(dQ);for i=1:Tmaxfprintf('第%d次迭代\n',i);%雅可比矩阵的计算if(mode==1)J=Jacobi( n,m,U,cita,B,G,Pi,Qi );disp('雅可比矩阵');disp(J);end%求解节点电压修正量if(mode==1)[dU,dcita]=Correct( n,m,U,dP,dQ,J );else[dU,dcita]=PQ_LJ( n,m,dP,dQ,U,B );enddisp('电压、相角修正量：');disp(dU);disp(rad2deg(dcita));%修正节点电压U=U+dU;cita=cita+dcita;disp('节点电压幅值：');disp(U);disp('节点电压相角：');disp(rad2deg(cita));%计算功率不平衡量[dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita );disp('有功不平衡量：');disp(dP);disp('无功不平衡量：');disp(dQ);if (max(abs(dP))<limit && max(abs(dQ))<limit )break;end%ifend%for%迭代结束，判断收敛if (max(abs(dP))<limit && max(abs(dQ))<limit )disp('计算收敛');elsedisp('计算不收敛或未达到要求精度');end%打印功率fprintf('迭代总次数：%d\n', i);disp('节点电压幅值：');disp(U);disp('节点电压相角：');disp(rad2deg(cita));disp('有功计算结果：');disp(Pi);disp('无功计算结果：');disp(Qi);子程序一% filename:Unbalanced.m% author: 山东科技大学罗江% function: 计算功率不平衡量function [ dP,dQ,Pi,Qi ] = Unbalanced( n,m,P,Q,U,G,B,cita )%计算ΔPi有功的不平衡量for i=1:nfor j=1:nPn(j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));endPi(i)=sum(Pn);enddP=P(1:n-1)-Pi(1:n-1); %dP有n-1个%计算ΔQi无功的不平衡量for i=1:nfor j=1:nQn(j)=U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));endQi(i)=sum(Qn);enddQ=Q(1:m)-Qi(1:m); %dQ有m个end%func子程序二% filename:Jacobi.m% author:山东科技大学罗江% function: 计算雅可比矩阵function [ J ] = Jacobi( n,m,U,cita,B,G,Pi,Qi )%雅可比矩阵的计算%分块H N K L%i!=j时for i=1:n-1for j=1:n-1H(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));endendfor i=1:n-1for j=1:mN(i,j)=-U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));endendfor i=1:mfor j=1:n-1K(i,j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));endendfor i=1:mfor j=1:mL(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));endend%i==j时for i=1:n-1H(i,i)=U(i).^2*B(i,i)+Qi(i);endfor i=1:mN(i,i)=-U(i).^2*G(i,i)-Pi(i);endfor i=1:mK(i,i)=U(i).^2*G(i,i)-Pi(i);endfor i=1:mL(i,i)=U(i).^2*B(i,i)-Qi(i);end%合成雅可比矩阵J=[H N;K L];end子程序三% filename:Correct.m% author:山东科技大学罗江% function:修正节点电压function [ dU,dcita ] = Correct( n,m,U,dP,dQ,J )%求解节点电压修正量for i=1:mUd2(i,i)=U(i);enddPQ=[dP dQ]';dUcita=(-inv(J)*dPQ)';dcita=dUcita(1:n-1);dcita=[dcita 0];dU=(Ud2*dUcita(n:n+m-1)')';dU=[dU zeros(1,n-m)];end子程序四% filename:PQ_LJ.m% author:山东科技大学罗江% function:使用PQ分解法计算电压修正量function [ dU,dcita ] = PQ_LJ( n,m,dP,dQ,U,B )dP_U=dP./U(1:n-1);dQ_U=dQ./U(1:m);dUdcita=(-inv(B(1:n-1,1:n-1))*dP_U')';dcita=dUdcita./U(1:n-1);dU=(-inv(B(1:m,1:m))*dQ_U')';dU=[dU zeros(1,n-m)];dcita=[dcita 0];%补零end （使用时此括号删去。
matlab潮流程序
clcdisp('此程序为牛拉法解潮流')nPQ=input('请输入PQ节点的个数：');nPV=input('请输入PV节点的个数：');n=nPQ+nPV+1;Ps=[0;-0.5;0.2];Qs=[0;-0.3];Us=[1.0+j*0;1.0+j*0;1.05+j*0;1.05+j*0];% nl nr R X Bl Br zdata=[1 2 0 0.1880 -0.6815 0.6040;1 3 0.1302 0.2479 0.0129 0.0129;1 4 0.1736 0.3306 0.0172 0.0172;3 4 0.2603 0.4959 0.0259 0.0259];% nx Bxdata=[2 0.05];dPQU=1;%计算导纳矩阵nl=zdata(:,1);nr=zdata(:,2);R=zdata(:,3);X=zdata(:,4);Bl=zdata(:,5);Br=zdata(:,6);nx=xdata(:,1);Bx=xdata(:,2);nbr=length(zdata(:,1));nbrx=length(xdata(:,1));Z=R+j*X;y=ones(nbr,1)./Z;Y=zeros(n,n);%计算非对角元素for ii=1:nbrY(nl(ii),nr(ii))= Y(nl(ii),nr(ii))-y(ii);Y(nr(ii),nl(ii))= Y(nl(ii),nr(ii));end%计算对角元素for ii=1:nfor jj=1:nbrif nl(jj)==ii|nr(jj)==iiY(ii,ii)=Y(ii,ii)+y(jj);endendendfor ii=1:nbrY(nl(ii),nl(ii))= Y(nl(ii),nl(ii))+j*Bl(ii);Y(nr(ii),nr(ii))= Y(nr(ii),nr(ii))+j*Br(ii);endfor ii=1:nbrxY(nx(ii),nx(ii))=Y(nx(ii),nx(ii))+j*Bx(ii);end%分离G、BG=real(Y);B=imag(Y);disp('导纳矩阵：');Ye=real(Us);f=imag(Us);k=0;while dPQU>0.00001%求dPdP=zeros(n-1,1);for ii=1:n-1t=0;for jj=1:nt=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));enddP(ii)=Ps(ii)-real(t*(e(ii)+j*f(ii)));end%求dQdQ=zeros(nPQ,1);for ii=1:nPQt=0;for jj=1:nt=t+conj(Y(ii,jj))*(e(jj)-j*f(jj));enddQ(ii)=Qs(ii)-imag(t*(e(ii)+j*f(ii)));end%求dU^2dU2=zeros(nPV,1);ii=1:nPV;dU2(ii)=abs(Us(ii+nPQ))^2-abs(e(ii+nPQ)+j*f(ii+nPQ))^2;dS=[dP;dQ;dU2];dPQU=max(abs(dS));if(dPQU>0.00001)k=k+1%形成雅克比行列式Jacob=zeros(2*(n-1),2*(n-1));%P部分for ii=1:n-1mid1=0;mid2=0;for jj=1:nif ii~=jj&&jj<nJacob(ii,2*jj-1)=-(G(ii,jj)*e( ii)+B(ii,jj)*f(ii));Jacob(ii,2*jj)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);endmid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj );mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj );endJacob(ii,2*ii-1)=-mid2-G(ii,ii)*e(ii)-B(ii,ii) *f(ii);Jacob(ii,2*ii)=-mid1+B(ii,ii)*e(ii)-G(ii,ii)*f (ii);end%Q部分for ii=1:nPQmid1=0;mid2=0;for jj=1:nif ii~=jj&&jj<nJacob(ii+n-1,2*jj-1)=B(ii,jj)* e(ii)-G(ii,jj)*f(ii);Jacob(ii+n-1,2*jj)=G(ii,jj)*e( ii)+B(ii,jj)*f(ii);endmid1=mid1+G(ii,jj)*f(jj)+B(ii,jj)*e(jj );mid2=mid2+G(ii,jj)*e(jj)-B(ii,jj)*f(jj );endJacob(ii+n-1,2*ii-1)=mid1+B(ii,ii)*e(ii)-G(ii, ii)*f(ii);Jacob(ii+n-1,2*ii)=-mid2+G(ii,ii)*e(ii)+B(ii,i i)*f(ii);end%U2部分for ii=nPQ+1:n-1Jacob(ii+n-1,2*ii-1)=-2*e(ii);Jacob(ii+n-1,2*ii)=-2*f(ii);enddU=-inv(Jacob)*dS;de=zeros(n-1,1);df=zeros(n-1,1);ii=1:n-1;de(ii)=dU(2*ii-1);df(ii)=dU(2*ii);e(ii)=e(ii)+de(ii);f(ii)=f(ii)+df(ii);endend%迭代结束U=e+j*f;%计算PV节点的QP=zeros(n,1);Q=zeros(n,1);for ii=1:nPVt=0;for jj=1:nt=t+conj(Y(ii+nPQ,jj))*(e(jj)-j*f(jj));endQ(ii+nPQ)=imag(t*(e(ii+nPQ)+j*f(ii+nPQ)));end%计算平衡节点t=0;for jj=1:nt=t+conj(Y(n,jj))*(e(jj)-j*f(jj));endP(n)=real(t*(e(n)+j*f(n)));Q(n)=imag(t*(e(n)+j*f(n)));ii=1:n-1;P(ii)=Ps(ii);ii=1:nPQ;Q(ii)=Qs(ii);%计算线路潮流Sij=zeros(nbr,1);Sji=zeros(nbr,1);dSij=zeros(nbr,1);for ii=1:nbrSij(ii)=U(nl(ii))*(conj(U(nl(ii)))*(-j*Bl(ii))+(conj(U((nl(ii) )))-conj(U((nr(ii)))))*conj(y(ii)));Sji(ii)=U(nr(ii))*(conj(U(nr(ii)))*(-j*Br(ii))+(conj(U((nr(ii))))-conj(U((nl(ii)))))*conj(y(ii)));dSij(ii)=Sij(ii)+Sji(ii);endnn=[1:n]';disp(' n e f P Q');Display1=[nn e f P Q]disp(' nl nrSij SjidSij');Display2=[nl nr Sij Sji dSij]。
牛顿—拉夫逊法潮流计算MATLAB程序
牛顿一拉夫逊法潮流计算程序By Yuluo%牛顿--拉夫逊法进行潮流计算n=i nput(' 请输入节点数：n=');n1=i nput('请输入支路数：n仁')；isb=i nput(' 请输入平衡母线节点号：isb=');pr=i nput('请输入误差精度：pr=');B1=input('请输入由支路参数形成的矩阵：B1=');B2=input('请输入各节点参数形成的矩阵：B2=');X=input('请输入由节点参数形成的矩阵：X=');Y=zeros( n);e=zeros(1, n);f=zeros(1, n);V=seros(1, n);O=zeros(1, n);S1=zeros( n1); for i=1: nif X(i,2)~=0;p=X(i,1);丫(p,p)=1./X(i,2);endendfor i=1: n1if B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);end丫(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5);Y(p,q)=Y(p,q);Y(p,q)=Y(q,q)+1./(B1(i,3)*B1(i,5F2)+B1(i,4)./2;丫(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;end % 求导纳矩阵G=real(Y);B=imag(Y);for i=1: ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfor i=1: nS(i)=B2(i,1)-B2(i,2);B(i,i)=B(i,i)+B2(i,5);endP=rea(S);Q=imag(S);ICT1=0;IT2=1;NO=2* n;N=NO+1;a=0;while IT2~=0IT2=0;a=a+1;for i=1: n;C(i)=0;D(i)=0;for j1=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);endP仁C(i)*e(i)+f(i)*D(i);Q仁f(i)*C(i)-D(i)*e(i); % 求'P,Q'V2=e(i)A2+f(i)A2;if B2(i,6)~=3DP=P(i)-P1;DQ=Q(i)-Q1;for j1=1: nif j1~=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; end end else DP=P(i)-P1;DV=V(i)~2-V2;for j1=1: nif j1~=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;elseif j1==i&j1~=isbX仁-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; end end end endend % 求雅可比矩阵for k=3:N0k1=k+1;N 1=N;for k2=k1:N1J(k,k2)=J(k,k2)./J(k,k);endJ( k,k)=1;k4=k-1;for k3=3:k4for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J (k, k2);endJ(k3,k)=0;endendfor k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J (k, k2);endJ(k3,k)=0;endendendfor k=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J (k,N);k1=k+1;f(L)=f(L)-J(k1,N);endfor k=3:N0DET=abs (J(k, N));if DET>=prIT2=IT2+1endendICT2(a)=IT2ICT1=ICT1+1;for k=1: ndy(k)=sqrt(e(k)A2+f(k)A2);endfor i=1: nDy(k)=sqrt(e(k)A2+f(kF2);endfor i=1: nDy(ICT1,i)=dy(i);endend % 用高斯消去法解“ w=-J*V”disp('迭代次数');disp(ICTI);disp('没有达到精度要求的个数');disp(ICT2);for k=1: nV(k)=sqrt(e(k)A2+f(k)A2);O(k)=ata n(f(k)./e(k))*180./pi;endE=e+f*j;disp('各节点的实际电压标么值E为(节点号从小到大的排列)：’)；disp(E);disp('各节点的电压大小V为(节点号从小到大的排列)：’)；disp(V);disp('各节点的电压相角0为(节点号从小到大的排列)：’)；disp(O);for p=1: nC(p)=0;for q=1: nC(p)=C(p)+conj(丫(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp('各节点的功率S为(节点号从小到大排列)：’)；disp(S);disp('各条支路的首端功率Si为(顺序同您输入B1时一样)：‘)；for i=1: n1if B1 ( i,6)==0p=B1(i,1);q=B1(i,2);else p=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)))*con j(1./(B1(i,3)*B1(i,5))));disp(Si(p.q));enddisp('各条支路的末端功率Sj为(顺序同您的输入B1时一样)：‘)；for i=1: n1if B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(xonj(E(q)./B1(i,5))-conj(E(p)))*xo nj(1./(B1(i,3)*B1(i,5))));disp(Sj(q,p));enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一样)：'；for i=1: n1if B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endDS(i)=Si(p,q)+Sj(q,p);disp(DS(i));endfor i=1:ICT1Cs(i)=i;enddisp('以下是每次迭代后各节点的电压值(如图所示) ‘)；plot(Cs,Dy),xlabel('迭代次数')，ylabel('电压'),title(' 电压迭代次数曲线');。
基于MATLAB的牛顿拉夫逊迭代法计算潮流(附加短路计算)
这个程序可以适用于三机九节点系统（参数见主程序），本来是想编一个通用各种结构的程序的，但是因为鄙人比较懒，老师留作业时候又没说要通用，就没改完。
惭愧啊。
不过大同小异啦。
有兴趣的慢慢改吧。
使用方法：按后文中给出的代码建立.m文件放在一个文件夹里面。
先运行Powerflow_main.m计算算例系统的潮流；然后运行ShortcircuitCalc.m计算算例系统三相短路电流；程序说明详见各.m文件注释部分，写的已经很详细了，慢慢看吧。
Powerflow_main.m文件代码如下：clear%牛顿拉夫逊迭代法计算潮流format short %规定参数数据显示精度%节点参数矩阵%第一列为节点编号%第二列表示有功注入P%第三列表示无功注入Q%第四列表示电压幅值U%第五列表示电压角度θ%第六列表示发电机x′%第七列表示发电机E′%第八列表示节点类型（2表示平衡节点，1表示PV节点，0表示PQ节点）Node_p=[ 1, 0, 0, 1.04 , 0, 0.3, 1.137, 2;2, 1.63, 0, 1.025, 0, 0.3, 1.211, 1;3, 0.85, 0, 1.025, 0, 0.3, 1.047, 1;4, 0, 0, 1.0, 0, 0, 0, 0;5, -1.25, -0.5, 1.0, 0, 0, 0, 0;6, -0.9, -0.3, 1.0, 0, 0, 0, 0;7, 0, 0, 1.0, 0, 0, 0, 0;8, -1, -0.35, 1.0, 0, 0, 0, 0;9, 0, 0, 1.0, 0, 0, 0, 0];count_s=0;countPV=0;for k=1:size(Node_p,1)if Node_p(k,8)==1countPV=countPV+1;else if Node_p(k,8)==2count_s=count_s+1;end;end;end;countPV;count_s;countPQ=size(Node_p,1)-1-countPV;%显示节点参数disp('节点参数如下：')disp(Node_p)%支路参数%第一列为首节点，第二列为末节点，第三列表示R，第四列表示X，第五列表示B/2 %第六列表示支路类型（1为变比为1的变压器元件；2为输电线元件;0为接地支路）Branch_p =[ 1, 4, 0 , 0.0576, 0 , 1;2, 7, 0 , 0.0625, 0 , 1;3, 9, 0 , 0.0586, 0 , 1;4, 5, 0.01 , 0.085 , 0.088 , 2;4, 6, 0.017 , 0.092 , 0.079 , 2;5, 7, 0.032 , 0.161 , 0.153 , 2;6, 9, 0.039 , 0.17 , 0.179 , 2;7, 8, 0.0085, 0.072 , 0.0745, 2;8, 9, 0.0119, 0.1008, 0.1045, 2];%显示支路参数disp('支路参数如下：')disp(Branch_p)%设置节点初值U=Node_p(:,4);e_ang=Node_p(:,5);P=Node_p(:,2);Q=Node_p(:,3);save data.matformat long%计算结果数据显示精度%显示节点导纳矩阵admi();disp('节点导纳矩阵Y');sparseYKmax=10; %设置最大迭代次数kaccuracy=10^-7;%设置迭代精度k=0;%迭代次数初始化为零for k1=1:Kmax[dP,dQ,y]=getY(U,e_ang);if max(abs(y))<accuracybreak;end;J=jacob(U,e_ang,dP,dQ);x=-inv(J)*y;de_ang=[0;x(1:8)];dU=x(9:14);u1=U(2:9)*dU.';u=diag(u1);U(4:9)=U(4:9)+u;e_ang=e_ang+de_ang;k=k+1;end;e_ang=e_ang/pi*180;save result.matdisp('迭代次数：')kdisp('4号节点至9号节点电压幅值如下：') disp(U(4:9))disp('2号节点至9号节点电压相角如下：') disp(e_ang(2:9))admi.m文件代码如下：%节点导纳矩阵的形成format long %规定数据格式NI=size(Branch_p,1);k_t=1;Y=zeros(NI);for m1=1:NI;I=Branch_p(m1,1);J=Branch_p(m1,2);R=Branch_p(m1,3);X=Branch_p(m1,4);b=Branch_p(m1,5);Style=Branch_p(m1,6);if Style==1 %判断为变压器元件Y(I,I)=Y(I,I)+1/(R+1j*X);Y(J,J)=Y(J,J)+1/(R+1j*X)/k_t/k_t;Y(I,J)=Y(I,J)-1/(R+1j*X)/k_t;Y(J,I)=Y(J,I)-1/(R+1j*X)/k_t;else if Style==0 %判断为母线接地支路元件Y(I,J)=Y(I,J)+1/(R+1j*X);else %判断为输电线元件Y(I,I)=Y(I,I)+1j*b+1/(R+1j*X);Y(J,J)=Y(J,J)+1j*b+1/(R+1j*X);Y(I,J)=Y(I,J)-1/(R+1j*X);Y(J,I)=Y(J,I)-1/(R+1j*X);end;end;end;Y;G=real(Y);B=imag(Y);sparseY=sparse(Y);save data.mat;getY.m文件代码如下：（线性方程组常写作AX=Y形式，故此处命名为get_Y）%计算ΔP和ΔQ的函数function [dP,dQ,y]=getY(U,e_ang)load data P Q G B Node_p countPV count_s;dP=zeros(size(Node_p,1),1); sum1=zeros(size(Node_p,1),1);for i=1:size(Node_p,1)for j=1:size(Node_p,1)sum1(i)=sum1(i)+U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j))+B(i,j)*sin(e_ang(i)-e_ang(j)));end;dP(i)=P(i)-U(i)*sum1(i);end;dQ=zeros(size(Node_p,1),1);sum2=zeros(size(Node_p,1),1);for i=1:size(Node_p,1)for j=1:size(Node_p,1)sum2(i)=sum2(i)+U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j))-B(i,j)*cos(e_ang(i)-e_ang(j)));end;dQ(i)=Q(i)-U(i)*sum2(i);end;y=[dP((count_s+1):9);dQ((countPV+count_s+1):9)]; %拼接ΔP和ΔQ构成方程-J*x=y的向量yjacob.m文件的代码如下：%形成雅克比矩阵的函数function J=jacob(U,e_ang,dP,dQ)load data B G Q P Node_p countPV count_s;size_Y=size(Node_p,1);%H矩阵H1=zeros(size_Y);for i=1:size_Yfor j=1:size_Yif j==i;H1(i,j)= U(i)*U(i)*B(i,j)+Q(i)-dQ(i);elseH1(i,j)=-U(i)*U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j))-B(i,j)*cos(e_ang(i)-e_ang(j)));end;endendH=H1(2:size_Y,2:size_Y);%N矩阵N1=zeros(size_Y);for i=1:size_Yfor j=1:size_Yif j==i;N1(i,j)= -U(i)*U(i)*G(i,j)-(P(i)-dP(i));elseN1(i,j)=-U(i)*U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j))+B(i,j)*sin(e_ang(i)-e_ang(j)));end;end;end;N=N1(2:size_Y,(countPV+count_s+1):size_Y);%M矩阵M1=zeros(size_Y);for i=1:size_Yfor j=1:size_Yif j==i;M1(i,j)=U(i)*U(i)*G(i,j)-(P(j)-dP(j));elseM1(i,j)=U(i)*U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j))+B(i,j)*sin(e_ang(i)-e_ang(j)));end;end;end;M=M1((countPV+count_s+1):size_Y,2:size_Y);%L矩阵L1=zeros(size_Y);for i=1:size_Yfor j=1:size_Yif j==i;L1(i,j)= U(i)*U(i)*B(i,j)-(Q(i)-dQ(i));elseL1(i,j)= -U(i)*U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j))-B(i,j)*cos(e_ang(i)-e_ang(j)));end;endendL=L1((countPV+count_s+1):size_Y,(countPV+count_s+1):size_Y);J=[H N;M L];%拼接构成雅克比矩阵ShortcircuitCalc.m文件的代码如下：clear%三相短路计算format long;%对YN进行修正，形成包括发电机内阻抗和负荷阻抗的节点导纳矩阵re_admi();Z=inv(rY);%计算节点阻抗矩阵disp('4节点发生金属短路')f=4;%短路点为4节点%输出节点阻抗矩阵的短路点所在列disp('节点阻抗矩阵的短路点所在列Z(:,f)=');Z(:,f)%短路电流If计算load result U e_ang;e_angf=e_ang(f);Uf=U(f)*(cos(e_angf/180*pi)+1j*sin(e_angf/180*pi));Zff=Z(f,f);zf=0;If=Uf/(Zff+zf);i_angf=angle(If)*180/pi;If=abs(If);%短路电流计算结果显示disp('短路电流幅值:')Ifdisp('短路电流相角(单位为°)')i_angf%短路时各节点电压计算U_k=U;for x1=1:size(Node_p,1)U_k(x1)=U(x1)*(cos(e_ang(x1)*pi/180)+1j*sin(e_ang(x1)*pi/180))-Z(x1,f)*If*(cos(i_angf*pi/180)+ 1j*sin(i_angf*pi/180));end;%Uk为短路时各节点电压幅值Uk=abs(U_k);%uk_ang为短路时各节点电压相角（单位为°）uk_ang=angle(U_k)*180/pi;%输出短路时各节点电压disp('短路时1-9节点电压：')disp('幅值：')Ukdisp('相角：(单位为°)')uk_ang%计算短路时各支路电流%输出矩阵初始化%第一列为首节点i%第二列为末节点j%第三列为Ii幅值%第四列为Ii相角(单位为°)%第五列为Ij幅值%第六列为Ij相角（单位为°）Ik=[1 4 1 0 1 0;2 7 1 0 1 0;3 9 1 0 1 0;4 5 1 0 1 0;4 6 1 0 1 0;5 7 1 0 1 0;6 9 1 0 1 0;7 8 1 0 1 0;8 9 1 0 1 0];load data Branch_p k_t;%其中变压器变比为k_tNI=size(Branch_p,1);for m1=1:NI;I=Branch_p(m1,1);J=Branch_p(m1,2);R=Branch_p(m1,3);X=Branch_p(m1,4);b=Branch_p(m1,5);Style=Branch_p(m1,6);if Style~=1 %判断为非变压器支路Ik(m1,3)=(U_k(I)-U_k(J))/(R+1j*X)+U_k(I)*1j*b;Ik(m1,4)=angle(Ik(m1,3))*180/pi;Ik(m1,3)=abs(Ik(m1,3));Ik(m1,5)=(U_k(J)-U_k(I))/(R+1j*X)+U_k(J)*1j*b;Ik(m1,6)=angle(Ik(m1,5))*180/pi;Ik(m1,5)=abs(Ik(m1,5));else %否则为变压器支路Ik(m1,3)=(U_k(I)-U_k(J))/(R+1j*X)/k_t;Ik(m1,4)=angle(Ik(m1,3))*180/pi;Ik(m1,3)=abs(Ik(m1,3));Ik(m1,5)=(U_k(J)-U_k(I))/(R+1j*X)/k_t;Ik(m1,6)=angle(Ik(m1,5))*180/pi;Ik(m1,5)=abs(Ik(m1,5));endend%输出短路时各节点电流disp('第一列为首节点i')disp('第二列为末节点j')disp('第三列为Ii幅值')disp('第四列为Ii相角(单位为°)')disp('第五列为Ij幅值')disp('第六列为Ij相角（单位为°）')format shortIksave result_k%形成包括发电机内阻抗、负荷阻抗的节点导纳矩阵format long %规定数据显示精度load data P Q Y Node_p ;load result U;rY=Y;for n1=1:size(Node_p,1)Xd=Node_p(n1,6);Ed=Node_p(n1,7);style=Node_p(n1,8);if style~=0 %判断是否为发电机节点rY(n1,n1)=rY(n1,n1)+1/(1j*Xd);%YN中与发电机节点对应的对角线元素增加发电机导纳else %其他节点(包含负荷节点）rY(n1,n1)=rY(n1,n1)+(-P(n1)+1j*Q(n1))/(U(n1)*U(n1));%与负荷节点对应的对角线元素增加负荷导纳%对于非负荷节点Y矩阵元素不做修改但仍满足上式end;end;rY;。
直角坐标牛顿-拉夫逊法潮流计算matlab程序(仅供参考)
%该程序仅针对《电力系统分析下》何仰赞P61的4节点算例。
%节点电压用直角坐标表示时的牛顿-拉夫逊法潮流计算（matlab程序，可能有小错误，仅供学习之用，如果想要通用的程序，可自己动手改或再到pudn、csdn等网站搜索更好的）%南昌大学电力061李圣涛2009年5月编写，2012年5月上传clc%清空command windowsclear all%清空workspace%为了提高可移植性、可读性、通用性，设置以下变量N=4;%独立节点数NPQ=2;%PQ节点数NPV=1;%PV节点数K=0;%迭代次数%请输入最大迭代次数Kmax。
可从0开始，以观察第Kmax次迭代的结果Kmax=input('\n\n请输入最大迭代次数后回车(可从零开始) Kmax=\n');small=10^(-5);%ε不能太小%i为节点标号，其中1号……NPQ号为PQ节点，（NPQ+1）%号……（N-1）号为PV节点，N号节点为平衡节点%节点导纳矩阵Y的实部G=[1.042093-0.5882350-0.453858;-0.5882351.0690050-0.480769;0000;-0.453858-0.48076900.934627 ];%节点导纳矩阵Y的虚部B=[ -8.2428762.3529413.6666671.891074;2.352941-4.72737702.403846;3.6666670-3.33333330;1.8910742.4038460-4.261590 ];%Y矩阵Y=complex(G,B);%给定PQ节点的Pnode、Qnode，PV节点的Pnode、Vnode。
（Vnode为节点电压的幅值）Pnode=[ -0.3-0.550.5];%PQ、PV节点的初值PQnode=[-0.18-0.130];%PQ节点的初值QVnode=[ 001.10];%PV节点的初值V%迭代初值e=[ 1.01.01.11.05];f=[ 0000];%利用for循环来实现多次迭代。
基于牛顿拉夫逊法潮流计算的matlab实验报告(含源程序和结果)
基于牛顿拉夫逊法潮流计算的matlab实验报告一、实验目的和要求1.学习掌握matlab的基本用法2.应用MATLAB语言编写具有一定通用性的牛顿-拉夫逊法潮流计算程序。
要求：（1）潮流计算方法为牛顿-拉夫逊法。
（2）编程语言为MATLAB。
（3）程序具有较强通用性。
二、程序流程图1.程序流程图开始形成节点导纳矩阵输入原始数据，节点重新编号设节点电压初值(0)(0)i ie f，i=1,2…,n,i≠s置迭代次数P=0置节点号i=1计算雅克比矩阵元素按公式计算PQ节点的()k i P∆，()kiQ∆，PV节点的()kiP∆，()2kiU∆求解修正方程式，得()kie∆，()kif∆雅克比矩阵是否已全部形成？求()max||ke∆，()max||kf∆迭代次数P=P+1i=i+1计算各节点电压的新值：(1)()()k k kie e e+=+∆(1)()()k k kif f f+=+∆三、求解问题及其结果1.求解问题：IEEE-美国新英格兰10机39节点测试系统1）系统单线图2）系统参数1）系统容量基准值为100MV A。
2) 负荷数据见表D-1表D-1 负荷数据3）发电机数据见表D-24）线路参数见表D-3LN35: BUS-4接有并联电容器，B 4＝1.0000 LN36: BUS-5接有并联电容器，B 4＝2.00005)变压器参数见表D-4%IEEE-美国新英格兰10机39节点测试系统% 1 2 3 4 5 6% bus volt angle p q typebus=[ 1 1.0000 0.00 0.00 0.00 12 1.0000 0.00 0.00 0.00 13 1.0000 0.00 -3.22 -0.024 14 1.0000 0.00 -5.00 -1.84 15 1.0000 0.00 0.00 0.00 16 1.0000 0.00 0.00 0.00 17 1.0000 0.00 -2.338 -0.84 18 1.0000 0.00 -5.22 -1.76 19 1.0000 0.00 0.00 0.00 110 1.0000 0.00 0.00 0.00 111 1.0000 0.00 0.00 0.00 112 1.0000 0.00 -0.085 -0.88 113 1.0000 0.00 0.00 0.00 114 1.0000 0.00 0.00 0.00 115 1.0000 0.00 -3.20 -1.53 116 1.0000 0.00 -3.29 -0.323 117 1.0000 0.00 0.00 0.00 118 1.0000 0.00 -1.58 -0.30 119 1.0000 0.00 0.00 0.00 120 1.0000 0.00 -6.80 -1.03 121 1.0000 0.00 -2.74 -1.15 122 1.0000 0.00 0.00 0.00 123 1.0000 0.00 -2.475 -1.15 124 1.0000 0.00 -3.08 -0.922 125 1.0000 0.00 -2.24 -0.472 126 1.0000 0.00 -1.39 -0.17 127 1.0000 0.00 -2.81 -0.755 128 1.0000 0.00 -2.06 -0.276 129 1.0000 0.00 -2.835 -0.269 130 1.0475 0.00 2.50 0.00 231 1.0000 0.00 0.00 0.00 332 1.0000 0.00 6.50 1.759 133 1.0000 0.00 6.32 1.0335 134 1.0123 0.00 5.08 0.00 235 1.0493 0.00 6.50 0.00 236 1.0000 0.00 5.60 0.9688 137 1.0278 0.00 5.40 0.00 238 1.0265 0.00 8.30 0.00 239 1.0300 0.00 -1.04 0.00 2];% 1 2 3 4 5 6 7 % line: from bus to bus R, X, G, B/2 Kline=[ 2 1 0.00350 0.04110 0 0.34935 0;39 1 0.00100 0.02500 0 0.37500 0;3 2 0.00130 0.01510 0 0.12860 0;25 2 0.00700 0.00860 0 0.07300 0;4 3 0.00130 0.02130 0 0.11070 0;18 3 0.00110 0.01330 0 0.10690 0;5 4 0.00080 0.01280 0 0.06710 0;14 4 0.00080 0.01290 0 0.06910 0;6 5 0.00020 0.00260 0 0.02170 0;8 5 0.00080 0.01120 0 0.07380 0;7 6 0.00060 0.00920 0 0.05650 0;11 6 0.00070 0.00820 0 0.06945 0;8 7 0.00040 0.00460 0 0.03900 0;9 8 0.00230 0.03630 0 0.19020 0;39 9 0.00100 0.02500 0 0.60000 0;11 10 0.00040 0.00430 0 0.03645 0;13 10 0.00040 0.00430 0 0.03645 0;14 13 0.00090 0.01010 0 0.08615 0;15 14 0.00180 0.02170 0 0.18300 0;16 15 0.00090 0.00940 0 0.08550 0;17 16 0.00070 0.00890 0 0.06710 0;19 16 0.00160 0.01950 0 0.15200 0;21 16 0.00080 0.01350 0 0.12740 0;24 16 0.00030 0.00590 0 0.03400 0;18 17 0.00070 0.00820 0 0.06595 0;27 17 0.00130 0.01730 0 0.16080 0;22 21 0.00080 0.01400 0 0.12825 0;23 22 0.00060 0.00960 0 0.09230 0;24 23 0.00220 0.03500 0 0.18050 0;26 25 0.00320 0.03230 0 0.25650 0;27 26 0.00140 0.01470 0 0.11980 0;28 26 0.00430 0.04740 0 0.39010 0;29 26 0.00570 0.06250 0 0.51450 0;29 28 0.00140 0.01510 0 0.12450 0;4 0 0 0 0 1.0000 0;5 0 0 0 0 2.0000 0;11 12 0.00160 0.04350 0 0 100.60000/100;13 12 0.00160 0.04350 0 0 100.60000/100;30 2 0.00000 0.01810 0 0 102.50000/100 ;31 6 0.00000 0.02500 0 0 107.00000/100 ;32 10 0.00000 0.02000 0 0 107.00000/100 ;34 20 0.00090 0.01800 0 0 100.90000/100 ;33 19 0.00070 0.01420 0 0 107.00000/100 ;35 22 0.00000 0.01430 0 0 102.50000/100 ;36 23 0.00050 0.02720 0 0 100.00000/100 ;37 25 0.00060 0.02320 0 0 102.50000/100 ;38 29 0.00080 0.01560 0 0 102.50000/100 ;20 19 0.00070 0.01380 0 0 106.00000/100] ;计算结果牛顿－拉夫逊法潮流计算结果节点计算结果：n节点节点电压节点相角（角度）节点注入功率1 1.049185 -8.874991 0.000000 + j 0.0000002 1.053167 -6.367180 0.000000 + j 0.0000003 1.041493 -9.207297 -3.220000 + j -0.0240004 1.036574 -10.042585 -5.000000 + j -1.8400005 1.044652 -8.959237 0.000000 + j 0.0000006 1.043883 -8.293104 0.000000 + j 0.0000007 1.032645 -10.342431 -2.338000 + j -0.8400008 1.031177 -10.811816 -5.220000 + j -1.7600009 1.042715 -10.595648 0.000000 + j 0.00000010 1.046426 -6.010476 0.000000 + j 0.00000011 1.044322 -6.792462 0.000000 + j 0.00000012 1.030736 -6.795388 -0.085000 + j -0.88000013 1.042351 -6.675491 0.000000 + j 0.00000014 1.036310 -8.232337 0.000000 + j 0.00000015 1.018517 -8.519794 -3.200000 + j -1.53000016 1.025492 -7.051856 -3.290000 + j -0.32300017 1.032750 -8.077118 0.000000 + j 0.00000018 1.034779 -8.936485 -1.580000 + j -0.30000019 1.044862 -2.382169 0.000000 + j 0.00000020 0.988148 -3.811032 -6.800000 + j -1.03000021 1.024926 -4.596980 -2.740000 + j -1.15000022 1.042650 -0.070512 0.000000 + j 0.00000023 1.032952 -0.245457 -2.475000 + j -1.15000024 1.021125 -6.906503 -3.080000 + j -0.92200025 1.060163 -4.952002 -2.240000 + j -0.47200026 1.052697 -6.205207 -1.390000 + j -0.17000027 1.037683 -8.217337 -2.810000 + j -0.75500028 1.050444 -2.695196 -2.060000 + j -0.27600029 1.050163 0.063077 -2.835000 + j -0.26900030 1.004392 1.594781 6.500000 + j 1.75900031 0.991632 2.892572 6.320000 + j 1.03350032 1.050539 7.797786 5.600000 + j 0.96880033 1.047500 -3.957598 2.500000 + j 1.21117434 1.012300 1.385774 5.080000 + j 1.82635935 1.049300 4.925324 6.500000 + j 2.63756636 1.027800 1.819476 5.400000 + j -0.10822437 1.026500 7.125579 8.300000 + j 0.21422538 1.030000 -10.390696 -1.040000 + j -2.29163939 1.000000 0.000000 5.628660 + j 1.384403线路计算结果：n节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)2 1 1.178698 + j -0.360055 -1.174311 + j -0.360481 0.004386 + j -0.720536 39 1 6.405845 + j -2.096152 -6.361848 + j 2.408287 0.043997 + j 0.3121353 2 -3.633961 + j -0.542613 3.649983 + j 0.446577 0.016021 + j -0.096036 25 2 2.370242 + j -1.109311 -2.328681 + j 0.997356 0.041562 + j -0.1119554 3 -0.750370 + j -0.307172 0.751094 + j 0.080014 0.000724 + j -0.227159 18 3 0.337560 + j -0.663855 -0.337133 + j 0.438599 0.000427 + j -0.225256 5 4 1.635254 + j 0.499000 -1.633054 + j -0.609119 0.002200 + j -0.110119 14 4 2.621711 + j -0.216428 -2.616576 + j 0.150777 0.005135 + j -0.065651 6 5 4.826035 + j -0.675350 -4.821682 + j 0.684607 0.004353 + j 0.0092578 5 -3.178130 + j -1.041836 3.186428 + j 0.998989 0.008297 + j -0.042847 7 6 -4.249274 + j -0.969559 4.259899 + j 1.010657 0.010625 + j 0.04109811 6 3.465003 + j -0.270003 -3.457273 + j 0.209136 0.007730 + j -0.0608668 7 -1.909893 + j -0.196732 1.911274 + j 0.129559 0.001381 + j -0.0671739 8 0.132235 + j 0.116464 -0.131977 + j -0.521432 0.000258 + j -0.404968 39 9 7.617154 + j -1.902126 -7.557438 + j 2.142687 0.059717 + j 0.24056111 10 -3.483660 + j -0.203064 3.488121 + j 0.171352 0.004461 + j -0.03171213 10 -3.008372 + j -0.730489 3.011879 + j 0.688680 0.003508 + j -0.04180914 13 -2.934129 + j -0.411463 2.941429 + j 0.307264 0.007300 + j -0.10419915 14 -0.311115 + j -0.998556 0.312417 + j 0.627891 0.001303 + j -0.37066516 15 2.896296 + j 0.430232 -2.888885 + j -0.531444 0.007411 + j -0.10121217 16 -2.048841 + j 0.950740 2.052282 + j -1.049122 0.003441 + j -0.098383 19 16 4.542969 + j 0.681545 -4.511670 + j -0.625873 0.031300 + j 0.05567221 16 3.324778 + j -0.302389 -3.316338 + j 0.177006 0.008440 + j -0.125383 24 16 0.410793 + j -0.811601 -0.410571 + j 0.744757 0.000222 + j -0.066844 18 17 -1.917560 + j 0.363855 1.920087 + j -0.475208 0.002527 + j -0.111353 27 17 -0.128621 + j 0.132648 0.128754 + j -0.475531 0.000133 + j -0.342884 22 21 6.093176 + j 1.070437 -6.064778 + j -0.847611 0.028398 + j 0.22282623 22 -0.406149 + j -1.116063 0.406824 + j 0.928040 0.000675 + j -0.18802424 23 -3.490793 + j -0.110399 3.516516 + j 0.138837 0.025723 + j 0.02843826 25 -0.771398 + j -0.442881 0.773189 + j -0.111580 0.001791 + j -0.55446027 26 -2.681379 + j -0.887648 2.691475 + j 0.731900 0.010096 + j -0.15574828 26 1.416063 + j -0.565082 -1.408178 + j -0.210747 0.007885 + j -0.77583029 26 1.921038 + j -0.679443 -1.901899 + j -0.248272 0.019138 + j -0.927715 29 28 3.491624 + j -0.395924 -3.476063 + j 0.289082 0.015561 + j -0.106842 4 0 0.000000 + j -1.074485 0.000000 + j 0.000000 0.000000 + j -1.0744855 0 0.000000 + j -2.182596 0.000000 + j 0.000000 0.000000 + j -2.18259611 12 0.018656 + j 0.473066 -0.018327 + j -0.464126 0.000329 + j 0.00894013 12 0.066943 + j 0.423225 -0.066673 + j -0.415874 0.000270 + j 0.00735130 2 7.897633 + j -0.731582 -7.897633 + j 1.860277 0.000000 + j 1.12869531 6 7.506817 + j 1.371343 -7.506817 + j 0.109153 0.000000 + j 1.48049632 10 12.260592 + j 5.296517 -12.260592 + j -2.064007 0.000000 + j 3.23250934 20 5.080000 + j 1.826359 -5.054406 + j -1.314473 0.025594 + j 0.51188633 19 -1.716763 + j 5.348910 1.736896 + j -4.940504 0.020133 + j 0.40840535 22 6.500000 + j 2.637566 -6.500000 + j -1.998477 0.000000 + j 0.63908936 23 1.402814 + j -0.195113 -1.401865 + j 0.246763 0.000949 + j 0.05165037 25 9.586236 + j 0.419689 -9.533808 + j 1.607517 0.052428 + j 2.02720638 29 -12.165903 + j 2.106593 12.280860 + j 0.135062 0.114957 + j 2.24165520 19 -1.745594 + j 0.284473 1.747837 + j -0.240265 0.002242 + j 0.044208结果分析：此程序的运行结果和试验程序给出的结果是一致的。
基于MATLAB的直角坐标下牛顿拉夫逊法潮流计算
节点的给定功率设为 Pis 和 Qis （称为注入功率）。 假定系统中的第 1、2、…、m 节点为 PQ 节点，对其中每一个节点的 N-R
法表达式 F(x)=0[如 Si 0 、 Pi 0 、 Qi 0 ]形式有些下列方程：
n
n
Pi Pis Pi Pis ei (Gije j Bij f j ) fi (Gij f j Bije j ) 0
j 1
Bij f j ) Giiei
Bii fii
J
ii
U
2 i
ei
2ei
U
2 i
fi
2 fi
（2-5）
当 j i 时,矩阵非对角元素为：
Pi ei
Qiij
Jij
Pi f j
2.2 牛顿——拉夫逊法潮流计算计算公式
把牛顿法用于潮流计算，采用直角坐标形式表示的如式（2-2）所示的形式。
其中电压和支路导纳可表示为：
Ui ei jfi Yij Gij jBij
U j ej jf j
Y ij Gij jBij
（2-1）
将上述表示式（2-1）代入（1-1）式的右端，展开并分出实部和虚部，便
Key words: Electric power system; flow calculation; MATLAB
matpower中的牛顿拉夫逊算法代码
matpower中的牛顿拉夫逊算法代码Matpower是一个用于电力系统潮流和稳定性分析的开源软件包。
牛顿拉夫逊算法是Matpower中用于潮流计算的核心算法之一。
在本文中，我们将一步一步回答有关Matpower中牛顿拉夫逊算法的问题，并深入探讨其原理和应用。
第一部分：Matpower简介和潮流计算1.1 Matpower简介Matpower是一款由美国康奈尔大学和加州大学伯克利分校共同开发的Matlab工具箱。
它提供了一套用于电力系统静态分析的功能齐全的工具，包括潮流计算、最优潮流计算、潮流控制和稳定性分析等。
Matpower已经成为电力系统工程师进行各种电力系统分析的首选软件之一。
1.2 潮流计算潮流计算是电力系统中用于计算系统中各节点电压和功率的方法。
它是电力系统故障分析和计算负荷需求等重要的预处理步骤。
潮流计算的目标是求解电压相角和电压幅值，以及电网中的潮流分布。
第二部分：牛顿拉夫逊算法原理2.1 概述牛顿拉夫逊算法是一种迭代方法，用于求解非线性方程组的数值解。
在潮流计算中，我们需要使用牛顿拉夫逊算法来求解非线性的功率方程组。
2.2 第一步：建立功率方程组牛顿拉夫逊算法的第一步是建立纳入功率平衡方程和节点潮流约束条件的功率方程组。
这个方程组通常表示为“F(x)=0”的形式，其中F是一个向量函数，x是一个向量，它包含了待求的电压相角和电压幅值。
2.3 第二步：线性化功率方程组牛顿拉夫逊算法的第二步是线性化功率方程组。
这可以通过Taylor级数展开来实现。
通过在初始点附近进行一阶近似，我们可以将非线性功率方程组转化为一个线性方程组。
2.4 第三步：求解线性方程组线性化的功率方程组可以写成“J(x)Δx=-F(x)”的形式，其中J 是功率方程组的雅可比矩阵，Δx是一个向量，表示电压相角和电压幅值的变化量。
解这个线性方程组的目的是找到Δx的近似解。
2.5 第四步：更新变量牛顿拉夫逊算法的第四步是更新变量。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
%本程序的功能是用牛顿——拉夫逊法进行潮流计算% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳% 5、支路的变比；6、支路首端处于K侧为1，1侧为0% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；% 3为PV节点；clear;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);sida=zeros(1,n);S1=zeros(nl);% % %---------------------------------------------------for i=1:nl%支路数if B1(i,6)==0%左节点处于1侧p=B1(i,1);q=B1(i,2);else %左节点处于K侧p=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;%对角元K侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;%对角元1侧end%求导纳矩阵disp('导纳矩阵Y=');disp(Y)%----------------------------------------------------------G=real(Y);B=imag(Y);%分解出导纳阵的实部和虚部for i=1:n%给定各节点初始电压的实部和虚部e(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);%PV节点电压给定模值endfor i=1:n%给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量end%===================================================================P=real(S);Q=imag(S); %分解出各节点注入的有功和无功功率ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; %迭代次数ICT1、a；不满足收敛要求的节点数IT2while IT2~=0 % N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列IT2=0;a=a+1;for i=1:nif i~=isb%非平衡节点C(i)=0;D(i)=0;for j1=1:nC(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)endP1=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)%求i节点有功和无功功率P',Q'的计算值V2=e(i)^2+f(i)^2;%电压模平方%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素=========if B2(i,6)~=3%非PV节点DP=P(i)-P1;%节点有功功率差DQ=Q(i)-Q1; %节点无功功率差%=============== 以上为除平衡节点外其它节点的功率计算=================%================= 求取Jacobi矩阵===================for j1=1:nif j1~=isb&j1~=i%非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/deX3=X2; % X2=dp/df X3=dQ/deX4=-X1; % X1=dP/de X4=dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1; % X3=dQ/de J(p,N)=DQ节点无功功率差J(m,q)=X1;J(m,N)=DP;q=q+1; % X1=dP/de J(m,N)=DP节点有功功率差J(p,q)=X4;J(m,q)=X2; % X4=dQ/df X2=dp/dfelseif j1==i&j1~=isb%非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△PJ(m,q)=X2;endendelse%=============== 下面是针对PV节点来求取Jacobi矩阵的元素===========DP=P(i)-P1;% PV节点有功误差DV=V(i)^2-V2;% PV节点电压误差for j1=1:nif j1~=isb&j1~=i%非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/dfX5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; % PV节点电压误差m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; % PV节点有功误差J(m,q)=X2;elseif j1==i&j1~=isb%非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; % PV节点电压误差m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; % PV节点有功误差J(m,q)=X2;endendendendend%========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差=====================for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N; % N=N0+1 即N=2*n+1扩展列△P、△Q 或△Ufor k2=k1:N1% 从k+1列的Jacobi元素到扩展列的△P、△Q 或△UJ(k,k2)=J(k,k2)./J(k,k);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化endJ(k,k)=1; % 对角元规格化K行K列对角元素赋1%==================== 回代运算=======================================if k~=3 % 不是第三行k > 3k4=k-1;for k3=3:k4% 用k3行从第三行开始到当前行的前一行k4行消去for k2=k1:N1% k3行后各行上三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）end %用当前行K2列元素减去当前行k列元素乘以第k行K2列元素J(k3,k)=0; %当前行第k列元素已消为0endif k==N0 %若已到最后一行break;end%================== 前代运算==================================for k3=k1:N0 % 从k+1行到2*n最后一行for k2=k1:N1 % 从k+1列到扩展列消去k+1行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算end %用当前行K2列元素减去当前行k列元素乘以第k行K2列元素J(k3,k)=0; %当前行第k列元素已消为0endelse %是第三行k=3%====================== 第三行k=3的前代运算======================== for k3=k1:N0 %从第四行到2n行（最后一行）for k2=k1:N1 %从第四列到2n+1列（即扩展列）J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）end %用当前行K2列元素减去当前行3列元素乘以第三行K2列元素J(k3,k)=0; %当前行第3列元素已消为0endendend%====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====for k=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N); %修改节点电压实部k1=k+1;f(L)=f(L)-J(k1,N); %修改节点电压虚部end%------修改节点电压-----------for k=3:N0DET=abs(J(k,N));if DET>=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足要求的节点数加1endendICT2(a)=IT2; %不满足要求的节点数ICT1=ICT1+1; %迭代次数end%用高斯消去法解"w=-J*V"disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2); %计算各节点电压的模值sida(k)=atan(f(k)./e(k))*180./pi; %计算各节点电压的角度E(k)=e(k)+f(k)*j; %将各节点电压用复数表示end%=============== 计算各输出量===========================disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E); %显示各节点的实际电压标幺值E用复数表示disp('-----------------------------------------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V); %显示各节点的电压大小V的模值disp('-----------------------------------------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida); %显示各节点的电压相角for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q)); %计算各节点的注入电流的共轭值endS(p)=E(p)*C(p); %计算各节点的功率S = 电压X 注入电流的共轭值enddisp('各节点的功率S为(节点号从小到大排列)：');disp(S); %显示各节点的注入功率disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Si(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))));Siz(i)=Si(p,q);elseSi(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))));Siz(i)=Si(p,q);enddisp(Si(p,q));SSi(p,q)=Si(p,q);ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);elseSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);enddisp(Sj(q,p));SSj(q,p)=Sj(q,p);ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);disp('-----------------------------------------------------');enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i));DDS(i)=DS(i);ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);disp('-----------------------------------------------------');endfigure(1);subplot(1,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(1,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;figure(2);subplot(2,2,1);P=real(S);Q=imag(S);bar(P);xlabel('节点号');ylabel('节点注入有功');grid on;subplot(2,2,2);bar(Q);xlabel('节点号');ylabel('节点注入无功');grid on;P1=real(Siz);Q1=imag(Siz);bar(P1);xlabel('支路号');ylabel('支路首端注入有功'); grid on;subplot(2,2,4);bar(Q1);xlabel('支路号');ylabel('支路首端注入无功'); grid on;。