1、此程序经40余同学使用检验，无误。这是一个电气狗熬两个礼拜图书馆的成果，根据华中科技大学电力系统分析中原理编写，可用牛 顿-拉夫逊和PQ分解法计算给定标幺值条件的潮流。本人水平有限，仅供参考，欢迎一起找Bug2018/07/06说明：由于本人变压器建模与 PSASP不同，本人使用模型如下图，参数输入时请按该 模型计算。Zt 口| | i -1=1_GD j2018/06/18主程序更新：增加补偿电容参数主程序% file name:chaoliu_lj.m% auther:山东科技大学罗江% function:使用牛顿-拉夫逊法、PQ分解法计算潮流% updata:2018/6/18 13:2
2、2增加补偿电容参数%节点类型标号%PQ节点1%PV节点2%slack 节点 3%能计算给定标幺值网络，有且仅有一个平衡节点的潮流%注意：母线标号顺序要求：PQ节点-PV节点-平衡节点%若某元件不存在，其导纳为0,阻抗为inf clear %清除工作空间变量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;%变压器参数：
3、I侧母线J侧母线 阻抗 变比变压器阻抗归算到I侧Trans=2 3 0.0137+0.2881i 0.9504;%加接地电容器补偿：母线导纳Cap=2 0.5i;%发电机参数：母线节点类型P V/U 0Gen=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(:
4、1)=Trans(:,3)./Trans(:,4);Zt(:,2)=Trans(:,3)./(1-Trans(:,4);Zt(:,3)=Trans(:,3)./(Trans(:,4).A2-Trans(:,4);Trans_pi=Trans(:,1:2) Zt(:,1) 1./Zt(:,2) 1./Zt(:,3);n=numel(Bus); %总节点数m=n-1; %PQ节点数fo门=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;enden
5、d%PQ节点包含浮游节点，其 PQ=0%提取PQ,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)=Gen(i,4);endif Gen(i,2)=2 %PV 节点 P(Gen(i,1)=Gen(i,3);U(Gen(i,1)=Gen(
6、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)
7、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
8、/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)
9、 %与变压器联结的对地导纳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);
10、G=real(Y); %电导矩阵B=imag(Y); %fe 纳矩阵% %以上为基础数据整理%接下来是牛拉法的大循环% %计算功率不平衡量dP,dQ,Pi,Qi=Unbalanced( n,m,P,Q,U,G,B,cita );disp(有功不平衡量：);disp(dP);disp(无功不平衡量：,);disp(dQ);for i=1:Tmaxfprintf(第 次迭代 n,i);%雅可比矩阵的计算if(mode=1)J=Jacobi( n,m,U,cita,B,G,Pi,Qi );disp(雅可比矩阵,);disp(J);end%求解节点电压修正量if(mode=1)dU,dcita=Cor
11、rect( n,m,U,dP ,dQ,J );elsedU,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 (
12、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(迭代总次数：%dn, i);disp(节点电压幅值：,);disp(U);disp(节点电压相角:);disp(rad2deg(cita);disp(有功计算结果：);disp(Pi);disp(无功计算结果：,);disp(Qi);子程序一% filename:Unbalan
13、ced.m% author:山东科技大学罗江% function:计算功率不平衡量function dP,dQ,Pi,Qi = Unbalanced( n,m,P,Q,U,G,B,cita )%计算A 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 个%计算A Qi无功的不平衡量for i=1:nfor j=1:nQn(j)=U(i)*U(j
14、)*(G(i,j)*sin(cita(i)-cita(j)-B(i,j)*cos(cita(i)-cita(j); endQi(i)=sum(Qn);end dQ=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(
15、i)-cita(j)-B(i,j)*cos(cita(i)-cita(j); end end for i=1:n-1 for 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); end end for i=1:m for 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); end end for i=1:mfor j=1:mL(i,j)=-U(i)*U(j)*(G(i,j)*sin(
16、cita(i)-cita(j)-B(i,j)*cos(cita(i)-cita(j); end end %i=j 时 for i=1:n-1H(i,i)=U(i).A2*B(i,i)+Qi(i);end for i=1:mN(i,i)=-U(i).A2*G(i,i)-Pi(i);end for i=1:mK(i,i)=U(i).A2*G(i,i)-Pi(i);end for i=1:mL(i,i)=U(i)A2*B(i,i)-Qi(i);end%合成雅可比矩阵J=H N;K L;end子程序三% filename:Correct.m% author:山东科技大学 罗江% function:修正
17、节点电压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 （使用时此括号删去。版权所有人：罗江）
