此程序经40余同学使用检验，无误。
这是一个电气狗熬两个礼拜图书馆的成果，根据华中科技大学《电力系统分析》中原理编写，可用牛顿-拉夫逊和PQ分解法计算给定标幺值条件的潮流。
本人水平有限，仅供参考，欢迎一起找Bug。
2018/07/06 说明：由于本人变压器建模与PSASP不同，本人使用模型如下图，参数输入时请按该模型计算。
2018/06/18 主程序更新：增加补偿电容参数主程序% file name:chaoliu_lj.m% auther: 山东科技大学罗江% function:使用牛顿-拉夫逊法、PQ分解法计算潮流% updata:2018/6/18 13:22 增加补偿电容参数%节点类型标号%PQ节点 1%PV节点 2%slack节点 3%能计算给定标幺值网络，有且仅有一个平衡节点的潮流%注意：母线标号顺序要求：PQ节点-PV节点-平衡节点%若某元件不存在，其导纳为0，阻抗为infclear %清除工作空间变量clc %清屏%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%数据输入（标幺值）SB=100; %基准容量，单位MVA%母线基准电压Bus=[115 10.5 115 115];%交流线参数：I侧母线J侧母线阻抗1/2接地导纳Line=[4 1 0.06125+0.09527i 0;4 3 0.08469+0.12703i 0;1 3 0.13989+0.15501i 0];%变压器参数：I侧母线J侧母线阻抗变比%变压器阻抗归算到I侧Trans=[2 3 0.0137+0.2881i 0.9504];%加接地电容器补偿: 母线导纳Cap=[2 0.5i];%发电机参数：母线节点类型P V/U θGen=[4 3 1 0];%负荷参数：母线节点类型P Q%按参考方向，发电机发出功率（正值），负荷消耗功率（负值）Load=[1 1 -0.18 -0.06;2 1 -0.32 -0.12];mode=1; %1-极坐标下牛拉法，2-PQ分解法Tmax=50; %最大迭代次数limit=1.0e-4; %要求精度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%变压器π型等效阻抗参数Zt=zeros(size(Trans,1),3);Zt(:,1)=Trans(:,3)./Trans(:,4);Zt(:,2)=Trans(:,3)./(1-Trans(:,4));Zt(:,3)=Trans(:,3)./(Trans(:,4).^2-Trans(:,4));Trans_pi=[Trans(:,1:2) Zt(:,1) 1./Zt(:,2) 1./Zt(:,3)];n=numel(Bus); %总节点数m=n-1; %PQ节点数for i=1:size(Gen,1)%数组行数if Gen(i,2)==2 %除去PV节点就是PQ节点m=m-1;endendfor i=1:size(Load,1)if Load(i,2)==2m=m-1;endend%PQ节点包含浮游节点，其PQ=0%提取P,Q,U向量P=zeros(1,n); %P,Q为原始数据，Pi,Qi为计算结果Q=zeros(1,n);U=ones(1,n); %电压初始值由此确定cita=zeros(1,n); %此处未知节点皆设为1.0∠0 %注意：此处角度单位为度，提取后再转换成弧度，后面计算使用弧度for i=1:size(Gen,1)if Gen(i,2)==1 %PQ节点P(Gen(i,1))=Gen(i,3);Q(Gen(i,1))=Gen(i,4);endif Gen(i,2)==2 %PV节点P(Gen(i,1))=Gen(i,3);U(Gen(i,1))=Gen(i,4);endif Gen(i,2)==3 %slack节点U(Gen(i,1))=Gen(i,3);cita(Gen(i,1))=Gen(i,4);endendfor i=1:size(Load,1)if Load(i,2)==1 %PQ节点P(Load(i,1))=Load(i,3);Q(Load(i,1))=Load(i,4);endif Load(i,2)==2 %PV节点P(Load(i,1))=Load(i,3);U(Load(i,1))=Load(i,4);endif Load(i,2)==3 %slack节点U(Load(i,1))=Load(i,3);cita(Load(i,1))=Load(i,4);endenddisp('初始条件：')disp('各节点有功：')disp(P);disp('各节点无功：')disp(Q);disp('各节点电压幅值：')disp(U);cita=(deg2rad(cita)); %角度转换成弧度disp('各节点电压相角(度)：')disp(rad2deg(cita)); %显示依然使用角度%节点导纳矩阵的计算Y=zeros(n); %新建节点导纳矩阵y=zeros(n); %网络中的真实导纳%计算y(i,j)for i=1:size(Line,1) %与交流线联结的真实导纳ii=Line(i,1); jj=Line(i,2);y(ii,jj)=1/Line(i,3);y(jj,ii)=y(ii,jj);endfor i=1:size(Trans_pi,1) %与变压器联结的真实导纳ii=Trans_pi(i,1); jj=Trans_pi(i,2);y(ii,jj)=1/Trans_pi(i,3);y(jj,ii)=y(ii,jj);end%计算y(i,i)for i=1:size(Line,1) %与交流线联结的对地导纳ii=Line(i,1); jj=Line(i,2);y(ii,ii)=y(ii,ii)+Line(i,4);y(jj,jj)=y(jj,jj)+Line(i,4);endfor i=1:size(Trans_pi,1) %与变压器联结的对地导纳ii=Trans_pi(i,1); jj=Trans_pi(i,2);y(ii,ii)=y(ii,ii)+Trans_pi(i,4);y(jj,jj)=y(jj,jj)+Trans_pi(i,5);end%算上补偿电容for i=1:size(Cap,1)ii=Cap(i,1);y(ii,ii)=y(ii,ii)+Cap(i,2);end%由y计算Yysum=sum(y,1); %每一行求和for i=1:nfor j=1:nif i==jY(i,j)=ysum(i);elseY(i,j)=-y(i,j);endendenddisp('节点导纳矩阵：');disp(Y);G=real(Y); %电导矩阵B=imag(Y); %电纳矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%以上为基础数据整理%接下来是牛拉法的大循环%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%计算功率不平衡量[dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita );disp('有功不平衡量：');disp(dP);disp('无功不平衡量：');disp(dQ);for i=1:Tmaxfprintf('第%d次迭代\n',i);%雅可比矩阵的计算if(mode==1)J=Jacobi( n,m,U,cita,B,G,Pi,Qi );disp('雅可比矩阵');disp(J);end%求解节点电压修正量if(mode==1)[dU,dcita]=Correct( n,m,U,dP,dQ,J );else[dU,dcita]=PQ_LJ( n,m,dP,dQ,U,B );enddisp('电压、相角修正量：');disp(dU);disp(rad2deg(dcita));%修正节点电压U=U+dU;cita=cita+dcita;disp('节点电压幅值：');disp(U);disp('节点电压相角：');disp(rad2deg(cita));%计算功率不平衡量[dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita );disp('有功不平衡量：');disp(dP);disp('无功不平衡量：');disp(dQ);if (max(abs(dP))<limit && max(abs(dQ))<limit )break;end%ifend%for%迭代结束，判断收敛if (max(abs(dP))<limit && max(abs(dQ))<limit )disp('计算收敛');elsedisp('计算不收敛或未达到要求精度');end%打印功率fprintf('迭代总次数：%d\n', i);disp('节点电压幅值：');disp(U);disp('节点电压相角：');disp(rad2deg(cita));disp('有功计算结果：');disp(Pi);disp('无功计算结果：');disp(Qi);子程序一% filename:Unbalanced.m% author: 山东科技大学罗江% function: 计算功率不平衡量function [ dP,dQ,Pi,Qi ] = Unbalanced( n,m,P,Q,U,G,B,cita )%计算ΔPi有功的不平衡量for i=1:nfor j=1:nPn(j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));endPi(i)=sum(Pn);enddP=P(1:n-1)-Pi(1:n-1); %dP有n-1个%计算ΔQi无功的不平衡量for i=1:nfor j=1:nQn(j)=U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));endQi(i)=sum(Qn);enddQ=Q(1:m)-Qi(1:m); %dQ有m个end%func子程序二% filename:Jacobi.m% author:山东科技大学罗江% function: 计算雅可比矩阵function [ J ] = Jacobi( n,m,U,cita,B,G,Pi,Qi )%雅可比矩阵的计算%分块H N K L%i!=j时for i=1:n-1for j=1:n-1H(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));endendfor i=1:n-1for j=1:mN(i,j)=-U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));endendfor i=1:mfor j=1:n-1K(i,j)=U(i)*U(j)*(G(i,j)*cos(cita(i)-cita(j))+B(i,j)*sin(cita(i)-cita(j)));endendfor i=1:mfor j=1:mL(i,j)=-U(i)*U(j)*(G(i,j)*sin(cita(i)-cita(j))-B(i,j)*cos(cita(i)-cita(j)));endend%i==j时for i=1:n-1H(i,i)=U(i).^2*B(i,i)+Qi(i);endfor i=1:mN(i,i)=-U(i).^2*G(i,i)-Pi(i);endfor i=1:mK(i,i)=U(i).^2*G(i,i)-Pi(i);endfor i=1:mL(i,i)=U(i).^2*B(i,i)-Qi(i);end%合成雅可比矩阵J=[H N;K L];end子程序三% filename:Correct.m% author:山东科技大学罗江% function:修正节点电压function [ dU,dcita ] = Correct( n,m,U,dP,dQ,J )%求解节点电压修正量for i=1:mUd2(i,i)=U(i);enddPQ=[dP dQ]';dUcita=(-inv(J)*dPQ)';dcita=dUcita(1:n-1);dcita=[dcita 0];dU=(Ud2*dUcita(n:n+m-1)')';dU=[dU zeros(1,n-m)];end子程序四% filename:PQ_LJ.m% author:山东科技大学罗江% function:使用PQ分解法计算电压修正量function [ dU,dcita ] = PQ_LJ( n,m,dP,dQ,U,B )dP_U=dP./U(1:n-1);dQ_U=dQ./U(1:m);dUdcita=(-inv(B(1:n-1,1:n-1))*dP_U')';dcita=dUdcita./U(1:n-1);dU=(-inv(B(1:m,1:m))*dQ_U')';dU=[dU zeros(1,n-m)];dcita=[dcita 0];%补零end （使用时此括号删去。
