基于MATLAB的牛顿拉夫逊迭代法计算潮流(附加短路计算).doc

这个程序可以适用于三机九节点系统（参数见主程序），本来是想编一个通用各种结构的程序的，但是因为鄙人比较懒，老师留作业时候又没说要通用，就没改完。惭愧啊。不过大同小异啦。有兴趣的慢慢改吧。使用方法：按后文中给出的代码建立.m文件放在一个文件夹里面。先运行Powerflow_main.m计算算例系统的潮流；然后运行ShortcircuitCalc.m计算算例系统三相短路电流；程序说明详见各.m文件注释部分，写的已经很详细了，慢慢看吧。Powerflow_main.m文件代码如下：clear %牛顿拉夫逊迭代法计算潮流 format short %规定参数数据显示精度%节点参数矩阵 %第一列为节点编号 %第二列表示有功注入P %第三列表示无功注入Q %第四列表示电压幅值U %第五列表示电压角度 %第六列表示发电机x %第七列表示发电机E %第八列表示节点类型（2表示平衡节点，1表示PV节点，0表示PQ节点） Node_p= 1, 0, 0, 1.04 , 0, 0.3, 1.137, 2; 2, 1.63, 0, 1.025, 0, 0.3, 1.211, 1; 3, 0.85, 0, 1.025, 0, 0.3, 1.047, 1; 4, 0, 0, 1.0, 0, 0, 0, 0; 5, -1.25, -0.5, 1.0, 0, 0, 0, 0; 6, -0.9, -0.3, 1.0, 0, 0, 0, 0; 7, 0, 0, 1.0, 0, 0, 0, 0; 8, -1, -0.35, 1.0, 0, 0, 0, 0; 9, 0, 0, 1.0, 0, 0, 0, 0 ; count_s=0; countPV=0; for k=1:size(Node_p,1) if Node_p(k,8)=1 countPV=countPV+1; else if Node_p(k,8)=2 count_s=count_s+1; end; end; end; countPV; count_s; countPQ=size(Node_p,1)-1-countPV; %显示节点参数 disp(节点参数如下：) disp(Node_p) %支路参数%第一列为首节点，第二列为末节点，第三列表示R，第四列表示X，第五列表示B/2%第六列表示支路类型（1为变比为1的变压器元件；2为输电线元件;0为接地支路）Branch_p = 1, 4, 0 , 0.0576, 0 , 1; 2, 7, 0 , 0.0625, 0 , 1; 3, 9, 0 , 0.0586, 0 , 1; 4, 5, 0.01 , 0.085 , 0.088 , 2; 4, 6, 0.017 , 0.092 , 0.079 , 2; 5, 7, 0.032 , 0.161 , 0.153 , 2; 6, 9, 0.039 , 0.17 , 0.179 , 2; 7, 8, 0.0085, 0.072 , 0.0745, 2; 8, 9, 0.0119, 0.1008, 0.1045, 2 ;%显示支路参数 disp(支路参数如下：) disp(Branch_p)%设置节点初值 U=Node_p(:,4); e_ang=Node_p(:,5); P=Node_p(:,2); Q=Node_p(:,3); save data.mat format long%计算结果数据显示精度%显示节点导纳矩阵admi();disp(节点导纳矩阵Y); sparseY Kmax=10; %设置最大迭代次数k accuracy=10-7;%设置迭代精度 k=0;%迭代次数初始化为零 for k1=1:Kmax dP,dQ,y=getY(U,e_ang); if max(abs(y)accuracy break; end; J=jacob(U,e_ang,dP,dQ); x=-inv(J)*y; de_ang=0;x(1:8); dU=x(9:14); u1=U(2:9)*dU.; u=diag(u1); U(4:9)=U(4:9)+u; e_ang=e_ang+de_ang; k=k+1; end; e_ang=e_ang/pi*180; save result.mat disp(迭代次数：) k disp(4号节点至9号节点电压幅值如下：) disp(U(4:9) disp(2号节点至9号节点电压相角如下：)disp(e_ang(2:9)admi.m文件代码如下：%节点导纳矩阵的形成format long %规定数据格式 NI=size(Branch_p,1); k_t=1; Y=zeros(NI); for m1=1:NI; I=Branch_p(m1,1); J=Branch_p(m1,2); R=Branch_p(m1,3); X=Branch_p(m1,4); b=Branch_p(m1,5); Style=Branch_p(m1,6); if Style=1 %判断为变压器元件 Y(I,I)=Y(I,I)+1/(R+1j*X); Y(J,J)=Y(J,J)+1/(R+1j*X)/k_t/k_t; Y(I,J)=Y(I,J)-1/(R+1j*X)/k_t; Y(J,I)=Y(J,I)-1/(R+1j*X)/k_t; else if Style=0 %判断为母线接地支路元件 Y(I,J)=Y(I,J)+1/(R+1j*X); else %判断为输电线元件 Y(I,I)=Y(I,I)+1j*b+1/(R+1j*X); Y(J,J)=Y(J,J)+1j*b+1/(R+1j*X); Y(I,J)=Y(I,J)-1/(R+1j*X); Y(J,I)=Y(J,I)-1/(R+1j*X); end; end; end; Y; G=real(Y); B=imag(Y); sparseY=sparse(Y); save data.mat; getY.m文件代码如下：（线性方程组常写作AX=Y形式，故此处命名为get_Y）%计算P和Q的函数function dP,dQ,y=getY(U,e_ang)load data P Q G B Node_p countPV count_s;dP=zeros(size(Node_p,1),1); sum1=zeros(size(Node_p,1),1); for i=1:size(Node_p,1) for j=1:size(Node_p,1) sum1(i)=sum1(i)+U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j)+B(i,j)*sin(e_ang(i)-e_ang(j); end; dP(i)=P(i)-U(i)*sum1(i); end; dQ=zeros(size(Node_p,1),1);sum2=zeros(size(Node_p,1),1); for i=1:size(Node_p,1) for j=1:size(Node_p,1) sum2(i)=sum2(i)+U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j)-B(i,j)*cos(e_ang(i)-e_ang(j); end; dQ(i)=Q(i)-U(i)*sum2(i); end; y=dP(count_s+1):9);dQ(countPV+count_s+1):9); %拼接P和Q构成方程 -J*x=y的向量yjacob.m文件的代码如下： %形成雅克比矩阵的函数 function J=jacob(U,e_ang,dP,dQ) load data B G Q P Node_p countPV count_s; size_Y=size(Node_p,1); %H矩阵 H1=zeros(size_Y); for i=1:size_Y for j=1:size_Y if j=i; H1(i,j)= U(i)*U(i)*B(i,j)+Q(i)-dQ(i); else H1(i,j)= -U(i)*U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j)-B(i,j)*cos(e_ang(i)-e_ang(j); end; end end H=H1(2:size_Y,2:size_Y); %N矩阵 N1=zeros(size_Y); for i=1:size_Y for j=1:size_Y if j=i; N1(i,j)= -U(i)*U(i)*G(i,j)-(P(i)-dP(i); else N1(i,j)= -U(i)*U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j)+B(i,j)*sin(e_ang(i)-e_ang(j); end; end; end; N=N1(2:size_Y,(countPV+count_s+1):size_Y); %M矩阵 M1=zeros(size_Y); for i=1:size_Y for j=1:size_Y if j=i; M1(i,j)=U(i)*U(i)*G(i,j)-(P(j)-dP(j); else M1(i,j)=U(i)*U(j)*(G(i,j)*cos(e_ang(i)-e_ang(j)+B(i,j)*sin(e_ang(i)-e_ang(j); end; end; end; M=M1(countPV+count_s+1):size_Y,2:size_Y); %L矩阵 L1=zeros(size_Y); for i=1:size_Y for j=1:size_Y if j=i; L1(i,j)= U(i)*U(i)*B(i,j)-(Q(i)-dQ(i); else L1(i,j)= -U(i)*U(j)*(G(i,j)*sin(e_ang(i)-e_ang(j)-B(i,j)*cos(e_ang(i)-e_ang(j); end; end end L=L1(countPV+count_s+1):size_Y,(countPV+count_s+1):size_Y); J=H N;M L;%拼接构成雅克比矩阵ShortcircuitCalc.m文件的代码如下：clear%三相短路计算format long;%对YN进行修正，形成包括发电机内阻抗和负荷阻抗的节点导纳矩阵re_admi();Z=inv(rY);%计算节点阻抗矩阵disp(4节点发生金属短路)f=4;%短路点为4节点%输出节点阻抗矩阵的短路点所在列disp(节点阻抗矩阵的短路点所在列Z(:,f)=);Z(:,f)%短路电流If计算load result U e_ang;e_angf=e_ang(f);Uf=U(f)*(cos(e_angf/180*pi)+1j*sin(e_angf/180*pi);Zff=Z(f,f);zf=0;If=Uf/(Zff+zf);i_angf=angle(If)*180/pi;If=abs(If);%短路电流计算结果显示disp(短路电流幅值:)Ifdisp(短路电流相角(单位为 )i_angf%短路时各节点电压计算U_k=U;for x1=1:size(Node_p,1) U_k(x1)=U(x1)*(cos(e_ang(x1)*pi/180)+1j*sin(e_ang(x1)*pi/180)-Z(x1,f)*If*(cos(i_angf*pi/180)+1j*sin(i_angf*pi/180);end;%Uk为短路时各节点电压幅值Uk=abs(U_k);%uk_ang为短路时各节点电压相角（单位为 ）uk_ang=angle(U_k)*180/pi;%输出短路时各节点电压disp(短路时1-9节点电压：)disp(幅值：)Ukdisp(相角：(单位为 )uk_ang%计算短路时各支路电流%输出矩阵初始化%第一列为首节点i %第二列为末节点j%第三列为Ii幅值%第四列为Ii相角(单位为 ) %第五列为Ij幅值%第六列为Ij相角（单位为 ）Ik=1 4 1 0 1 0; 2 7 1 0 1 0; 3 9 1 0 1 0; 4 5 1 0 1 0; 4 6 1 0 1 0; 5 7 1 0 1 0; 6 9 1 0 1 0; 7 8 1 0 1 0; 8 9 1 0 1 0;load data Branch_p k_t;%其中变压器变比为k_tNI=size(Branch_p,1);for m1=1:NI; I=Branch_p(m1,1); J=Branch_p(m1,2); R=Branch_p(m1,3); X=Branch_p(m1,4); b=Branch_p(m1,5); Style=Branch_p(m1,6); if Style=1 %判断为非变压器支路 Ik(m1,3)=(U_k(I)-U_k(J)/(R+1j*X)+U_k(I)*1j*b; Ik(m1,4)=angle(Ik(m1,3)*180/pi; Ik(m1,3)=abs(Ik(m1,3); Ik(m1,5)=(U_k(J)-U_k(I)/(R+1j*X)+U_k(J)*1j*b; Ik(m1,6)=angle(Ik(m1,5)*180/pi; Ik(m1,5)=abs(Ik(m1,5); else %否则为变压器支路 Ik(m1,3)=(U_k(I)-U_k(J)/(R+1j*X)/k_t; Ik(m1,4)=angle(Ik(m1,3)*180/pi; Ik(m1,3