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

%本程序的功能是用牛拉法进行潮流计算%原理介绍详见鞠平著《电气工程》%默认数据为鞠平著《电气工程》例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）是一种常用于计算潮流的数值求解方法。

它是基于潮流计算的功率流方程的非线性特性而设计的，通过迭代求解来逼近潮流计算的稳态解。

下面将介绍牛顿拉夫逊法计算潮流的基本步骤。

首先，我们需要明确潮流计算的目标，即确定电力系统中各节点的电压相角和幅值。

这些节点是电力系统中的发电机、负荷和交流输电线路的连接点。

通过潮流计算，我们可以得到各节点的电压相角和幅值，从而分析系统的功率分布、电压稳定性等运行特性。

接下来，我们需要建立电力系统的潮流计算模型。

这个模型中，我们需要考虑发电机的注入功率、负荷的吸收功率、线路的传输损耗等因素。

通过利用功率流方程，我们可以将这些因素表示为电压、功率和导纳之间的方程。

然后，我们需要进行初始化操作。

在进行牛顿拉夫逊法迭代计算之前，我们需要对电力系统的各节点进行初始电压值的设定。

这些初始值可以根据经验或者历史数据来得到，但需要满足物理约束条件，如一致性、电压幅值在合理范围内等。

接下来，我们进入迭代计算的过程。

首先，我们需要对系统的节点进行编号，然后选择某一节点作为基准节点，其他节点相对于基准节点的电压相角进行计算。

然后，我们根据节点注入功率和导纳矩阵的关系，得到节点注入电流。

接着，我们根据节点注入电流和电压相角的关系，计算各节点的电压相角和幅值的改变量。

在计算改变量后，我们需要对节点电压进行更新。

更新后，我们判断系统是否达到收敛条件。

如果满足收敛条件，则停止迭代，得到最终的潮流计算结果；如果不满足收敛条件，则继续进行下一轮迭代计算。

最后，我们对潮流计算结果进行分析和验证。

通过比较计算得到的结果和实际运行数据进行对比，我们可以评估潮流计算的准确性。

同时，我们还可以通过故障分析、电压稳定性评估等手段对电力系统进行优化和改进。

总而言之，牛顿拉夫逊法是一种常用的求解潮流计算问题的方法。

它通过迭代求解潮流计算的功率流方程，逼近潮流计算的稳态解。

摘要电力系统潮流计算是研究电力系统稳态运行的一种重要方法,它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态,包括各母线的电压、线路的功率分布以及功率损耗等等。

潮流计算主要用于电网规划和静态安全分析，它可为扩建电力网络，以达到规划周期内所需要的输电能力提供依据；也可以对预想事故进行模拟和分析，校核预想事故下的电力系统安全性。

