电力系统潮流计算牛顿拉夫迅法与PQ分解法通用MATLAB计算程序

1、此程序经40余同学使用检验，无误。 这是一个电气狗熬两个礼拜图书馆的成果，根据华中科技大学电力系统分析中原理编写，可用牛 顿-拉夫逊和 PQ分解法计算给定标幺值条件的潮流。本人水平有限，仅供参考，欢迎一起找 Bug。 2018/07/06 模型计算。 说明：由于本人变压器建模与 PSASP不同，本人使用模型如下图，参数输入时请按该 2018/06/18 主程序 主程序更新：增加补偿电容参数 % file n ame:chaoliu_lj.m % auther:山东科技大学罗江 % function:使用牛顿-拉夫逊法、PQ分解法计算潮流 % updata:2018/6Z18 13:22增加补偿

2、电容参数 标号 1 2 3 %节点类型 %PQ节点 %PV节点 %slack节点 %能计算给定标幺值网络，有且仅有一个平衡节点的潮流 %注意：母线标号顺序要求：PQ节点-PV节点-平衡节点 %若某元件不存在，其导纳为0，阻抗为inf clear %清除工作空间变量 cic %清屏 % % %数据输入（标幺值） SB=100; %基准容量，单位 MVA %母线基准电压 Bus=115 10.5 115 115; %交流线参数：I侧母线J侧母线阻抗1/2接地导纳 Lin e=4 1 0.06125+0.09527i 0; 4 3 0.08469+0.12703i 0; 1 3 0.13989+0.

3、15501i 0; %变压器参数：I侧母线J侧母线 阻抗 变比变压器阻抗归算到I侧 Tran s=2 3 0.0137+0.2881i 0.9504; %加接地电容器补偿：母线导纳 Cap=2 0.5i; %发电机参数：母线节点类型P V/U 0 Gen=4 3 1 0; %负荷参数：母线节点类型 P Q %按参考方向，发电机发出功率(正值)，负荷消耗功率(负值) Load=1 1 -0.18 -0.06; 2 1 -0.32 -0.12; mode=1; %1-极坐标下牛拉法，2-PQ分解法 Tmax=50; %最大迭代次数 limit=1.0e-4; %要求精度 % % %变压器n型等效阻

