牛拉法潮流计算
ieee333节点牛拉法潮流计算结果
ieee333节点牛拉法潮流计算结果潮流计算是电力系统分析中的一项重要工作,用于确定系统中各节点的电压幅值和相角的分布情况。
本文将以IEEE 333节点系统为例,使用牛拉法潮流计算方法,对该系统进行潮流计算,并给出计算结果。
IEEE 333节点系统是一个中等规模的电力系统,包含333个节点。
在进行潮流计算之前,我们需要确定系统中的各个节点的发电机有功和无功注入,以及负载的有功和无功消耗。
注入和消耗的功率值可以通过实际测量或者根据电力系统数据获得。
假设我们已经获取了这些信息,下面将进行潮流计算。
潮流计算的主要目标是确定系统中各节点的电压幅值和相角。
潮流计算可分为以下几个步骤:1.建立雅可比矩阵潮流计算的第一步是建立雅可比矩阵。
雅可比矩阵描述了节点电压和注入功率之间的关系。
在IEEE 333节点系统中,节点电压表示为复数形式,即幅值和相角。
雅可比矩阵的大小由系统的节点数决定,对于333节点系统,雅可比矩阵的大小为333x333。
2.初始化节点电压和功率不平衡在开始潮流计算之前,需要初始化节点电压和功率不平衡。
初始化时,可以假设节点电压的幅值为1,相角为0度。
同时,初始化功率不平衡为初始负荷值。
3.迭代计算节点电压和功率不平衡通过迭代计算的方式,逐步更新节点电压和功率不平衡,直到收敛为止。
在每一次迭代计算中,通过雅可比矩阵和牛拉法方程来更新节点电压和功率不平衡。
4.收敛判断和结果分析在迭代计算过程中,需要判断潮流计算是否收敛。
通常使用节点电压和功率不平衡的变化情况来判断收敛性。
当节点电压和功率不平衡的变化小于预定的阈值时,可以认为潮流计算已经收敛。
此时,可以得到系统中各节点的电压幅值和相角。
通过对IEEE 333节点系统进行潮流计算,可以得到系统中各节点的电压幅值和相角分布情况。
这些结果对电力系统的运行和规划具有重要意义,可以用于判断系统的稳定性和对系统进行优化。
值得注意的是,潮流计算是一项复杂而繁琐的工作,需要进行大量的计算和数据处理。
牛拉法潮流计算
%本程序的功能是用牛拉法进行潮流计算%原理介绍详见鞠平著《电气工程》%默认数据为鞠平著《电气工程》例8.4所示数据%B1是支路参数矩阵%第一列和第二列是节点编号。
节点编号由小到大编写%对于含有变压器的支路,第一列为低压侧节点编号,第二列为高压侧节点编号%第三列为支路的串列阻抗参数,含变压器支路此值为变压器短路电抗%第四列为支路的对地导纳参数,含变压器支路此值不代入计算%第五烈为含变压器支路的变压器的变比,变压器非标准电压比%第六列为变压器是否是否含有变压器的参数,其中“1”为含有变压器,“0”为不含有变压器%B2为节点参数矩阵%第一列为节点注入发电功率参数%第二列为节点负荷功率参数%第三列为节点电压参数%第四列%第五列%第六列为节点类型参数,“1”为平衡节点,“2”为PQ节点,“3”为PV节点参数%X为节点号和对地参数矩阵%第一列为节点编号%第二列为节点对地参数clear;clc;num=input('是否采用默认数据?(1-默认数据;2-手动输入)');if num==1n=4;n1=4;isb=4;pr=0.00001;B1=[1 2 0.1667i 0 0.8864 1;1 3 0.1302+0.2479i 0.0258i 1 0;1 4 0.1736+0.3306i 0.0344i 1 0;3 4 0.2603+0.4959i 0.0518i 1 0];B2=[0 0 1 0 0 2;0 -0.5-0.3i 1 0 0 2;0.2 0 1.05 0 0 3;0 -0.15-0.1i 1.05 0 0 1];X=[1 0;2 0.05i;3 0;4 0];elsen=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:X=');endTimes=1; %迭代次数%创建节点导纳矩阵Y=zeros(n);for i=1:n1if B1(i,6)==0 %不含变压器的支路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(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else %含有变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-B1(i,5)/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+B1(i,5)/B1(i,3)+(1-B1(i,5))/B1(i,3);Y(q,q)=Y(q,q)+B1(i,5)/B1(i,3)+(B1(i,5)*(B1(i,5)-1))/B1(i,3);endendfor i=1:n1Y(i,i)=Y(i,i)+X(i,2); %计及补偿电容电纳enddisp('导纳矩阵为:');disp(Y); %显示导纳矩阵%初始化OrgS、DetaSOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%创建OrgS,用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2 %不是平衡点&是PQ点h=h+1;for j=1:n%公式8-74%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j ,3)));OrgS(2*h,1) =OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:n %对PV节点的处理,注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3 %不是平衡点&是PV点h=h+1;for j=1:n%公式8-75-a%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j ,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendend%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endend%创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); %delPiDetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1); %delQi endendt=0;for i=1:n %对PV节点的处理,注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1); %delPiDetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2; %delUiendend% DetaS%创建I,用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));%conj求共轭endend%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理if B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1) %将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0;endendendendendk=0;for i=1:n %对PV节点的处理if B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0;endendendendenddisp('初始雅可比矩阵为:');disp(Jacbi);%求解修正方程,获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %inv矩阵求逆% DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理if B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n %对PV节点的处理if B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endend% B2%开始循环********************************************************************** while abs(max(DetaU))>prOrgS=zeros(2*n-2,1);h=0;j=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j ,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:nif i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j ,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendend% OrgS%创建DetaSh=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,6)==3h=h+1;t=t+1;% DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;endend% DetaS%创建Ii=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endend% I%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendend% JacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;% DetaU%修正节点电压j=0;for i=1:nif B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:nif B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endend% B2Times=Times+1; %迭代次数加1enddisp('迭代次数为:');disp(Times);disp('收敛时电压修正量为::');disp(DetaU);for k=1:nE(k)=B2(k,3);e(k)=real(E(k));f(k)=imag(E(k));V(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;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); %显示各节点的注入功率Sline=zeros(n1,5);disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致):');for i=1:n1p=B1(i,1);q=B1(i,2);Sline(i,1)=B1(i,1);Sline(i,2)=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))*((1-B1(i,5))/B1(i,3))+(conj(E(p))-conj(E(q)))*(B1(i,5)/B1(i,3)));Siz(i)=Si(p,q);endSSi(p,q)=Si(p,q);Sline(i,3)=Siz(i);ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);enddisp('-----------------------------------------------------');disp('各条支路的末端功率Sj为(顺序同您输入B1时一致):');for i=1:n1p=B1(i,1);q=B1(i,2);if B1(i,6)==0Sj(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))*((B1(i,5)*(B1(i,5)-1))/B1(i,3))+(conj(E(q))-conj(E(p)))*(B1(i,5)/B1(i,3)));Sjy(i)=Sj(q,p);endSSj(q,p)=Sj(q,p);Sline(i,4)=Sjy(i);ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);enddisp('-----------------------------------------------------');disp('各条支路的功率损耗DS为(顺序同您输入B1时一致):');for i=1:n1p=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);DDS(i)=DS(i);Sline(i,5)=DS(i);ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);enddisp('-----------------------------------------------------');disp('各支路首端编号末端编号首端功率末端功率线路损耗');disp(Sline);。
牛顿拉夫逊法计算潮流步骤
牛顿拉夫逊法计算潮流步骤牛顿拉夫逊法(Newton-Raphson Method)是一种常用于计算潮流的数值求解方法。
它是基于潮流计算的功率流方程的非线性特性而设计的,通过迭代求解来逼近潮流计算的稳态解。
下面将介绍牛顿拉夫逊法计算潮流的基本步骤。
首先,我们需要明确潮流计算的目标,即确定电力系统中各节点的电压相角和幅值。
这些节点是电力系统中的发电机、负荷和交流输电线路的连接点。
通过潮流计算,我们可以得到各节点的电压相角和幅值,从而分析系统的功率分布、电压稳定性等运行特性。
接下来,我们需要建立电力系统的潮流计算模型。
这个模型中,我们需要考虑发电机的注入功率、负荷的吸收功率、线路的传输损耗等因素。
通过利用功率流方程,我们可以将这些因素表示为电压、功率和导纳之间的方程。
然后,我们需要进行初始化操作。
在进行牛顿拉夫逊法迭代计算之前,我们需要对电力系统的各节点进行初始电压值的设定。
这些初始值可以根据经验或者历史数据来得到,但需要满足物理约束条件,如一致性、电压幅值在合理范围内等。
接下来,我们进入迭代计算的过程。
首先,我们需要对系统的节点进行编号,然后选择某一节点作为基准节点,其他节点相对于基准节点的电压相角进行计算。
然后,我们根据节点注入功率和导纳矩阵的关系,得到节点注入电流。
接着,我们根据节点注入电流和电压相角的关系,计算各节点的电压相角和幅值的改变量。
在计算改变量后,我们需要对节点电压进行更新。
更新后,我们判断系统是否达到收敛条件。
如果满足收敛条件,则停止迭代,得到最终的潮流计算结果;如果不满足收敛条件,则继续进行下一轮迭代计算。
最后,我们对潮流计算结果进行分析和验证。
通过比较计算得到的结果和实际运行数据进行对比,我们可以评估潮流计算的准确性。
同时,我们还可以通过故障分析、电压稳定性评估等手段对电力系统进行优化和改进。
总而言之,牛顿拉夫逊法是一种常用的求解潮流计算问题的方法。
它通过迭代求解潮流计算的功率流方程,逼近潮流计算的稳态解。
牛拉法潮流计算
自动化07-1班段佳function nl;%------------------------------------------------------------------------%===================================================================%======================牛顿——拉夫逊法==============================%===========================潮流计算=================================%===================================================================%-----------------------------------------------------------------------% % %---------------使用说明部分---------------------------display('% %本程序的功能是用牛顿——拉夫逊法进行潮流计算');display('% %本程序要求用户按照一定的格式将电力系统的参数制成excel表格,系统运行时将从excel中加载这些参数,随后后即可进行潮流计算');display('% %为了方便运算,用户再给系统节点进行编号时,请按照先PQ节点,再PV节点,最后平衡节点的顺序从小到大编号');display('% %电力系统潮流计算excel格式——支路参数:1、支路首端号;2、末端号;3、支路阻抗;4、支路对地电纳;5、支路的变比K:1;6、支路首端处于K侧为1,1侧为0'); display('% %电力系统潮流计算excel格式——节点参数:1、节点号;2、电压大小;3、相位角;4、发电机有功;5、发电机无功;6、负载有功;7、负载无功;8、节点类型'); %===================================================================%==============================数据准备==============================%===================================================================% %---------------------电力系统数据加载部分-----------------------------------------------clearx=0;Branch=0;%支路参数Note=0;%节点参数[filename, pathname] = uigetfile('*.xls', 'please choose the excel file with your powersystem parameters ');%从外部excel导入电力系统潮流计算相关参数tryif filename ~= 0x=xlsread([pathname,filename],'sheet1', 'A3:F3');Branch=xlsread([pathname,filename],'sheet1', 'A5:G10');%读支路参数Note=xlsread([pathname,filename],'sheet1', 'A15:H19');%读节点参数endcatch%进行出错处理errmsg = lasterr;errordlg(errmsg,'Save as Error');rethrow(lasterror);end% %---------------------支路参数初始化部分-----------------------------------------------SB=100;UB=220;n=1;m=1;pr=0.0001;SB=x(5);%功率基准值UB=x(6);%电压基准值n=x(1);%节点数nl=x(2);%支路数m=x(3);%PQ节点的个数pr=x(4);;%误差精度B1(:,1)=Branch(:,1);%1、支路首端号B1(:,2)=Branch(:,2);%2、末端号B1(:,3)=Branch(:,3)+Branch(:,4)*i;%3、支路阻抗B1(:,4)=Branch(:,5)*i;%4、支路对地电纳B1(:,5)=Branch(:,6);%5、支路的变比K:1;B1(:,6)=Branch(:,7);%6、支路首端处于K侧为1,1侧为0'% %% %---------------------节点参数初始化部分--------------------------------------------------U=ones(n,1);a=zeros(n,1);Ps=zeros(n,1);Qs=zeros(n,1);P=zeros(n,1);Q=zeros(n,1);detp=zeros(n-1,1);detq=zeros(m,1);deta=zeros(n-1,1);detu=zeros(m,1);k=0;%迭代次数U=Note(:,2);%各节点电压初始值(标幺值)a=Note(:,3);%各节点电压相位初始值(弧度)Gp=Note(:,4);%各节点发电机有功功率初始值(标幺值)Gq=Note(:,5);%各节点发电机无功功率初始值(标幺值)Lp=Note(:,6);%各节点负载有功功率初始值(标幺值)Lq=Note(:,7);%各节点负载无功功率初始值(标幺值)type=Note(:,8);%节点类型,PQ节点=1 ,PV节点=2 ,平衡节点=3for h=1:nPs(h)=Gp(h)-Lp(h);%各节点注入的有功功率Qs(h)=Gq(h)-Lq(h);%各节点注入的无功功率end% % %---------------------导纳矩阵计算部分-----------------------------------------------------Y=zeros(n);for h=1:nl %支路数if B1(h,6)==0 %左节点处于低压侧(6、支路首端处于K侧为1,1侧为0)p=B1(h,1);q=B1(h,2); %1、支路首端号;2、末端号;Y(p,q)=Y(p,q)-1./B1(h,3); %非对角元 3、支路阻抗;4、支路对地电纳;5、支路的变比;Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1./B1(h,3)+B1(h,4);Y(q,q)=Y(q,q)+1./B1(h,3)+B1(h,4);elsep=B1(h,1);q=B1(h,2); %1、支路首端号;2、末端号;Y(p,q)=Y(p,q)-1./(B1(h,3)*B1(h,5));%非对角元 3、支路阻抗;4、支路对地电纳;5、支路的变比;Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1./B1(h,3)+B1(h,4);Y(q,q)=Y(q,q)+1./(B1(h,3)*B1(h,5)^2)+B1(h,4);endend%导纳矩阵显示disp('导纳矩阵 Y=');disp(Y)% % %-------------OK,至此潮流计算所需的数据已经准备好了----------------%===================================================================%==============================潮流计算==============================%===================================================================%u(i)=e(i)+jf(i);Y(ij)=G(ij)+jB(ij);G=real(Y);B=imag(Y);%分解出导纳阵的实部和虚部%============================计算失配功率初始值detp\detq==========================for h=1:n-1s=0;for j=1:ns=s+U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));endP(h)=U(h)*s;endfor h=1:n-1s=0;for j=1:ns=s+U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));endQ(h)=U(h)*s;endfor h=1:n-1detp(h)=Ps(h)-P(h);endfor h=1:mdetq(h)=Qs(h)-Q(h);end%============================不满足精度要求则进入循环========================== while(max(abs(detp))>=pr|max(abs(detq))>=pr)%%不满足精度要求则循环%=================================求取Jacobi矩阵===============================H=zeros(n-1,n-1);N=zeros(n-1,m);K=zeros(m,n-1);L=zeros(m,m);for h=1:n-1for j=1:n-1if h==jH(h,j)=U(h)^2*B(h,j)+Q(h);elseH(h,j)=-U(h)*U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));endendendfor h=1:n-1for j=1:mif h==jN(h,j)=-U(h)^2*G(h,j)-P(h);elseN(h,j)=-U(h)*U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));endendendfor h=1:mfor j=1:n-1if h==jK(h,j)=U(h)^2*G(h,j)-P(h);elseK(h,j)=U(h)*U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));endendendfor h=1:mfor j=1:mif h==jL(h,j)=U(h)^2*B(h,j)-Q(h);elseL(h,j)=-U(h)*U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));endendend%========================解修正方程,得到修正量detu,deta============================Jacobi=[H N;K L];display(Jacobi);dets=[detp;detq];solutions=-inv(Jacobi)*dets;deta=solutions(1:n-1,:);detu=solutions(n:n-1+m,:);%==============================迭代过程中的电压====================================for h=1:n-1a(h)=a(h)+deta(h);endfor h=1:mU(h)=U(h)+detu(h);endk=k+1;fprintf('迭代次数k=%d\n',k);disp('节点电压大小(标幺值)');disp(U);disp('节点电压相位角(弧度)');disp(a);%===========================迭代过程中的失配功率detp\detq===========================for h=1:n-1s=0;for j=1:ns=s+U(j)*(G(h,j)*cos(a(h)-a(j))+B(h,j)*sin(a(h)-a(j)));endP(h)=U(h)*s;endfor h=1:n-1s=0;for j=1:ns=s+U(j)*(G(h,j)*sin(a(h)-a(j))-B(h,j)*cos(a(h)-a(j)));endQ(h)=U(h)*s;endfor h=1:n-1detp(h)=Ps(h)-P(h);endfor h=1:mdetq(h)=Qs(h)-Q(h);enddisp('迭代过程中的有功失配功率(标幺值)');disp(detp);disp('迭代过程中的无功失配功率(标幺值)');disp(detq);end% % %-------------OK,至此潮流计算已经完成了----------------%===================================================================%==============================计算结果输出到工作区========================== %===================================================================%=================================迭代次数、各节点电压和视在功率==============================disp('计算结果');fprintf('总的迭代次数k=%d\n',k);disp('-----------------------------------------------------');disp('各节点电压大小(标幺值)为(节点号从小到大排列)');disp(U);disp('各节点电压相位角(角度)为(节点号从小到大排列)');A=a*180/pi;disp(A);disp('-----------------------------------------------------');disp('各节点视在功率(标幺值)为(节点号从小到大排列)');S=P+Q*i;disp(S);%=============================各条支路功率损耗和总损耗========================= ZSH=0;DS=zeros(nl,1);for h=1:nlp=B1(h,1);q=B1(h,2);DS(h)=S(p)-S(q);ZSH=ZSH+DS(h);DDS(h)=DS(h)*SB;ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(h)),' (MVA) 标么值:',num2str(DS(h))];disp(ZF);enddisp('-----------------------------------------------------');disp(['总损耗为:ZSH=',num2str(ZSH*SB),' (MVA) 标么值:',num2str(ZSH)]);%==============================结果输出到原excel========================== result0=U;%电压result1=A;%相位result2=P;%节点有功result3=Q;%节点无功result4=real(DS);%线路有功损耗result5=imag(DS);%线路无功损耗result6=real(ZSH);%系统总有功损耗result7=imag(ZSH);%系统总无功损耗[filename1, pathname1] = uiputfile('*.xls', 'put the result into the excel with your powersystem parameters ');%从外部excel导入电力系统潮流计算相关参数tryif filename1 ~= 0xlswrite([pathname1,filename1],result0 , 'sheet1', 'J3');xlswrite([pathname1,filename1],result1 , 'sheet1', 'K3');xlswrite([pathname1,filename1],result2 , 'sheet1', 'L3');xlswrite([pathname1,filename1],result3 , 'sheet1', 'M3');xlswrite([pathname1,filename1],result4 , 'sheet1', 'N3');xlswrite([pathname1,filename1],result5 , 'sheet1', 'O3');xlswrite([pathname1,filename1],result6 , 'sheet1', 'P3');xlswrite([pathname1,filename1],result7 , 'sheet1', 'Q3');endcatch%进行出错处理errmsg = lasterr;errordlg(errmsg,'Save as Error');rethrow(lasterror);end%==============================打开excel查看计算结果========================== winopen([pathname1,filename1]);% % %-------------OK,至此潮流计算已经全部完成----------------% % %-------------O(∩_∩)O哈!----------------% %本程序的功能是用牛顿——拉夫逊法进行潮流计算% %本程序要求用户按照一定的格式将电力系统的参数制成excel表格,系统运行时将从excel中加载这些参数,随后后即可进行潮流计算% %为了方便运算,用户再给系统节点进行编号时,请按照先PQ节点,再PV节点,最后平衡节点的顺序从小到大编号% %电力系统潮流计算excel格式——支路参数:1、支路首端号;2、末端号;3、支路阻抗;4、支路对地电纳;5、支路的变比K:1;6、支路首端处于K侧为1,1侧为0% %电力系统潮流计算excel格式——节点参数:1、节点号;2、电压大小;3、相位角;4、发电机有功;5、发电机无功;6、负载有功;7、负载无功;8、节点类型导纳矩阵 Y=6.3110 -20.4022i -3.5587 +11.3879i -2.7523 + 9.1743i 0 0-3.5587 +11.3879i 8.5587 -31.0093i -5.0000 +15.0000i 0 + 4.9889i-2.7523 + 9.1743i -5.0000 +15.0000i 7.7523 -28.7757i 0 0 + 4.9889i0 0 + 4.9889i 0 0 - 5.2493i0 0 0 + 4.9889i 0 0 - 5.2493iJacobi =-20.5622 11.3879 9.1743 0 -6.3110 3.5587 2.752311.3879 -31.3768 15.0000 4.9889 3.5587 -8.5587 5.00009.1743 15.0000 -29.1632 0 2.7523 5.0000 -7.75230 4.9889 0 -4.9889 0 0 06.3110 -3.5587 -2.7523 0 -20.2422 11.3879 9.1743-3.5587 8.5587 -5.0000 0 11.3879 -30.6418 15.0000-2.7523 -5.0000 7.7523 0 9.1743 15.0000 -28.3882迭代次数k=1节点电压大小(标幺值)1.00361.02971.03271.00001.0000节点电压相位角(弧度)-0.0900-0.0577-0.06120.0425迭代过程中的有功失配功率(标幺值)0.0193-0.0059-0.0007-0.0140迭代过程中的无功失配功率(标幺值)-0.0148-0.0387-0.0270Jacobi =-21.0649 11.6431 9.4218 0 -5.5312 4.0561 3.1247 11.8809 -32.9618 15.9694 5.1115 3.2952 -9.0811 5.2601 9.5859 15.9316 -30.6597 0 2.5776 5.3736 -8.2678 0 5.1115 0 -5.1115 0 -0.5140 0 7.1808 -4.0561 -3.1247 0 -20.0304 11.6431 9.4218 -3.2952 9.0693 -5.2601 -0.5140 11.8809 -32.7992 15.9694 -2.5776 -5.3736 8.2664 0 9.5859 15.9316 -30.7138迭代次数k=2节点电压大小(标幺值)0.99371.02071.02401.00001.0000节点电压相位角(弧度)-0.0904-0.0587-0.06200.0397迭代过程中的有功失配功率(标幺值)0.0019-0.0010-0.0008-0.0002迭代过程中的无功失配功率(标幺值)0.0046-0.0041-0.0041Jacobi =-20.6812 11.4294 9.2518 0 -5.4239 3.9739 3.0648 11.6585 -32.4210 15.6951 5.0675 3.2410 -8.9173 5.1741 9.4110 15.6605 -30.1704 0 2.5340 5.2777 -8.1299 0 5.0675 0 -5.0675 0 -0.5002 0 7.0387 -3.9739 -3.0648 0 -19.6079 11.4294 9.2518 -3.2410 8.9154 -5.1741 -0.5002 11.6585 -32.1891 15.6951 -2.5340 -5.2777 8.1284 0 9.4110 15.6605 -30.1786迭代次数k=3节点电压大小(标幺值)0.99351.02031.02371.00001.0000节点电压相位角(弧度)-0.0905-0.0587-0.06200.0397迭代过程中的有功失配功率(标幺值)1.0e-004 *0.6037-0.2358-0.3406-0.0364迭代过程中的无功失配功率(标幺值)1.0e-003 *0.1808-0.1108-0.1555Jacobi =-20.6709 11.4238 9.2471 0 -5.4240 3.9719 3.0632 11.6527 -32.4023 15.6839 5.0657 3.2395 -8.9101 5.1705 9.4062 15.6495 -30.1528 0 2.5328 5.2739 -8.1234 0 5.0657 0 -5.0657 0 -0.5000 0 7.0351 -3.9719 -3.0632 0 -19.6066 11.4238 9.2471 -3.2395 8.9101 -5.1705 -0.5000 11.6527 -32.1625 15.6839 -2.5328 -5.2739 8.1233 0 9.4062 15.6495 -30.1531迭代次数k=4节点电压大小(标幺值)0.99351.02031.02361.00001.0000节点电压相位角(弧度)-0.0905-0.0587-0.06200.0397迭代过程中的有功失配功率(标幺值)1.0e-005 *0.1199-0.0264-0.0861-0.0075迭代过程中的无功失配功率(标幺值)1.0e-005 *0.3521-0.1563-0.3785计算结果总的迭代次数k=4-----------------------------------------------------各节点电压大小(标幺值)为(节点号从小到大排列)0.99351.02031.02361.00001.0000各节点电压相位角(角度)为(节点号从小到大排列)-5.1825-3.3648-3.55382.2723-----------------------------------------------------各节点视在功率(标幺值)为(节点号从小到大排列)-0.8055 - 0.5320i0.0000 - 0.1200i0.0000 + 0.0000i0.5000 + 0.1837iDS(1,2)=-80.5501-41.2005i (MVA) 标么值:-0.8055-0.41201iDS(1,3)=-80.5502-53.2007i (MVA) 标么值:-0.8055-0.53201iDS(2,3)=-5.97263e-005-12.0002i (MVA) 标么值:-5.9726e-007-0.12i DS(4,2)=50+30.3696i (MVA) 标么值:0.5+0.3037iDS(5,3)=-8.6117e-005-0.i (MVA) 标么值:-8.6117e-007-3.7855e-006i -----------------------------------------------------总损耗为:ZSH=-111.1005-76.03228i (MVA) 标么值:-1.111-0.76032i。
牛顿拉夫逊法计算潮流步骤
牛顿拉夫逊法计算潮流步骤牛顿拉夫逊法(Newton-Raphson method)是一种用于求解非线性方程组的迭代方法,它可以用来计算电力系统潮流的解。
潮流计算是电力系统规划和运行中的重要任务,它的目标是求解电力系统中各节点的电压幅值和相角,以及线路的功率流向等参数,用于分析电力系统的稳定性和安全性,以及进行电力系统规划和调度。
下面是使用牛顿拉夫逊法计算潮流的一般步骤:步骤1:初始化首先,需要对电力系统的各个节点(包括发电机节点和负荷节点)的电压幅值和相角进行初始化,一般可以使用其中一种估计值或者历史数据作为初始值。
步骤2:建立潮流方程根据电力系统的潮流计算模型,可以建立节点电压幅值和相角的平衡方程,一般采用节点注入功率和节点电压的关系来表示。
潮流方程一般是一个非线性方程组,包含了各个节点之间的复杂关系。
步骤3:线性化方程组将潮流方程组进行线性化处理,一般采用泰勒展开的方法,将非线性方程组变为线性方程组。
线性化的过程需要计算雅可比矩阵,即方程组中的系数矩阵。
步骤4:求解线性方程组利用线性方程组的求解方法,比如高斯消元法或LU分解法等,求解线性方程组,得到电压幅值和相角的修正量。
步骤5:更新节点电压根据线性方程组的解,更新各个节点的电压幅值和相角,得到新的节点电压。
步骤6:检查收敛性判断节点电压的修正量是否小于设定的收敛阈值,如果满足收敛条件,则停止迭代;否则,返回步骤3,循环进行线性化方程组和线性方程组的求解。
步骤7:输出结果当潮流计算收敛时,输出最终的节点电压幅值和相角,以及线路的功率流向等参数。
牛顿拉夫逊法是一种高效、快速且收敛性良好的方法,广泛应用于电力系统潮流计算。
在实际应用中,可能会遇到多次迭代或者收敛性不好的情况,此时可以采用退火技术或其他优化算法进行改进。
此外,牛顿拉夫逊法的计算也可以并行化,利用多核处理器或者分布式计算集群来加速计算过程。
总之,牛顿拉夫逊法是一种重要的潮流计算方法,通过迭代计算逼近非线性方程组的解,可以得到电力系统中各节点的电压幅值和相角,用于分析电力系统的稳定性和安全性。
牛顿-拉夫逊法潮流计算
目录摘要11.设计意义与要求2 1.1设计意义21.2设计要求32.牛顿—拉夫逊算法3 2.1牛顿算法数学原理:32.2 直角坐标系下牛顿法潮流计算的原理43 详细设计过程10 3.1节点类型103.2待求量103.3导纳矩阵103.4潮流方程113.5修正方程124.程序设计15 4.1 节点导纳矩阵的形成154.2 计算各节点不平衡量164.3 雅克比矩阵计算- 19 -4.4 LU分解法求修正方程- 22 -4.5 计算网络中功率分布- 25 -5.结果分析- 25 -6.小结- 29 -参考文献- 30 -附录:- 31 -摘要潮流计算是电力网络设计及运行中最基本的计算,对电力网络的各种设计方案及各种运行方式进行潮流计算,可以得到各种电网各节点的电压,并求得网络的潮流及网络中各元件的电力损耗,进而求得电能损耗。
在数学上是多元非线性方程组的求解问题,求解的方法有很多种。
牛顿—拉夫逊法是数学上解非线性方程式的有效方法,有较好的收敛性。
将牛顿法用于潮流计算是以导纳矩阵为基础的,由于利用了导纳矩阵的对称性、稀疏性及节点编号顺序优化等技巧,使牛顿法在收敛性、占用存、计算速度等方面都达到了一定的要求。
本文以一个具体例子分析潮流计算的具体方法,并运用牛顿—拉夫逊算法求解线性方程关键词:电力系统潮流计算牛顿—拉夫逊算法1.设计意义与要求1.1设计意义潮流计算是电力系统分析中的一种最基本的计算,他的任务是对给定运行条件确定系统运行状态,如各母线上的电压(幅值及相角)、网络中的功率分布及功率损耗等。
潮流计算的结果是电力系统稳定计算和故障分析的基础。
具体表现在以下方面:(1)在电网规划阶段,通过潮流计算,合理规划电源容量及接入点,合理规划网架,选择无功补偿方案,满足规划水平的大、小方式下潮流交换控制、调峰、调相、调压的要求。
(2)在编制年运行方式时,在预计负荷增长及新设备投运基础上,选择典型方式进行潮流计算,发现电网中薄弱环节,供调度员日常调度控制参考,并对规划、基建部门提出改进网架结构,加快基建进度的建议。
牛拉法潮流计算例题
牛拉法潮流计算例题首先,牛拉法潮流计算是一种用于电力系统稳态分析的方法,它可以用来计算电力系统中各个节点的电压幅值和相角,以及各个支路的电流大小和相角。
下面是一个牛拉法潮流计算的例题。
假设有一条简单的电力系统,由三个节点和两条支路组成。
节点1和节点2之间连接一条1欧姆的电阻,节点2和节点3之间连接一条0.5欧姆的电阻。
节点1的电压幅值为1.05千伏,相角为0度,节点3的电压幅值为1千伏,相角为-120度。
现在需要计算节点2的电压幅值和相角,以及两条支路的电流大小和相角,假设电力系统中各个元件均为纯电阻。
首先,我们可以列出节点间的导纳矩阵,其中导纳元素为各个支路的导纳值,节点1和节点2之间的导纳为1欧姆的导纳,节点2和节点3之间的导纳为0.5欧姆的导纳,对角线元素为各自节点所连支路的导纳之和。
接下来,我们需要选择一个节点作为参考节点,假设我们选择节点1作为参考节点。
然后,我们可以将节点电压表示为复数形式,即V1=1.05∠0度,V3=1∠-120度。
由于节点1的电压已知,我们可以将其表示为参考电压,即V1=1∠0度=1+j0。
然后,我们可以利用导纳矩阵和节点电压,求解未知节点的电压和支路电流。
具体地,我们可以列出节点2的电压方程式:I12=(V1-V2)/1I23=(V2-V3)/0.5I12=-I23其中,I12和I23分别是支路12和支路23的电流。
将节点电压表示为复数形式,并带入上式,得到:(V1-V2)/1=(1+j0-V2)/1(V2-V3)/0.5=(V2-1∠-120度)/0.5I12=I23化简上式,可得:V2=1.045-j0.2558I12=0.0045-j0.2558I23=0.0045+j0.1279因此,节点2的电压幅值为1.056千伏,相角为-14.34度,支路12的电流大小为0.2558安,相角为-83.66度,支路23的电流大小为0.1279安,相角为29.74度,计算完成。
牛顿拉夫逊法潮流计算
牛顿拉夫逊法潮流计算牛顿-拉夫逊法(Newton-Raphson method)是一种用于求解非线性方程的数值方法。
它通过迭代逼近根的方式,将非线性方程转化为一系列的线性方程来求解。
在电力系统中,潮流计算用于确定电力网中节点的电压幅值和相角。
潮流计算是电力系统分析的重要基础,可以用于计算电力系统的潮流分布、功率损耗、节点电压稳定度等参数,为电力系统的规划、运行和控制提供参考依据。
牛顿-拉夫逊法是一种常用的潮流计算方法,它的基本思想是通过不断迭代来逼近电网的潮流分布,直到满足一定的收敛条件。
下面将对牛顿-拉夫逊法的具体步骤进行详细介绍。
首先,我们需要建立电力网络的节点潮流方程,即功率方程。
对于每一个节点i,其节点功率方程可以表示为:Pi - Vi * (sum(Gij * cos(θi - θj)) - sum(Bij * sin(θi -θj))) = 0Qi - Vi * (sum(Gij * sin(θi - θj)) + sum(Bij * cos(θi -θj))) = 0其中,Pi和Qi分别为节点i的有功功率和无功功率,Vi和θi分别为节点i的电压幅值和相角,Gij和Bij分别为节点i和节点j之间的导纳和电纳。
接下来,我们需要对每个节点的电压幅值和相角进行初始化。
一般情况下,可以将电压幅值设置为1,相角设置为0。
然后,我们可以开始进行迭代计算。
在每一轮迭代中,我们需要计算每个节点的雅可比矩阵和功率残差,然后更新电压幅值和相角。
雅可比矩阵可以通过对节点功率方程进行求导得到,具体如下:dPi/dVi = -sum(Vj * (Gij * sin(θi - θj) + Bij * cos(θi - θj)))dPi/dθi = sum(Vj * (Gij * Vi * cos(θi - θj) - Bij * Vi * sin(θi - θj)))dQi/dVi = sum(Vj * (Gij * cos(θi - θj) - Bij * sin(θi - θj)))dQi/dθi = sum(Vj * (Gij * Vi * sin(θi - θj) + Bij * Vi * cos(θi - θj)))功率残差可以通过将节点功率方程代入得到,如下:RPi = Pi - Vi * (sum(Gij * cos(θi - θj)) - sum(Bij *sin(θi - θj)))RQi = Qi - Vi * (sum(Gij * sin(θi - θj)) + sum(Bij *cos(θi - θj)))最后,我们可以使用牛顿-拉夫逊法的迭代公式来更新电压幅值和相角,具体如下:Vi(new) = Vi(old) + ΔViθi(new) = θi(old) + Δθi其中,ΔVi和Δθi分别为通过求解线性方程组得到的电压幅值和相角的增量。
牛拉法潮流计算程序(附3机9节点结果对比)
摘要电力系统潮流计算是研究电力系统稳态运行的一种重要方法,它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态,包括各母线的电压、线路的功率分布以及功率损耗等等。
潮流计算主要用于电网规划和静态安全分析,它可为扩建电力网络,以达到规划周期内所需要的输电能力提供依据;也可以对预想事故进行模拟和分析,校核预想事故下的电力系统安全性。
本文简单介绍了牛顿-拉夫逊潮流计算的原理、模型与算法,然后用具体的实例,利用MATLAB对牛顿-拉夫逊法的算法进行了验证。
关键词:电力系统潮流计算牛顿-拉夫逊法 MATLAB一、牛拉法的数学模型对一个N 节点的电力网路,列写节点电压方程,即I =Y V(1.1)式中,I 为节点注入电流列相量,Y 为节点导纳矩阵,V 为节点电压列相量。
由于异地测量的两个电流缺少时间同步信息,以注入功率替换注入电流作为已知量。
即***1+niij j ij j i i i Y V V I V Q P ∙∙===∑(1.2)其中,Y ij =G ij +j B ij ,带入上式,得到有功功率和无功功率方程 P i =V i V j G ij cos θij +B ij sin θij n j =1 (1.3)Q i =V i V j G ij sin θij −B ij cos θij n j =1(1.4)大部分情况下,已知PQ ,求解V θ。
考虑到电网的功率平衡,至少选择一台发电机来平衡全网有功功率,即至少有一个平衡节点,常选择调频或出线较多的发电机作为平衡节点。
具有无功补偿的母线能保持电压幅值恒定,这类节点可作为PV 节点。
潮流计算中节点分类总结如下:已知电力系统有m 个PQ 节点,r 个PV 节点和1个平衡节点,则可以提取m+r 个有功功率方程和m 个无功功率方程,从而求解出m+r 个θ和m 个V ,其余节点的有功和无功可通过式(1.3)、(1.4)求得,这样就完成了潮流计算。
ieee333节点牛拉法潮流计算结果
IEEE 333节点牛拉法潮流计算结果分析一、潮流计算简介潮流计算是电力系统分析的基础之一,通过对电力系统各个节点的电压、功率以及电流等参数进行计算和分析,从而得到电力系统各个节点的电气特性。
潮流计算的结果对于电力系统的稳定运行、负荷分配、设备运行等方面具有重要的指导意义,因此潮流计算一直是电力系统研究和运行中的一个重要课题。
二、IEEE 333节点牛拉法潮流计算IEEE 333节点系统是一个经典的电力系统仿真模型,它包括了发电机、负荷、变压器、输电线路等多种设备,并具有典型的电力系统特性。
针对IEEE 333节点系统进行潮流计算能够充分考察电力系统在不同工作条件下的运行特性,对于电力系统的研究和分析具有重要的参考价值。
在IEEE 333节点系统中,采用了牛拉法潮流计算方法,该方法通过对电力系统各个节点的功率平衡方程和节点电压平衡方程进行求解,从而得到电力系统各个节点的电压、相角、有功和无功功率等参数。
这些计算结果可以直观地反映出电力系统在不同工况下的运行状况,为电力系统的分析和设计提供了重要的数据支持。
三、IEEE 333节点潮流计算结果分析1. 电压分布通过对IEEE 333节点系统进行潮流计算,可以得到不同节点的电压值。
电压是电力系统中非常重要的参数,它直接关系到负载的供电质量和设备的安全运行。
潮流计算结果表明,在IEEE 333节点系统中,各个节点的电压分布相对均匀,没有出现明显的电压偏差,表明该系统在静态稳定方面具有较好的特性。
2. 有功功率分布有功功率是电力系统的重要性能指标,它直接关系到发电机的输出能力和负载的用电需求。
通过潮流计算得到的有功功率分布结果显示,在IEEE 333节点系统中,各个节点的有功功率消耗相对均衡,未出现严重的功率不平衡现象,表明该系统在功率分配方面具有较好的平衡性。
3. 无功功率分布无功功率是电力系统的另一个重要性能指标,它与电力系统的稳定运行和无功补偿设备的运行有着密切的关系。
第三节牛顿 拉夫逊法潮流计算
∂P H11 = 1 = U1 U 2 ( −G12 sin δ12 + B12 cos δ12 ) ∂δ1 +U 3 ( −G13 sin δ13 + B13 cos δ13 ) + ... = −U1 ∑ U j (Gij sin δ ij − Bij cos δ ij )
j =2
PV节点:δi • 节点功率和支路功率(第二求解对象)
4-3 牛顿—拉夫逊法潮流计算
共2(m-1)+(n-m)=n+m-2个变量, 则需n+m-2个独立方程
节点注入功率—电压实数方程组(极坐标形式)
对节点i:
~ & S i = Pi + jQ i = U i
∑
* * Yij U j j =1
~ Si = U i
n
∑ (G
j =1 ij
− jBij U j e
)
jδ ij
e
jδ ij
= cos δ ij + j sin δ ij
∑ (G
j =1
− jBij U j cos δ ij + j sin δ ij
) (
)
4-3 牛顿—拉夫逊法潮流计算
节点注入功率—电压实数方程组(极坐标形式)
j =1
n
n
(
)
)
Qi = U i ∑ U j Gij sin δ ij − Bij cos δ ij
j =1
(
(U,δ)不是真解
∆Pi (U, δ) = Pi − U i ∑U j Gij cosδ ij + Bij sin δ ij
j =1 n
j =1
电力系统牛拉法潮流计算
电力系统牛拉法潮流计算电力系统的潮流计算是电力系统运行中一项重要的工作,它用来确定电力系统中各节点的电压和功率的分布情况。
牛拉法(Newton-Raphson)方法是一种主要的潮流计算方法,它是基于牛顿迭代法的一种改进方法,可以用来求解非线性方程组,并被广泛应用于电力系统的潮流计算。
牛拉法潮流计算的基本原理是通过不断迭代求解节点电压和相应的功率,直到收敛为止。
具体步骤如下:1.建立潮流计算的数学模型。
电力系统的潮流计算可以被建模为一个非线性方程组,其中未知量为各节点的电压和功率,方程组的解表示系统的潮流分布情况。
2.初始化节点电压。
初始时,可以假设所有节点的电压为1,并根据负荷功率和潮流方向,计算各发电机节点的功率注入。
3.计算节点电压。
利用牛拉法迭代求解非线性方程组。
首先,根据电压相角和幅值的变化情况,更新节点电压;然后,利用更新的节点电压计算各发电机节点的功率注入,以及从节点注入到节点之间的功率传输;最后,根据功率平衡方程计算支路的功率。
4.判断迭代是否收敛。
判断迭代是否收敛的常用方法有两个:一是通过计算节点电压变化量来判断,如果变化量小于一定阈值,则认为计算收敛;二是通过计算功率平衡误差来判断,如果误差小于一定阈值,则认为计算收敛。
5.如果迭代未收敛,返回步骤3;如果迭代收敛,计算结束,得到系统的潮流分布情况。
牛拉法潮流计算的优点是能够处理复杂的非线性方程组,收敛速度快,并且适用于大规模电力系统的计算。
但是,牛拉法潮流计算也存在一些问题,比如可能出现发散情况,需要进行故障处理。
牛拉法潮流计算在电力系统调度和运行中起着重要的作用。
通过潮流计算,可以确保电力系统的稳定运行,优化电力系统的运行方式,提高系统的可靠性和经济性。
总结起来,牛拉法潮流计算是电力系统潮流计算的一种重要方法,通过迭代求解非线性方程组,可以得到电力系统各节点的电压和功率的分布情况。
它在电力系统调度和运行中具有重要的应用价值,可以帮助优化电力系统的运行方式,提高系统的稳定性和经济性。
潮流计算的基本算法及使用方法
潮流计算的基本算法及使用方法一、 潮流计算的基本算法1. 牛顿-拉夫逊法1.1 概述牛顿-拉夫逊法是目前求解非线性方程最好的一种方法。
这种方法的特点就是把对非线性方程的求解过程变成反复对相应的线性方程求解的过程,通常称为逐次线性化过程,就是牛顿-拉夫逊法的核心。
牛顿-拉夫逊法的基本原理是在解的某一邻域内的某一初始点出发,沿着该点的一阶偏导数——雅可比矩阵,朝减小方程的残差的方向前进一步,在新的点上再计算残差和雅可矩阵继续前进,重复这一过程直到残差达到收敛标准,即得到了非线性方程组的解。
因为越靠近解,偏导数的方向越准,收敛速度也越快,所以牛顿法具有二阶收敛特性。
而所谓“某一邻域”是指雅可比方向均指向解的范围,否则可能走向非线性函数的其它极值点,一般来说潮流由平电压即各母线电压(相角为0,幅值为1)启动即在此邻域内。
1.2 一般概念对于非线性代数方程组即 ()0,,,21=n i x x x f ()n i ,2,1= (1-1)在待求量x 的某一个初始计算值()0x附件,将上式展开泰勒级数并略去二阶及以上的高阶项,得到如下的线性化的方程组()()()()()0000=∆'+x x f x f (1-2)上式称之为牛顿法的修正方程式。
由此可以求得第一次迭代的修正量()()()[]()()0100x f x f x -'-=∆ (1-3)将()0x∆和()0x相加,得到变量的第一次改进值()1x。
接着再从()1x出发,重复上述计算过程。
因此从一定的初值()0x出发,应用牛顿法求解的迭代格式为()()()()()k k k x f x x f -=∆' (1-4)()()()k k k x x x ∆+=+1 (1-5)上两式中:()x f '是函数()x f 对于变量x 的一阶偏导数矩阵,即雅可比矩阵J ;k 为迭代次数。
由式(1-4)和式子(1-5)可见,牛顿法的核心便是反复形成求解修正方程式。
牛顿拉斐逊法潮流计算
牛顿拉斐逊法潮流计算牛顿拉夫逊法(Newton-Raphson method)是一种数值计算方法,用于解非线性方程。
其原理是通过迭代来逼近方程的根。
在电力系统中,牛顿拉夫逊法常用于求解潮流计算问题。
潮流计算是电力系统调度运行和规划的基础工作,其目的是确定电力系统各节点的电压幅值和相角,以及各支线上的功率和无功功率。
通过潮流计算可以有效地评估电力系统的稳定性和运行状态,并为电力系统的调度和规划提供参考依据。
牛顿拉夫逊法的核心思想是通过接近方程的根来求解非线性方程。
其基本步骤如下:1.初始化:选取一个初始点作为方程的近似解,通常选择电力系统的平衡状态作为初值。
2.构造雅可比矩阵:根据潮流方程的特点,建立牛顿拉夫逊法的雅可比矩阵。
雅可比矩阵描述了非线性方程的导数关系,用于迭代计算过程中的线性化。
3.迭代计算:利用雅可比矩阵和当前解向量,构建迭代格式,并计算得到新的解向量。
迭代格式中,包括牛顿方程和拉夫逊方程。
牛顿方程用于计算不平衡功率的校正量,而拉夫逊方程用于计算不平衡电压的校正量。
4.收敛判断:判断迭代计算得到的新解是否满足收敛条件。
通常使用误差向量的范数作为判断依据。
如果误差向量的范数小于预先设定的阈值,即可认为迭代已经收敛。
5.循环迭代:如果迭代计算得到的新解不满足收敛条件,继续进行迭代计算,直到达到收敛条件为止。
牛顿拉夫逊法的优点是收敛速度较快,尤其适用于求解非线性方程的问题。
然而,该方法也存在一些缺点。
首先,牛顿拉夫逊法需要提供一个合适的初始点,如果初始点选择不当,可能会导致迭代过程发散。
其次,构造雅可比矩阵和计算迭代格式的过程较为复杂,需要一定的数学基础和计算能力。
在电力系统潮流计算中,牛顿拉夫逊法广泛应用于求解节点电压和支路功率的平衡方程。
通过牛顿拉夫逊法,可以准确地计算出系统各节点的电压幅值和相角,指导电网的调度运营和规划工作。
总之,牛顿拉夫逊法是一种重要的数值计算方法,特别适用于求解非线性方程。
第三节牛顿拉夫逊法潮流计算
第三节牛顿拉夫逊法潮流计算牛顿-拉夫逊法(Newton-Raphson method)是一种数值计算方法,用于求解非线性方程和潮流计算问题。
它是基于牛顿迭代法和拉夫逊迭代法的结合,可高效地求解电力系统潮流计算问题。
潮流计算是电力系统运行分析中的重要环节,其目标是确定系统中每个节点的电压和相角,并计算各个支路的电流,以评估系统的功率传输和稳定性。
在传统的高压电力系统中,由于负荷、发电机和传输线等元件的非线性特性,潮流计算问题呈现为非线性的数学方程组,通常采用迭代方法求解。
牛顿-拉夫逊法的基本思想是通过对方程组的线性化近似,迭代求解线性方程组的解,以接近方程组的精确解。
它通过将非线性方程组转化为以下形式进行迭代:F(x)=0其中,F(x)是非线性方程组的向量函数,x是未知向量。
牛顿-拉夫逊法的迭代过程可通过以下步骤进行:1.初始化变量:根据系统的初始状态进行节点电压和相角的初始化。
2.计算雅可比矩阵:通过对非线性方程组进行偏导,得到雅可比矩阵。
雅可比矩阵描述了各个节点潮流量与节点电压和相角之间的关系。
3.迭代计算:通过牛顿迭代法进行迭代计算,直到达到指定的收敛条件。
具体步骤为:a.解线性方程组:根据雅可比矩阵和当前节点电压和相角,求解线性方程组,得到修正量。
b.更新变量:根据修正量和当前节点电压和相角,更新节点电压和相角的值。
c.判断收敛:判断修正量是否满足收敛条件,如果满足则结束迭代计算,否则返回步骤a。
牛顿-拉夫逊法的优点是收敛速度快,精度高。
然而,它的缺点是对于方程组的收敛性和初始值的选择要求较高,存在收敛到局部最小值的问题。
为了克服这些问题,可以采用改进的牛顿-拉夫逊法,如增加松弛因子或采用多起点迭代法等。
总之,牛顿-拉夫逊法是一种高效的求解非线性方程组和潮流计算问题的数值方法。
它在电力系统潮流计算中广泛应用,帮助分析和评估电力系统的稳定性和功率传输能力。
随着电力系统的规模和复杂性的增加,牛顿-拉夫逊法的进一步改进和优化仍然是一个研究的热点问题。
牛拉法潮流计算
牛拉法潮流计算%本程序的功能是用牛拉法进行潮流计算 %原理介绍详见鞠平著《电气工程》%默认数据为鞠平著《电气工程》例8.4所示数据±是支路参数矩阵%第一列和第二列是节点编号。
节点编号由小到大编写%对于含有变压器的支路,第一列为低压侧节点编号,第二列为高压侧节点编号 %第三列为支路的串列阻抗参数,含变压器支路此值为变压器短路电抗 %第四列为支路的对地导纳参数,含变压器支路此值不代入计算 %第五烈为含变压器支路的变压器的变比,变压器非标准电压比%第六列为变压器是否是否含有变压器的参数,其中“1”为含有变压器,“0”为不含有变压器2为节点参数矩阵%第一列为节点注入发电功率参数 %第二列为节点负荷功率参数 %第三列为节点电压参数 %第四列 %第五列%第六列为节点类型参数,“1”为平衡节点,“2”为PQ节点,“3”为PV节点参数%X为节点号和对地参数矩阵 %第一列为节点编号 %第二列为节点对地参数clear; clc;num=input('是否采用默认数据?(1-默认数据;2-手动输入)'); if num==1 n=4; n1=4; isb=4;pr=0.00001;B1=[1 2 0.1667i 0 0.8864 1;1 3 0.1302+0.2479i 0.0258i 1 0;1 40.1736+0.3306i 0.0344i 1 0;3 4 0.2603+0.4959i 0.0518i 1 0];B2=[0 0 1 0 0 2;0 -0.5-0.3i 1 0 0 2;0.2 0 1.05 0 0 3;0 -0.15-0.1i 1.05 0 0 1]; X=[1 0;2 0.05i;3 0;4 0]; elsen=input('请输入节点数:n='); n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb='); pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1='); B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:X='); endTimes=1; %迭代次数%创建节点导纳矩阵 Y=zeros(n); for i=1:n1if B1(i,6)==0 %不含变压器的支路 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(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4); else %含有变压器的支路p=B1(i,1); q=B1(i,2);Y(p,q)=Y(p,q)-B1(i,5)/B1(i,3); Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+B1(i,5)/B1(i,3)+(1-B1(i,5))/B1(i,3);Y(q,q)=Y(q,q)+B1(i,5)/B1(i,3)+(B1(i,5)*(B1(i,5)-1))/B1(i,3); end end for i=1:n1Y(i,i)=Y(i,i)+X(i,2); %计及补偿电容电纳 enddisp('导纳矩阵为:');disp(Y); %显示导纳矩阵%初始化OrgS、DetaS OrgS=zeros(2*n-2,1); DetaS=zeros(2*n-2,1);%创建OrgS,用于存储初始功率参数 h=0; j=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2 %不是平衡点&是PQ点 h=h+1; forj=1:n%公式8-74%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej) %Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i ,j))*real(B2(j,3))); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3))); end endendfor i=1:n %对PV节点的处理,注意这时不可再将h初始化为0 ifi~=isb&B2(i,6)==3 %不是平衡点&是PV点 h=h+1; for j=1:n%公式8-75-a%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej) %Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i ,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3))); end end end%创建PVU 用于存储PV节点的初始电压 PVU=zeros(n-h-1,1); t=0; for i=1:nif B2(i,6)==3 t=t+1;PVU(t,1)=B2(i,3); end end%创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量 h=0;for i=1:n %对PQ节点的处理 if i~=isb&B2(i,6)==2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); TlPiDetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1); TlQi end end t=0;for i=1:n %对PV节点的处理,注意这时不可再将h初始化为0 ifi~=isb&B2(i,6)==3 h=h+1; t=t+1;DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1); TlPiDetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2; TlUi end end % DetaS%创建I,用于存储节点电流参数i=zeros(n-1,1); h=0; for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));%conj求共轭end end%创建Jacbi(雅可比矩阵) Jacbi=zeros(2*n-2); h=0; k=0;for i=1:n %对PQ节点的处理 if B2(i,6)==2 h=h+1; forj=1:n if j~=isb k=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1)); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3)); Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3)); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); endif k==(n-1) %将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0; end end end end end k=0;for i=1:n %对PV节点的处理 if B2(i,6)==3 h=h+1; forj=1:n if j~=isb k=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1)); Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3)); else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3)); Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3)); Jacbi(2*h,2*k-1)=0; Jacbi(2*h,2*k)=0; endif k==(n-1) %将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0; end end end end enddisp('初始雅可比矩阵为:'); disp(Jacbi);%求解修正方程,获取节点电压的不平衡量 DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %inv矩阵求逆 % DetaU%修正节点电压 j=0;for i=1:n %对PQ节点处理 if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end endfor i=1:n %对PV节点的处理 if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end end % B2%开始循环********************************************************************** whileabs(max(DetaU))>pr OrgS=zeros(2*n-2,1); h=0;感谢您的阅读,祝您生活愉快。
电力系统牛拉法潮流计算教案资料
18
牛顿-拉夫逊法潮流计算
牛顿法计算流程 1 初始化,形成节点导纳阵,给出初值 x ( 0 ) 2 令k=0 进入迭代循环
2.1 计算函数值 f ( x ( k ) ) ,判断是否收敛 f (x(k)) 2.2 计算Jacobian矩阵 f ( x(k ) ) 2.3 计算修正量 x(k) ( f(x(k))) 1f(x(k)) 2.4 对变量进行修正 x(k 1)x(k) x(k),k=k+1返回
牛顿法的历史 牛顿法基本原理
对于非线性方程 f (x) 0
给定初值 x ( 0 ) 用Talor级数展开,有:
f(x(0) x(0))f(x(0))f'(x(0)) x(0)f''(x(0)) x(0)L 2!
忽略高阶项,0 则有 f(x(0))f'(x(0)) x(0)0
17
牛顿-拉夫逊法潮流计算
已知量为:平衡节点的电压;除平衡节点外所 有节点的有功注入量;PQ节点的无功注入量; PV节点的电压辐值
直角坐标下和极坐标下有不同的处理方法
9
直角坐标下潮流方程
直角坐标下待求变量
e1
M
x
en f1
M f n
直角坐标下功率方程
P1
M
Pn Q1
f (x) M
2.1
3 输出计算结果
19
牛顿-拉夫逊法潮流计算
牛顿法可写成如下简单迭代格式
x ( k 1 ) x ( k ) ( J ( x ( k ) ) ) 1 f( x ( k ) ) ( x ( k ) )
i1,2,LN
Pi Vi jiVj(Gijcos ij BijsinBij)
牛拉法潮流计算
a= sum(real(Y(2,:)).*e0-imag(Y(2,:)).*f0)*e0(2)+sum(real(Y(2,:)).*f0+imag(Y(2,:)).*e0)*f0(2);
deta_P2=P2s-a;
a= sum(real(Y(2,:)).*e0-imag(Y(2,:)).*f0)*f0(2)-sum(real(Y(2,:)).*f0+imag(Y(2,:)).*e0)*e0(2);
P1s=-0.3;Q1s=-0.18;
P2s=-0.55;Q2s=-0.13; %PQ节点的初值给定
P3s=0.5;V3s=1.1;V4s=1.05; %P V 以及平衡节点的初值给定
k=0;
PQV=[P1s Q1s P2s Q2s P3s V4s];
e=1;%初始化
%进入潮流计算多次循环判定
end
k
e0
f0
v=zeros(1,4);sgm=zeros(1,4);
v=sqrt(e0.^2+f0.^2) %电压幅值
sgm=atan(f0./e0)/pi*180 %电压相位
if rem(i,2)==1&&rem(j,2)==1
J(i,j)=-real(Y((i+1)/2,(j+1)/2))*e0((i+1)/2)-imag(Y((i+1)/2,(j+1)/2))*f0((i+1)/2);
y=[0,0.01528i,0,0.01920j;
0.01528i,0,0,0.01413j;
0,0,0,0;
0.01920i,0.01413i,0,0]; %(各支路的导纳)
两机五节点网络潮流计算—牛拉法
两机五节点网络潮流计算—牛拉法牛拉法(Gauss-Seidel Method)是一种常用的迭代方法,用于解决电力系统的潮流计算问题。
在电力系统中,潮流计算是一项重要的工作,用于求解网络中各节点的电压和功率大小。
牛拉法是一种有效的求解方法,适用于小型电力系统,其基本思想是通过迭代来逼近最优解。
潮流计算问题可以抽象成求解非线性方程组的问题,即求解节点电压复数值的方程组。
具体来说,我们需要求解以下方程组:P_i = V_i * ( G_ii * cosθ_i + ∑(G_ij * cos(θ_i - θ_j)) - B_ii * sinθ_i - ∑(B_ij * sin(θ_i - θ_j)))Q_i = V_i * ( G_ii * sinθ_i + ∑(G_ij * sin(θ_i - θ_j)) + B_ii * cosθ_i + ∑(B_ij * cos(θ_i - θ_j)))其中,P_i和Q_i分别表示第i个节点的有功功率和无功功率,V_i表示第i个节点的电压幅值,θ_i表示第i个节点的电压相角,G_ij和B_ij分别表示节点i和节点j之间的导纳和电纳。
牛拉法的基本思路是通过迭代,逐步逼近节点电压的最优解。
假设我们需要求解的是一个两机五节点网络。
首先,我们可以随机初始化每个节点的电压幅值和相角值(也可以根据经验给定初始值)。
然后,根据上述方程组,计算每个节点的有功功率和无功功率。
接下来,我们采用牛拉法的迭代步骤来逼近节点电压的最优解。
具体步骤如下:1.选择一个初始节点(可以是任意节点),将其电压相角θ_i固定为0。
2.通过方程组计算该节点的电压幅值V_i。
3.将计算得到的电压幅值V_i和电压相角θ_i作为该节点的新的电压值。
4.对于其他节点,计算它们的电压相角θ_i和电压幅值V_i,并将其更新为新的电压值。
5.重复2-4步骤,直到收敛或满足收敛条件。
在每次迭代过程中,我们可以根据收敛准则来判断是否达到收敛,通常是通过计算两次迭代之间电压的变化量来判断。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
%本程序的功能是用牛拉法进行潮流计算%原理介绍详见鞠平著《电气工程》% 默认数据为鞠平著《电气工程》例8.4 所示数据%B1 是支路参数矩阵%第一列和第二列是节点编号。
节点编号由小到大编写 %对于含有变压器的支路,第一列为低压侧节点编号,第二列为高压侧节点编号 %第三列为支路的串列阻抗参数,含变压器支路此值为变压器短路电抗 %第四列为支路的对地导纳参数,含变压器支路此值不代入计算 %第五烈为含变压器支路的变压器的变比,变压器非标准电压比%第六列为变压器是否是否含有变压器的参数,其中“1”为含有变压器,“ 0”为不含有变压器%B2 为节点参数矩阵%第一列为节点注入发电功率参数 %第二列为节点负荷功率参数 %第三列为节点电压参数%第四列%第五列 %第六列为节点类型参数,“ 1”为平衡节点,“2”为 PQ 节点,“3”为 PV 节点参数%X 为节点号和对地参数矩阵%第一列为节点编号%第二列为节点对地参数clear;clc;num=input(' 是否采用默认数据? (1- 默认数据 ;2- 手动输入 )');if num==1n=4;n1=4;isb=4;pr=0.00001;B1=[1 2 0.1667i 0 0.8864 1;1 3 0.1302+0.2479i 0.0258i 1 0;1 4 0.1736+0.3306i0.0344i 1 0;3 4 0.2603+0.4959i 0.0518i 1 0];B2=[0 0 1 0 0 2;0 -0.5-0.3i 1 0 0 2;0.2 0 1.05 0 0 3;0 -0.15-0.1i 1.05 0 0 1];X=[1 0;2 0.05i;3 0;4 0];elsen=input(' 请输入节点数 :n='); n1=input(' 请输入支路数 :n1='); isb=input(' 请输入平衡节点号:isb='); pr=input(' 请输入误差精度:pr='); B1=input(' 请输入支路参数:B1=');B2=input(' 请输入节点参数 :B2='); X=input(' 节点号和对地参数 :X=');endTimes=1; % 迭代次数%创建节点导纳矩阵Y=zeros(n);for i=1:n1if B1(i,6)==0 % 不含变压器的支路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(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else % 含有变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-B1(i,5)/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+B1(i,5)/B1(i,3)+(1-B1(i,5))/B1(i,3);Y(q,q)=Y(q,q)+B1(i,5)/B1(i,3)+(B1(i,5)*(B1(i,5)-1))/B1(i,3);endendfor i=1:n1Y(i,i)=Y(i,i)+X(i,2); % 计及补偿电容电纳enddisp(' 导纳矩阵为: ');disp(Y); % 显示导纳矩阵% 初始化 OrgS 、DetaS OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);% 创建 OrgS ,用于存储初始功率参数h=0;j=0;for i=1:n % 对 PQ 节点的处理if i~=isb&B2(i,6)==2 % 不是平衡点 & 是 PQ 点h=h+1;for j=1:n% 公式 8-74%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej) OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2( j ,3)));OrgS(2*h,1)j,3))- =OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2( imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:n % 对PV节点的处理,注意这时不可再将 h初始化为0if i~=isb&B2(i,6)==3 % 不是平衡点 &是 PV 点h=h+1;for j=1:n% 公式 8-75-a%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej) %Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2( j ,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendend% 创建 PVU 用于存储 PV 节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endend% 创建 DetaS ,用于存储有功功率、无功功率和电压幅值的不平衡量 h=0;for i=1:n % 对 PQ 节点的处理if i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); %delPi DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1); %delQi endendt=0;for i=1:n % 对PV节点的处理,注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1); %delPiDetaS(2*h,1)=real(PVU(t,1))A2+imag(PVU(t,1))A2-real(B2(i,3))A2-imag(B2(i,3))A2; %delUi endend% DetaS%创建I,用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));%conj 求共轭endend% 创建 Jacbi( 雅可比矩阵 )Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n % 对 PQ 节点的处理if B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j % 对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else % 非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1) % 将用于内循环的指针置于初始值,以确保雅可比矩阵换行 k=0;endendendendendk=0;for i=1:n % 对 PV 节点的处理if B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==j % 对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3)); Jacbi(2*h,2*k)=2*real(B2(i,3));else % 非对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) % 将用于内循环的指针置于初始值,以确保雅可比矩阵换行 k=0;endendendendenddisp(' 初始雅可比矩阵为: '); disp(Jacbi);%求解修正方程,获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %inv 矩阵求逆% DetaU%修正节点电压j=0;for i=1:n % 对 PQ 节点处理if B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n % 对 PV 节点的处理if B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endend% B2% 开始循环 ********************************************************************** while abs(max(DetaU))>prOrgS=zeros(2*n-2,1);h=0;j=0;for i=1:nif i~=isb&B2(i,6)==2 h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))- imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2( j ,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2( imag(Y(i,j)) *imag(B2(j,3)))- real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:nif i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))- imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2( j ,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2( imag(Y(i,j)) *imag(B2(j,3)))- real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendend% OrgS% 创建 DetaSh=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,6)==3 h=h+1;t=t+1;% DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))A2+imag(PVU(t,1))A2-real(B2(i,3))A2- imag(B2(i,3)F2;endend% DetaS% 创建 Ii=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endend% I% 创建 JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=- j,3))- j,3))-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1)); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendend% JacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;% DetaU%修正节点电压j=0;for i=1:nif B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end endfor i=1:nif B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endend% B2Times=Times+1; % 迭代次数加 1enddisp(' 迭代次数为: ');disp(Times);disp(' 收敛时电压修正量为:: '); disp(DetaU);for k=1:nE(k)=B2(k,3); e(k)=real(E(k));f(k)=imag(E(k));V(k)=sqrt(e(kF2+f(k)A2);sida(k)=atan(f(k)./e(k))*180./pi;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 注入电流的共轭值 end disp(' 各节点的功率S 为(节点号从小到大排列 ):');disp(S); % 显示各节点的注入功率Sline=zeros(n1,5);disp(' --------------------------------- ');disp(' 各条支路的首端功率 Si 为(顺序同您输入 B1 时一致 ):'); for i=1:n1 p=B1(i,1);q=B1(i,2);Sline(i,1)=B1(i,1);Sline(i,2)=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))*((1-B1(i,5))/B1(i,3))+(conj(E(p))-conj(E(q)))*(B1(i,5)/B1(i,3)));Siz(i)=Si(p,q);endSSi(p,q)=Si(p,q);Sline(i,3)=Siz(i);ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))]; disp(ZF);enddisp(' --------------------------------- ');disp(' 各条支路的末端功率 Sj 为(顺序同您输入 B1 时一致 ): ');for i=1:n1p=B1(i,1);q=B1(i,2);if B1(i,6)==0Sj(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))*((B1(i,5)*(B1(i,5)-1))/B1(i,3))+(conj(E(q))-conj(E(p)))*(B1(i,5)/B1(i,3)));Sjy(i)=Sj(q,p);endSSj(q,p)=Sj(q,p);Sline(i,4)=Sjy(i);ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))]; disp(ZF);enddisp(' --------------------------------- ');disp(' 各条支路的功率损耗 DS 为 (顺序同您输入 B1 时一致 ): '); for i=1:n1 p=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);DDS(i)=DS(i);Sline(i,5)=DS(i); ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);enddisp(' --------------------------------- ');disp(' 各支路首端编号末端编号首端功率线路损耗 ');末端功率disp(Sline);。