本文简单介绍了牛顿-拉夫逊潮流计算的原理、模型与算法，然后用具体的实例，利用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 +jB ij ，带入上式，得到有功功率和无功功率方程 P i =V i ∑V j (G ij cos θij +B ij sin θij )n j=1 (1.3)Q i =V i ∑Vj (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）求得，这样就完成了潮流计算。

自动化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：输出结果当潮流计算收敛时，输出最终的节点电压幅值和相角，以及线路的功率流向等参数。

牛顿拉夫逊法是一种高效、快速且收敛性良好的方法，广泛应用于电力系统潮流计算。

在实际应用中，可能会遇到多次迭代或者收敛性不好的情况，此时可以采用退火技术或其他优化算法进行改进。

此外，牛顿拉夫逊法的计算也可以并行化，利用多核处理器或者分布式计算集群来加速计算过程。

总之，牛顿拉夫逊法是一种重要的潮流计算方法，通过迭代计算逼近非线性方程组的解，可以得到电力系统中各节点的电压幅值和相角，用于分析电力系统的稳定性和安全性。

附表1：计算机计算潮流程序：%本程序的功能是用牛顿——拉夫逊法进行潮流计算% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳% 5、支路的变比；6、支路首端处于K侧为1，1侧为0% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号clear;n=13;%input('请输入节点数：n=');nl=13;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%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); %－－－－－－－修改部分－－－－－－－－－－－－ym=0;SB=100;UB=220;%ym=input('您输入的参数是标么值？（若不是则输入一个不为零的数值）'); if ym~=0%SB=input('请输入功率基准值:SB=');%UB=input('请输入电压基准值:UB=');YB=SB./UB./UB;BB1=B1;BB2=B2;for i=1:nlB1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;enddisp('B1矩阵B1=');disp(B1)for i=1:nB2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;enddisp('B2矩阵B2=');disp(B2)end% % %---------------------------------------------------for i=1:nl %支路数if B1(i,6)==0 %左节点处于低压侧p=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; %对角元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;while IT2~=0IT2=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)%求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;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;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/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;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~=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;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;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、△Qfor k2=k1:N1 % 扩展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); % 非对角元规格化endJ(k,k)=1; % 对角元规格化if k~=3 % 不是第三行%======================================================== ====k4=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);%消去运算endJ(k3,k)=0;endif k==N0break;end%==========================================for k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endelsefor k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endendend%====上面是用线性变换方式将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);EE=E*UB;disp(EE);disp('-----------------------------------------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);VV=V*UB;disp(VV);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);enddisp('各节点的功率S为(节点号从小到大排列)：');disp(S);disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');SS=S*SB;disp(SS);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)*SB;ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);%disp(SSi(p,q));disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nlp=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))*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)*SB;ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);%disp(SSj(q,p));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)*SB;ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);%disp(DDS(i));disp('-----------------------------------------------------');endfigure(1);subplot(2,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(2,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;subplot(2,2,3);bar(real(S));ylabel('节点注入有功');grid on;subplot(2,2,4);bar(Siz);ylabel('支路首端无功');grid on;1.冬季最大运行方式潮流计算结果:计算机运行的B1,B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.88+0.545*i 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.88+0.545*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.2+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：节点号电压值相角值支路标号首端功率末端功率支路功率损耗1 242.0000 0 1-2 148.9391+3.327637i -143-27.45476i 5.93913-24.1271i2 231.0000 -3.0341 1-11 89.1956+37.193i -88.19803-61.4687i 0.997519-24.2757i3 227.4154 -4.4383 11-3 88.19803+61.4687i -88-54.5i 0.19803+6.9687i4 226.3304 -4.9004 1-12 77.8364+36.4186i -77.2404-55.7062i 0.596039-19.2876i5 229.8318 -3.9547 12-4 77.2404+55.7062i -77-47.72i 0.24036+7.9862i6 217.7226 -8.2720 1-5 63.5438+12.2204i -61.8773-32.0321i 1.66656-19.8116i7 234.9326 -3.1047 5-6 77.2597+56.3501i -77-47.72i 0.25974+8.6301i8 226.0991 -6.2429 5-7 15.3825-24.3181i 15.5354+0.274108i 0.152961-24.044i9 231.0000 0.5851 7-8 88.20085+61.55009i -88-54.5i 0.20085+7.0501i10 231.0000 3.4950 1-7 40.806-6.05737i -40.0647-22.0555i 0.741264-28.1128i11 236.1931 -1.3348 7-13 -63.6716-39.7687i 64.0965+17.5832i 0.424934-22.1855i12 237.9346 -0.8892 13-9 -64.0965-17.5832i 64.2+21.0187i 0.10348+3.4355i13 238.1868 -2.2028 1-10 79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果图形：2.冬季最小运行方式潮流计算结果：计算机运行的B1B2矩阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.616+0.3817*i 1 0 0 20 0.539+0.3817*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.1+0.06197*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]电压调整措施：变电所1、4变压器变比：+2.5% 水电厂变压器变比：-2.5%5 234.2241 -2.0051 (12--4) 54.0234+42.2706i -53.9-38.17i 0.12342+4.1006i6 226.2159 -4.8577 (1--5) 34.1669+4.22856 -33.6427-28.276i 0.524165-24.0474i7 235.1128 -0.4737 (5--6) 54.0179+37.3169i -53.9-33.4i 0.11789+3.9169i8 223.9941 -2.4603 (5--7) -20.3752-9.04094i 20.5371-15.4558i 0.161947-24.4968i9 231.0000 3.4429 (1--7) 11.0013+1.45186i -10.8258-31.4513i 0.175442-29.9995i10 231.0000 3.9321 (7--8) 53.9768+36.0957i -53.9-33.4i 0.076797+2.6957i11 238.4187 -0.9561 (7--13) -63.6881+10.8115i 64.0874-32.7763i 0.39932-21.9648i12 239.1742 -0.6185 (13--9) -64.0874+32.7763i 64.2-29.0386i 0.11258+3.7377i13 234.9468 0.6942 (1--10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机运行结果的图形：3.夏季最大运行方式计算机计算结果：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.776+0.481*i 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.776+0.481*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：10 231.0000 3.4950 (7,8) 77.7585+53.6634i -77.6-48.1i0.1585+5.5634i11 237.0938 -1.1858 (7,13) -113.4984+9.476406i 114.6926-28.38885i 1.19422-18.9124i12 239.1742 -0.6185 (13,9) -114.6926+28.38885 115-18.18369i0.307384+10.2052i13 233.3912 3.2303 (1,10) -79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果如图：4.夏季最小运行方式：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.543+0.336*i 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.543+0.336*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]10 231.0000 3.9321 (7,8) 54.3735+36.2509i -54.3-33.67i 0.073528+2.5809i11 239.0099 -0.8518 (7,13) -113.4428+24.39707i 114.6754-43.79301i 1.23255-19.3959i12 239.8726 -0.5689 (13,9) -114.6754+43.79301i 115-33.01613i 0.324605+10.7769i13 235.9565 4.4510 (1,10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机计算结果如图：5.夏季故障运行状态:调压及无功补偿措施如下:变电所3的变压器变比为-2.5%,无功补偿容量为20Mvar。

系统的潮流及三相短路电流计算班级：电气1班学号：94姓名：杨鹏摘要潮流计算，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。

是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。

通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。

待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。

牛顿-拉夫逊法是电力系统潮流计算的常用算法之一，收敛性好，迭代次数少。

本文基于MATLAB 计算系统的潮流及三相短路电流。

关键词：潮流计算 matlab 牛顿-拉夫逊1电力系统的潮流计算电力系统常规潮流计算的任务是根据给定电网机构、发电计划及负荷分布情况，求出整个电网的运行状态，其中包括各节点母线电压、相角、线路传输的有功功率和无功功率等。

在电网的潮流计算中，一般给定的运行参数有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。

待求的参数包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。

节点的功率方程对n节点电力系统，节点i注入的有功功率Si：极坐标形式的节点功率方程：直角坐标形式的节点功率方程：节点分类:根据实际运行条件，节点可分成三类：PQ节点、PV节点和平衡节点PQ节点：节点注入的有功P和无功Q皆为给定量的节点。

一般负荷节点，联络节点和给定有功和无功的发电机节点在潮流计算中都视作PQ节点，PQ节点的节点电压（其幅值U和相角θ，或其实部e和虚部f）为待求变量。

PV节点：节点注入的有功P和无功Q皆为给定量的节点。

一般负荷节点，联络节点和给定有功和无功的发电机节点在潮流计算中都视作PQ节点，PQ节点的节点电压（其幅值U和相角θ，或其实部e和虚部f）为待求变量。

平衡节点：平衡节点的节点电压是给定值，对极坐标形式的节点功率方程，平衡节点的电压幅值一般情况下取作U=，相角取作00.θ=，对直角坐标形式的节点功率方程，平衡节点的实部和虚部一般分别取作e=和 f=。

I力系统分析课稈设计设计题目9节点电力网络潮流计算指导教师_______________________________________ 院（系、部）电气与控制工程学院_________ 专业班级_______________________________________ 学号 __________________________________ 姓名 __________________________________ 日期 ____________________________________电气工程系课程设计标准评分模板课程设计成绩评定表不及格标准：设计内容一项否决制，即5为不及格，整个设计不及格，其他4项否决；优、良、中、及格标准：以设计内容为主体，其他项超过三分之一为评定标准，否则评定为下一等级；如优秀评定，设计内容要符合5,其余九项要有4项符合才能评定为优，否则评定为良好,以此类推。

评定教师签字:最终成绩:目录1PSAS嗽件简介 (1)1.1PSASP平台的主要功能和特点 (1)1.2PSASP勺平台组成 (2)2牛顿拉夫逊潮流计算简介 (3)2.1牛顿—拉夫逊法概要 (3)2.2直角坐标下的牛顿—拉夫逊潮流计算 (5)2.3牛顿—拉夫逊潮流计算的方法 (6)3九节点系统单线图及元件数据 (7)3.1九节点系统单线图 (7)3.2系统各项元件的数据 (8)4潮流计算的结果 (10)4.1潮流计算后的单线图 (10)4.2潮流计算结果输出表格 (10)5结论 (14)电力系统分析课程设计任务书9节点系统单线图如下：基本数据如下：表1母线数据表2交流线数据表3两绕组变压器数据两绕组变压器数据（续表）表4发电机数据负荷数据表5区域数据1 PSASP 软件简介“电力系统分析综合程序” （Power System Analysis Software Package,PSA）SP 是一套历史悠久、功能强大、使用方便的电力系统分析程序，是高度集成和开发具有我国自主知识产权的大型软件包。

牛顿—拉夫逊法潮流计算一、 潮流计算的基本原理实际电力系统中的节点类型5二、实际电力系统中的节点类型123452s 3s 4s 过渡节点：PQ 为0的给定PQ 节点，如图的节点5网络中各节点的性质：负荷节点：给定功率P 、Q 如图中的3、4节点如图中的节点1，可能有两种情况：给定P 、Q 运行，给定P 、V 运行3. 负荷发电机混合节点：PQ 节点，如图中的节点2发电机节点负荷节点负荷节点混合节点过渡节点1. 负荷节点：2. 发电机节点：4.潮流计算中节点类型划分6三、潮流计算中节点类型的划分也称为松弛节点，摇摆节点123452s 3s 4s 平衡节点PQ 节点PQ 节点PV 节点PQ 节点PQ∈Ω1. PQ 节点：已知P 、Q负荷、过渡节点，PQ 给定的发电机节点，大部分节点PV ∈Ω给定PV 的发电机节点，具有可调电源的变电所，少量节点2.PV 节点：已知P 、V3. 平衡节点＋基准节点：已知V 、δ采用极坐标，节点电压表示为()cos sin i i i i i i V V V j δδδ=∠=+节点功率将写成⎪⎪⎭⎪⎪⎬⎫-=+=∑∑==n j ij ij ij ij j i i nj ij ij ij ij j i i B G V V Q B G V V P 11)cos sin ()sin cos (δδδδ (1) 式中，ij i j δδδ=-，是i 、j 两节点电压的相角差。

方程式把节点功率表示为节点电压的幅值和相角的函数。

在有n 个节点的系统中，假定第1~m 号节点为P Q 节点，第1~1m n +-号节点为PV 节点，第n 号节点为平衡节点。

n V 和n δ是给定的，PV 节点的电压幅值11~m n V V +-也是给定的。

因此，只剩下1n -个节点的电压相角121,,,n δδδ- 和m 个节点的电压幅值12,,,m V V V 是未知量。

实际上，对于每一个P Q 节点或每一个PV 节点都可以列写一个有功功率不平衡量方程式()1(cos sin )01,2,,1ni is i is i j ij ij ij ij j P P P P V V G B i n δδ=∆=-=-+==-∑ (2)而对于每一个P Q 节点还可以再列写一个无功功率不平衡量方程式()1(sin cos )01,2,,ni is i is i j ij ij ij ij j Q Q Q Q V V G B i m δδ=∆=-=--==∑ (3)式(2)和式(3)一共包含了1n m -+个方程式，正好同未知量的数目相等，而比直角坐标形式的方程少了1n m -+个。

摘要本文，首先简单介绍了基于在MALAB中行潮流计算的原理、意义，然后用具体的实例，简单介绍了如何利用MALAB去进行电力系统中的潮流计算。

众所周知，电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各线的电压、各元件中流过的功率、系统的功率损耗等等。

在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性和经济性。

此外，在进行电力系统静态及暂态稳定计算时，要利用潮流计算的结果作为其计算的基础；一些故障分析以及优化计算也需要有相应的潮流计算作配合；潮流计算往往成为上述计算程序的一个重要组成部分。

以上这些，主要是在系统规划设计及运行方式安排中的应用，属于离线计算范畴。

牛顿－拉夫逊法在电力系统潮流计算的常用算法之一，它收敛性好，迭代次数少。

本文介绍了电力系统潮流计算机辅助分析的基本知识及潮流计算牛顿－拉夫逊法，最后介绍了利用MTALAB程序运行的结果。

关键词：电力系统潮流计算，牛顿－拉夫逊法，MATLABABSTRACTThis article first introduces the flow calculation based on the principle of MALAB Bank of China, meaning, and then use specific examples,a brief introduction, how to use MALAB to the flow calculation in power systems.As we all know, is the study of power flow calculation of power system steady-state operation of a calculation, which according to the given operating conditions and system wiring the entire power system to determine the operational status of each part: the bus voltage flowing through the components power, system power loss and so on. In power system planning power system design and operation mode of the current study, are required to quantitatively calculated using the trend analysis and comparison of the program or run mode power supply reasonable, reliability and economy. In addition, during the power system static and transient stability calculation, the results of calculation to take advantage of the trend as its basis of calculation; number of fault analysis and optimization also requires a corresponding flow calculation for cooperation; power flow calculation program often become the an important part. These, mainly in the way of system design and operation arrangements in the application areas are off-line calculation.Newton - Raphson power flow calculation in power system is one commonly used method, it is good convergence of the iteration number of small, introduce the trend of computer-aided power system analysis of the basic knowledge and power flow Newton - Raphson method, introduced by the last matlab run results.Keywords：power system flow calculation, Newton – Raphson method, matlab目录1 绪论 (1)1.1 课题背景 (1)1.2 电力系统潮流计算的意义 (2)1.3 电力系统潮流计算的发展 (2)1.4 潮流计算的发展趋势 (4)2 潮流计算的数学模型 (5)2.1 电力线路的数学模型及其应用 (5)2.2 等值双绕组变压器模型及其应用 (6)2.3 电力网络的数学模型 (8)2.4 节点导纳矩阵 (9)2.4.1 节点导纳矩阵的形成 (9)2.4.2 节点导纳矩阵的修改 (10)2.5 潮流计算节点的类型 (11)2.6 节点功率方程 (12)7 潮流计算的约束条件 (13)2·3 牛顿－拉夫逊法潮流计算基本原理 (14)3.1 牛顿-拉夫逊法的基本原理 (14)3.2 牛顿-拉夫逊法潮流计算的修正方程 (17)3.3 潮流计算的基本特点 (20)3.4 节点功率方程 (21)4牛顿－拉夫逊法分解潮流程序 (22)41 牛顿－拉夫逊法分解潮流程序原理总框图 (22)·4.2 形成节点导纳矩阵程序框图及代码 (23)4.21 形成节点导纳矩阵程序框图 (23)。

牛拉法潮流计算%本程序的功能是用牛拉法进行潮流计算 %原理介绍详见鞠平著《电气工程》%默认数据为鞠平著《电气工程》例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;感谢您的阅读，祝您生活愉快。

牛顿拉夫逊潮流计算[整理版]牛顿拉夫逊潮流计算// 整个程序为: // 牛拉法解潮流程序#include<stdio.h> #include<math.h> #define N 4 // 节点数#define n_PQ 2 //PQ 节点数#define n_PV 1 //PV 节点数#define n_br 5 // 串联支路数void main(){void disp_matrix(float *disp_p,int disp_m,int disp_n); // 矩阵显示函float Us[2*N]={1.0,0,1.0,0,1.05,0,1.05,0}; // float Ps[N]={0,-0.5,0.2}; //float G[N][N],B[N][N]; //struct// 阻抗参数{int nl; // 左节点int nr; // 右节点float R; // 串联电阻值float X; // 串联电抗值float Bl; // 左节点并联电导float Br; // 右节点并联电纳电压初值有功初值float Qs[N]={0,-0.3}; // 无功初值各几点电导电纳}ydata[n_br]={{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},{2,2,0,0.05,0,0}};float Z2; //Z A2=R A2+X A2 各串联阻抗值的平方float e[N],f[N],dfe[2*(N-1)]; 〃e,f 存储电压的x 轴分量和y 轴分量,dfe 存储电压修正值float mid1[N],mid2[N],dS[2*(N-1)]; //midi 、mid2存储计算雅克比行列式对角线元素的中间值,dS 存储PQU勺不平衡量float Jacob[2*(N-1)][2*(N-1)],inv_J[2*(N-1)][2*(N-1)];// 雅克比行列式float dPQU=1.0; //PQU 不平衡量最大值int kk=0; // 迭代次数int i,j,k;float t;float Pij[n_br]; // 存储线路i->j 的有功的无功float Qij[n_br]; // 存储线路i->j的有功float Pji[n_br]; // 存储线路j->i的无功float Qji[n_br]; // 存储线路j->i的有功损耗float dPij[n_br]; // 存储线路i->jfloat dQij[n_br]; // 存储线路i->j 的无功损耗float AA,BB,CC,DD; // 存储线路潮流计算时的中间值// 形成导纳矩阵for(i=0;i<N;i++)for(j=0;j<N;j++){G[i][j]=0;B[i][j]=0;}for(i=0;i<n_br;i++){if(ydata[i].nl!=ydata[i].nr){Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X);// 串联阻抗等效导纳值// 非对角元素G[ydata[i].nl-1][ydata[i].nr-1]=(-ydata[i].R)/Z2;B[ydata[i].nl-1][ydata[i].nr-1]=ydata[i].X/Z2;G[ydata[i].nr-1][ydata[i].nl-1]=(-ydata[i].R)/Z2;B[ydata[i].nr-1][ydata[i].nl-1]=ydata[i].X/Z2;// 对角元素G[ydata[i].nl-1][ydata[i].nl-1]+=ydata[i].R/Z2;G[ydata[i].nr-1][ydata[i].nr-1]+=ydata[i].R/Z2;B[ydata[i].nl-1][ydata[i].nl-1]+=(-ydata[i].X/Z2);B[ydata[i].nr-1][ydata[i].nr-1]+=(-ydata[i].X/Z2);// 并联导纳等效导纳值B[ydata[i].nl-1][ydata[i].nl-1]+=ydata[i].Bl;B[ydata[i].nr-1][ydata[i].nr-1]+=ydata[i].Br;}else{G[ydata[i].nl-1][ydata[i].nr-1]+=ydata[i].R;B[ydata[i].nl-1][ydata[i].nr-1]+=ydata[i].X;}}printf("G=\n");disp_matrix(*G,N,N);printf("B=\n");disp_matrix(*B,N,N);// 分离e,ffor(i=0;i<N;i++){ e[i]=Us[2*i]; f[i]=Us[2*i+1];}// 主程序========================================== while(dPQU>0.00001) { // 计算功率不平衡量for(i=0;i<N-1;i++){mid1[i]=0;mid2[i]=0;for(j=0;j<N;j++){mid1[i]=mid1[i]+G[i][j]*e[j]-B[i][j]*f[j]; mid2[i]=mid2[i]+G[i][j]*f[j]+B[i][j]*e[j];} dS[2*i]=Ps[i]-(e[i]*mid1[i]+f[i]*mid2[i]);if(i<n_PQ) dS[2*i+1]=Qs[i]-(f[i]*mid1[i]-e[i]*mid2[i]);else dS[2*i+1]=Us[2*i]*Us[2*i]-(e[i]*e[i]+f[i]*f[i]);dPQU=0;for(i=0;i<2*(N-1);i++){if(dS[i]<0&&dPQU<-dS[i])dPQU=-dS[i];else if(dS[i]>0&&dPQU<dS[i])dPQU=dS[i];}if(dPQU>0.00001){kk++;// 形成雅克比行列式for(i=0;i<2*(N-1);i++) for(j=0;j<2*(N-1);j++) Jacob[i][j]=0; for(j=0;j<N-1;j++) { // 求H,Nfor(i=0;i<N-1;i++){if(i!=j)Jacob[2*i][2*j]=B[i][j]*e[i]-G[i][j]*f[i]; Jacob[2*i][2*j+1]=-G[i][j]*e[i]-B[i][j]*f[i];}else{Jacob[2*i][2*i]=B[i][i]*e[i]-G[i][i]*f[i]-mid2[i];Jacob[2*i][2*i+1]=-G[i][j]*e[i]-B[i][j]*f[i]-mid1[i]; }}// 求J,L for(i=0;i<n_PQ;i++){if(i!=j){Jacob[2*i+1][2*j]=G[i][j]*e[i]+B[i][j]*f[i]; Jacob[2*i+1][2*j+1]=B[i][j]*e[i]-G[i][j]*f[i]; }else{Jacob[2*i+1][2*i]=G[i][j]*e[i]+B[i][j]*f[i]-mid1[i];Jacob[2*i+1][2*i+1]=B[i][j]*e[i]-G[i][j]*f[i]+mid2[i];// 求R,Sfor(i=n_PQ;i<N-1;i++){if(i==j){Jacob[2*i+1][2*i]=-2*f[i];Jacob[2*i+1][2*i+1]=-2*e[i]; }}}// 雅克比行列式求逆// 初始化inv_J[N][N] for(i=0;i<2*(N-1);i++) for(j=0;j<2*(N-1);j++) {if(i!=j)inv_J[i][j]=0;elseinv_J[i][j]=1;}// 将原矩阵化简为对角阵for(i=0;i<2*(N-1);i++)for(j=0;j<2*(N-1);j++){if(i!=j){t=Jacob[j][i]/Jacob[i][i];for(k=0;k<2*(N-1);k++){Jacob[j][k]-=Jacob[i][k]*t; inv_J[j][k]-=inv_J[i][k]*t;}}}}// 原矩阵各对角元素化为1，画出逆矩阵for(i=0;i<2*(N-1);i++) if(Jacob[i][i]!=1){ t=Jacob[i][i];for(j=0;j<2*(N-1);j++)inv_J[i][j]=inv_J[i][j]/t;}// 求电压修正值for(i=0;i<2*(N-1);i++){dfe[i]=0;for(j=0;j<2*(N-1);j++) dfe[i]-=inv_J[i][j]*dS[j];}for(i=0;i<N-1;i++){e[i]+=dfe[2*i+1];f[i]+=dfe[2*i];}}elsebreak;}// 循环结束// 求平衡节点功率mid1[N-1]=0; mid2[N-1]=0;fo 「〒0 八A N j ++) 宀 mid 二 N ±ll m i d =N i +G 〔N ±M r e u 〒s N i u r f s m i d 2z ±ll m i d 2 运厶+GLN 厶曰*S +S N 厶曰*e s p s z t H e z 二 rmidlLNt+HNtrmidpN —二 Qs 运二ll H N t 「mid 1运二〒竺N t rmidpN —二 forTn —PQxN —lp 「inff(=kkH%2n=kk)八p 「inff(=PH=)八f o 「T o x N T +) p 「inff(=%9.4『「ps三)八 P「inff(-VIQH=)八f o 「T o x N T+)prinff(=%9.4『「QS§-p 「inff(-vleH=)八 f o 「Tox N T +)「printf("%9.4f",f[i]);printf("\n");// 求线路上的潮流// 计算S[i][j]for(i=0;i<n_br;i++){if(ydata[i].nl!=ydata[i].nr){Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X);AA=-f[ydata[i].nl-1]*ydata[i].Bl+(e[ydata[i].nl-1]-e[ydata[i].nr- 1])*ydata[i].R/Z2+(f[ydata[i].nl-1]-f[ydata[i].nr-1])*ydata[i].X/Z2;BB=-e[ydata[i].nl-1]*ydata[i].Bl-(f[ydata[i].nl-1]-f[ydata[i].nr- 1])*ydata[i].R/Z2+(e[ydata[i].nl-1]-e[ydata[i].nr-1])*ydata[i].X/Z2;Pij[i]=e[ydata[i].nl-1]*AA-f[ydata[i].nl-1]*BB;Qij[i]=e[ydata[i].nl-1]*BB+f[ydata[i].nl-1]*AA;printf("S[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nl,ydata[i].nr,Pij[i],Qij [i]);dPij[i]=Pij[i]+Pji[i];dQij[i]=Qij[i]+Qji[i];}}printf("\n");// 计算S[j][i]for(i=0;i<n_br;i++)if(ydata[i].nl!=ydata[i].nr){Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X);CC=-f[ydata[i].nr-1]*ydata[i].Br+(e[ydata[i].nr-1]-e[ydata[i].nl-1])*ydata[i].R/Z2+(f[ydata[i].nr-1]-f[ydata[i].nl-1])*ydata[i].X/Z2;DD=-e[ydata[i].nr-1]*ydata[i].Br-(f[ydata[i].nr-1]-f[ydata[i].nl-1])*ydata[i].R/Z2+(e[ydata[i].nr-1]-e[ydata[i].nl-1])*ydata[i].X/Z2; Pji[i]=e[ydata[i].nr-1]*CC-f[ydata[i].nr-1]*DD;Qji[i]=e[ydata[i].nr-1]*DD+f[ydata[i].nr-1]*CC;printf("S[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nr,ydata[i].nl,Pji[i],Qji[i]);}}printf("\n");// 计算dS[i][j]for(i=0;i<n_br;i++){if(ydata[i].nl!=ydata[i].nr){dPij[i]=Pij[i]+Pji[i];dQij[i]=Qij[i]+Qji[i];printf("dS[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nl,ydata[i].nr, dPij[i],dQij[i]);} printf("\n");}// 主程序结束// 矩阵显示函数============================================= void disp_matrix(float *disp_p,int disp_m,int disp_n){int i,j;for(i=0;i<disp_m;i++){ for(j=0;j<disp_n;j++)printf("%9.4f",*(disp_p+i*disp_m+j));printf("\n");}printf("\n");}。

clc;clear;close all;n1=9;n2=9;B1=load('b1.txt','ascii');B2=load('b2.txt','ascii');m=0;r=0;for i=1:n1if B1(i,8)==1m=m+1;endendr=n1-m-1;e=B1(1:n1,2);f=B1(1:n1,3);ps=B1(1:n1,4)-B1(1:n1,6);qs=B1(1:n1,5)-B1(1:n1,7);Y=zeros(n1);for j=1:n2p=B2(j,2);q=B2(j,3);J=zeros(2*n1);Y(p,q)=-1./((B2(j,4)+1i*B2(j,5))*B2(j,7));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B2(j,4)+1i*B2(j,5))+1i*B2(j,6);Y(q,q)=Y(q,q)+1/((B2(j,4)+1i*B2(j,5))*B2(j,7)^2)+1i*B2(j,6);end%导纳矩阵YH1=zeros(n1);N1=zeros(n1);J1=zeros(n1);L1=zeros(n1);G=real(Y);B=imag(Y);dp(1:(n1-1))=1;dq(1:m)=1;dx(1:(2*m+r))=1;while(abs(max(dp)>=1.0e-5)||abs(max(dq)>=1.0e-5)||abs(max(dx(1:(m+r)))>=1.0e-5)||abs(max(dx((m +r+1):(2*m+r)).*e(1:m))>=1.0e-5))for i=1:n1for j=1:n1si(i,j)=sin(f(i)-f(j));co(i,j)=cos(f(i)-f(j));if i==jH1(i,j)=e(i)^2*B(i,i)+qs(i);N1(i,j)=-e(i)^2*G(i,i)-ps(i);J1(i,j)=e(i)^2*G(i,i)-ps(i);L1(i,j)=e(i)^2*B(i,i)-qs(i);elseH1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));N1(i,j)=-e(i)*e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));J1(i,j)=e(i)*e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));L1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));endendendH=H1(2:n1,2:n1);N=N1(2:n1,r+2:n1);J=J1(r+2:n1,2:n1);L=L1(r+2:n1,r+2:n1);J=[H N;J L];%雅克比矩阵Jfor i=2:n1a(i)=0;b(i)=0;for j=1:n1a(i)=a(i)+e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));if i>=(r+2)b(i)=b(i)+e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));endenddp(i-1)=ps(i)-e(i)*a(i);if i>=(r+2)dq(i-(r+1))=qs(i)-e(i)*b(i);endendF=[dp';dq'];dx=-J\F;f(2:n1)=f(2:n1)+dx(1:(m+r));e1=e((r+2):n1)+dx((m+r+1):(2*m+r)).*e((r+2):n1); e2=e(1:(r+1));e=[e2;e1];enddisp('各节点电压');disp('e=');disp(e');disp('单位：度f=');disp(f'*180/3.1415926)%求出各节点电压幅值及相角a=0;for i=1:n1a=a+conj(Y(1,i))*(e(i)*cos(f(i))-1i*e(i)*sin(f(i))); enda=a*(e(1)*cos(f(1))+1i*e(1)*sin(f(1)));disp('平衡节点P+jQ=');disp(a);disp('输电线路有功无功:');for i=1:n2p=B2(i,2);q=B2(i,3);s1(p,q)=e(p)^2*(-1i*B2(i,6))+(e(p)*cos(f(p))+1i*e(p)*sin(f(p)))*(e(p)*cos(f(p))-1i*e(p)*sin(f(p))-e( q)*cos(f(q))+1i*e(q)*sin(f(q)))*conj(-Y(q,p));s1(q,p)=e(q)^2*(-1i*B2(i,6))+(e(q)*cos(f(q))+1i*e(q)*sin(f(q)))*(e(q)*cos(f(q))-1i*e(q)*sin(f(q))-e( p)*cos(f(p))+1i*e(p)*sin(f(p)))*conj(-Y(q,p));fprintf('S%d%d=',p,q);disp(s1(p,q));fprintf('S%d%d=',q,p);disp(s1(q,p));end%初值计算for j=1:n1s(j)=0;for i=1:n1s(j)=s(j)+conj(Y(j,i))*(e(i)*cos(f(i))-1i*e(i)*sin(f(i)));ends(j)=s(j)*(e(j)*cos(f(j))+1i*e(j)*sin(f(j)));endY1=conj(-s)./(e'.^2);%负荷等值导纳disp('各节点等值导纳');disp('Y1=');disp(Y1);disp('各节点功率');disp('S=');disp(s);Ra=[0 0 0];X_d=[0.0608 0.1198 0.1813];Xq=[0.0969 0.8645 1.2578];%Xq=X_d;%不计凸极效应v=e.*cos(f)+1i*e.*sin(f);for i=1:3I(i)=conj(s(i)/(e(i)*cos(f(i))+1i*e(i)*sin(f(i))));EQ(i)=e(i)*cos(f(i))+1i*e(i)*sin(f(i))+(Ra(i)+1i*Xq(i))*I(i);EQx(i)=real(EQ(i));EQy(i)=imag(EQ(i));delta(i)=atan(EQy(i)/EQx(i));enddisp('各发电机的暂态电动势，功角和输入机械功率初值');disp('delta=');disp(delta*180/3.1415926);for i=1:3Vx(i)=e(i)*cos(f(i));Vy(i)=e(i)*sin(f(i));Vdq=[sin(delta(i)) -cos(delta(i));cos(delta(i)) sin(delta(i))]*[Vx(i);Vy(i)];Vd(i)=Vdq(1);Vq(i)=Vdq(2);Ix(i)=real(I(i));Iy(i)=imag(I(i));Idq=[sin(delta(i)) -cos(delta(i));cos(delta(i)) sin(delta(i))]*[Ix(i);Iy(i)];Id(i)=Idq(1);Iq(i)=Idq(2);E_q(i)=Vq(i)+Ra(i)*Iq(i)+X_d(i)*Id(i);enddisp('E_q=');disp(E_q);for i=1:3Pm(i)=real(s(i))+(Ix(i)^2+Iy(i)^2)*Ra(i); enddisp('Pm=');disp(Pm);。