4、抗参数 Zt=zeros(size(Trans,1),3); Zt(:,1)=Trans(:,3)./Trans(:,4); Zt(:,2)=Trans(:,3)./(1-Trans(:,4); Zt(:,3)=Tra ns(:,3)./(Tra ns(:,4)A2-Tra ns(:,4); Trans_pi=Trans(:,1:2) Zt(:,1) 1./Zt(:,2) 1./Zt(:,3); n=numel(Bus); %总节点数 m=n-1; %PQ 节点数 for i=1:size(Ge n,1)% 数组行数 if Gen(i,2)=2 %除去PV节点就是PQ节点 m=m-1; end

5、 end for i=1:size(Load,1) if Load(i,2)=2 m=m-1; end end %PQ节点包含浮游节点，其 P Q=0 %提取 P,Q,U 向量 P=zeros(1,n); %PQ为原始数据，Pi,Qi为计算结果 Q=zeros(1,n); U=ones(1,n); %电压初始值由此确定 cita=zeros(1,n); %此处未知节点皆设为1.0 / 0 %注意：此处角度单位为度，提取后再转换成 弧度，后面计算使用弧度 for i=1:size(Gen,1) if Gen(i,2)=1 %PQ 节点 P(Gen(i,1)=Gen(i,3); Q(Gen(i,1

6、)=Gen(i,4); end if Gen (i,2)=2 %PV 节点 P(Gen(i,1)=Gen(i,3); U(Gen(i,1)=Gen(i,4); end if Gen( i,2)=3 %slack 节点 U(Gen(i,1)=Gen(i,3); cita(Gen(i,1)=Gen(i,4); end end for i=1:size(Load,1) if Load(i,2)=1 %PQ 节点 P(Load(i,1)=Load(i,3); Q(Load(i,1)=Load(i,4); end if Load(i,2)=2 %PV 节点 P(Load(i,1)=Load(i,3);

7、U(Load(i,1)=Load(i,4); end if Load(i,2)=3 %slack 节点 U(Load(i,1)=Load(i,3); cita(Load(i,1)=Load(i,4); end end disp(初始条件：) dis p(各节点有功：) disp(P); dis p(各节点无功：) disp(Q); ) disp (各节点电压幅值： disp(U); cita=(deg2rad(cita); %角度转换成弧度 disp(各节点电压相角(度)：) disp(rad2deg(cita); %显示依然使用角度 %节点导纳矩阵的计算 Y=zeros(n); %新建节点导

8、纳矩阵 y=zeros(n); %网络中的真实导纳 %计算 y(i,j) for i=1:size(Line,1) %与交流线联结的真实导纳 ii=Line(i,1); jj=Line(i,2); y(ii,jj)=1/Line(i,3); y(jj,ii)=y(ii,jj); end for i=1:size(Trans_pi,1) %与变压器联结的真实导纳 ii=Trans_pi(i,1); jj=Trans_pi(i,2); y(ii,jj)=1/Trans_pi(i,3); y(jj,ii)=y(ii,jj); end %计算 y(i,i) for i=1:size(Line,1) %

9、与交流线联结的对地导纳 ii=Line(i,1); jj=Line(i,2); y(ii,ii)=y(ii,ii)+Line(i,4); y(jj,jj)=y(jj,jj)+Line(i,4); end for i=1:size(Trans_pi,1) %与变压器联结的对地导纳 ii=Trans_pi(i,1); jj=Trans_pi(i,2); y(ii,ii)=y(ii,ii)+Trans_pi(i,4); y(jj,jj)=y(jj,jj)+Trans_pi(i,5); end %算上补偿电容 for i=1:size(Cap,1) ii=Cap(i,1); y(ii,ii)=y(ii

10、,ii)+Cap(i,2); end %由 y 计算 Y ysum=sum(y,1); %每一行求和 for i=1:n for j=1:n if i=j Y(i,j)=ysum(i); else Y(i,j)=-y(i,j); end end end disp( 节点导纳矩阵： ); disp(Y); G=real(Y); %电导矩阵 B=imag(Y); %电纳矩阵 % % %以上为基础数据整理 %接下来是牛拉法的大循环 % % %计算功率不平衡量 dP,dQ,Pi,Qi=Unbalanced( n,m,P,Q,U,G,B,cita ); disp( 有功不平衡量： ); disp(dP)

11、; disp( 无功不平衡量： ); disp(dQ); for i=1:Tmax fprintfC 第%d 次迭代 n,i); %雅可比矩阵的计算 if(mode=1) J=Jacobi( n,m,U,cita,B,G,Pi,Qi ); disp(雅可比矩阵); disp(J); end %求解节点电压修正量 if(mode=1) dU,dcita=Correct( n,m,U,dP,dQ,J ); else dU,dcita=PQ_LJ( n,m,dP,dQ,U,B ); end dis p(电压、相角修正量：)； disp(dU); disp(rad2deg(dcita); %修正节点电

12、压 U=U+dU; cita=cita+dcita; dis p(节点电压幅值：); disp(U); dis p(节点电压相角：); disp(rad2deg(cita); %计算功率不平衡量 dP,dQ,Pi,Qi=Unbalanced( n,m,P,Q,U,G,B,cita ); disp( 有功不平衡量： ); disp(dP); disp( 无功不平衡量： ); disp(dQ); if (max(abs(d P)limit en d%if en d%for %迭代结束，判断收敛 if (max(abs(d P)limit else dis p(计算不收敛或未达到要求精度); end

13、 %打印功率 fprintf(迭代总次数： disp(节点电压幅值： disp (U); dis p(节点电压相角： dis p( rad2deg(cita); disp(有功计算结果： dis p(Pi); disp(无功计算结果： disp (Qi); %dn, i); )； ); ); ); 子程序一 % file name: Un bala nced.m % author:山东科技大学罗江 % function:计算功率不平衡量 fun ctio n dP ,dQ, Pi,Qi = Un bala nced( n,m, P, Q,U,G,B,cita ) %计算 Pi有功的不平衡量 f

14、or i=1: n for j=1: n Pn (j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j)+B(i,j)*s in (cita(i)-cita(j); end Pi(i)=sum( Pn); end dP=P(1:n-1)-Pi(1:n-1); %dP 有 n-1 个 %计算 Qi无功的不平衡量 for i=1: n for j=1: n Qn (j)=U(i)*U(j)*(G(i,j)*si n(cita(i)-cita(j)-B(i,j)*cos(cita(i)-cita(j); end Qi(i)=sum(Q n); end dQ=Q(1:m)-Q

15、i(1:m); %dQ 有 m 个 en d%fu nc 子程序二 % file nameJacobi.m % author:山东科技大学罗江 % function:计算雅可比矩阵 fun ctio n J = Jacobi( n,m,U,cita,B,G, Pi,Qi ) %雅可比矩阵的计算 %分块H N KL %i!=j 时 for i=1: n-1 for j=1: n-1 H(i,j)=-U(i)*U(j)*(G(i,j)*s in (cita(i)-cita(j)-B(i,j)*cos(cita(i)-cita(j); end end for i=1: n-1 for j=1:m N

16、(i,j)=-U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j)+B(i,j)*s in (cita(i)-cita(j); end end for i=1:m for j=1: n-1 K(i,j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j)+B(i,j)*si n(cita(i)-cita(j); end end for i=1:m for j=1:m L(i,j)=-U(i)*U(j)*(G(i,j)*s in (cita(i)-cita(j)-B(i,j)*cos(cita(i)-cita(j); end end %i=j 时 f

17、or i=1: n-1 H(i,i)=U(i)A2*B(i,i)+Qi(i); end for i=1:m N(i,i)=-U(i)A2*G(i,i)-Pi(i); end for i=1:m K(i,i)=U(i)A2*G(i,i)-Pi(i); end for i=1:m L(i,i)=U(i)A2*B(i,i)-Qi(i); end %合成雅可比矩阵 J=H N;K L; end 子程序三 % file name:Correct.m % author:山东科技大学罗江 % function:修正节点电压 function dU,dcita = Correct( n,m,U,dP ,dQ,J ) %求解节点电压修正量 for i=1:m Ud2(i,i)=U(i); end dP Q=d P dQ; dUcita=(-i nv(J)*d PQ); dcita=dUcita(1: n-1); dcita=dcita 0; dU=(Ud2*dUcita( n:n+m-1); dU=dU zeros(1, n-m); end 子程序四 % file name: PQ_LJ.m % aut