《matlab潮流计算》由会员分享，可在线阅读，更多相关《matlab潮流计算(14页珍藏版)》请在人人文库网上搜索。

1、附录 1使用牛顿拉夫逊法进行潮流计算的Matlab 程序代码% 牛拉法计算潮流程序%-% B1 矩阵： 1、支路首端号； 2、末端号； 3、支路阻抗； 4、支路对地电纳%5、支路的变比； 6、支路首端处于 K 侧为 1，1 侧为 0% B2 矩阵： 1、该节点发电机功率； 2、该节点负荷功率； 3、节点电压初始值% 4、PV 节点电压 V 的给定值； 5、节点所接的无功补偿设备的容量% 6、节点分类标号： 1 为平衡节点(应为 1 号节点)；2 为 PQ 节点； 3 为 PV 节点；%-clear all;format long;n=input( 请输入节点数 :nodes=);nl=inpu。

2、t( 请输入支路数:lines=);isb=input( 请输入平衡母线节点号:balance=);pr=input( 请输入误差精度:precision=);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);%-for i=1:nl%支路数if B1(i,6)=0%左节点处于1 侧p=B1(i,1);q=B1(i,2);else% 左节点处于K 侧p=B1(i,2);q。

3、=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);%对角元 K 侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4);%对角元 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);专业文档供参考，如有帮。

4、助请下载。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; %迭代次数ICT1 、 a；不满足收敛要求的节点数IT2while IT2=0% N0=2*n雅可比矩阵的阶数；N=N0+1 扩展列IT2=0;a=a+1;for i=1:nif i=isb%非平衡节点C(i)=。

5、0;D(i)=0;for j1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);% (Gij*ej-Bij*f)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)%求 i 节点有功和无功功率P,Q的计算值V2=e(i。

6、)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/dfX3=dQ/deX4=-X1;% X1=dP。

7、/deX4=dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;% X3=dQ/deJ(p,N)=DQ节点无功功率差J(m,q)=X1;J(m,N)=DP;q=q+1;% X1=dP/deJ(m,N)=DP节点有功功率差J(p,q)=X4;J(m,q)=X2;% X4=dQ/dfX2=dp/dfelseif 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/df专业文档供参考，如有帮助请下载。X3=D。

8、(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。

9、=-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;% PV 节点电压误差m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;% PV 节点有功误差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/。

10、dfX5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV 节点电压误差m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;% PV 节点有功误差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、 Q 或 Ufor k2=k1:N1% 从 k。

11、+1 列的 Jacobi 元素到扩展列的P、 Q 或UJ(k,k2)=J(k,k2)./J(k,k);%用 K 行 K 列对角元素去除K 行 K 列后的非对角元素进行规格化专业文档供参考，如有帮助请下载。endJ(k,k)=1;% 对角元规格化K 行 K 列对角元素赋1%回代运算 -if k=3% 不是第三行k 3k4=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);%消去运算(当前行k 列元素消为0)end%用当前行K2 列元素减去当前行。

12、k 列元素乘以第k 行 K2 列元素J(k3,k)=0; % 当前行第k 列元素已消为0endif k=N0%若已到最后一行break;end%前代运算 -for k3=k1:N0% 从 k+1 行到 2*n 最后一行for k2=k1:N1%从 k+1 列到扩展列消去k+1 行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算end%用当前行K2 列元素减去当前行k 列元素乘以第k 行 K2 列元素J(k3,k)=0; % 当前行第k 列元素已消为0endelse%是第三行k=3%第三行 k=3 的前代运算 -for k3=k1:N0%从第四行到2。

13、n 行(最后一行)for k2=k1:N1%从第四列到2n+1 列(即扩展列)J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算(当前行3 列元素消为0)end%用当前行K2 列元素减去当前行3 列元素乘以第三行K2 列元素J(k3,k)=0; % 当前行第3 列元素已消为0endendend%上面是用线性变换方式高斯消去法将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:N。

14、0DET=abs(J(k,N);专业文档供参考，如有帮助请下载。if DET=pr%电压偏差量是否满足要求IT2=IT2+1; % 不满足要求的节点数加1endendICT2(a)=IT2;%不满足要求的节点数ICT1=ICT1+1;% 迭代次数end%用高斯消去法解w=-J*Vdisp( 迭代次数： );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)*1i;。

15、%将各节点电压用复数表示end%计算各输出量 -disp( 各节点的实际电压标幺值E 为： );disp(E);%显示各节点的实际电压标幺值E 用复数表示disp(-);disp( 各节点的电压大小V 为： );disp(V);%显示各节点的电压大小V 的模值disp(-);disp( 各节点的电压相角deg 为： );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 注入电流的共轭值。

16、enddisp( 各节点的功率S 为： );disp(S);% 显示各节点的注入功率disp(-);disp( 各条支路的首端功率Si 为： );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);else专业文档供参考，如有帮助请下载。Si(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)./2)+(conj(E(p)./B。

17、1(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);ZF=S(,num2str(p),num2str(q),)=,num2str(SSi(p,q); disp(ZF);disp(-);enddisp( 各条支路的末端功率Sj 为： );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)*con。

18、j(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);ZF=S(,num2str(q),num2str(p),)=,num2str(SSj(q,p); disp(ZF);disp(-);enddisp( 各条支路的功率损耗DS 为： );for i=1:nlp=B1(i,1);q=B。

19、1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i);DDS(i)=DS(i);ZF=DS(,num2str(p),num2str(q),)=,num2str(DDS(i);disp(ZF);disp(-);end专业文档供参考，如有帮助请下载。附录 2使用 PQ 分解法进行潮流计算的Matlab 程序代码%PQ 分解法潮流计算程序%本文中的实例数据如下：节点数为 9；支路数为 9；平衡母线节点号为 1；误差精度为 0.00001；PQ 节点数为 5；%主程序clear all;format long;n=input(请输入节点数： n=);nl=input( 请输入。

20、支路数： nl=);isb=input(请输入平衡母线节点号：isb=);pr=input(请输入误差精度： pr=);B1=input( 请输入由支路参数形成的矩阵：B1=); % 输入 B1B2=input( 请输入由支路参数形成的矩阵：B2=); % 输入 B2na=input(请输入 PQ 节点数 na=);Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n);for i=1:nlif B1(i,6)=0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1)。

21、;endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5);YI(p,q)=YI(p,q)-1./B1(i,3);Y(q,p)=Y(p,q);YI(q,p)=YI(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)2)+B1(i,4);YI(q,q)=YI(q,q)+1./B1(i,3);Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4);YI(p,p)=YI(p,p)+1./B1(i,3);end %求导纳矩阵disp(节点导纳矩阵为： );disp(Y);G=real(Y);B=imag(YI);BI=imag(Y);for i=1:。

22、nS(i)=B2(i,1)-B2(i,2);BI(i,i)=BI(i,i)+B2(i,5);endP=real(S);Q=imag(S);专业文档供参考，如有帮助请下载。for i=1:ne(i)=real(B2(i,3);f(i)=imag(B2(i,3);V(i)=B2(i,4);endfor i=1:nif B2(i,6)=2V(i)=sqrt(e(i)2+f(i)2);O(i)=atan(f(i)./e(i);endendfor i=2:nif i=nB(i,i)=1./B(i,i);else IC1=i+1;for j1=IC1:nB(i,j1)=B(i,j1)./B(i,i);en。

23、dB(i,i)=1./B(i,i);for k=i+1:nfor j1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1);endendendendp=0;q=0;for i=1:nif B2(i,6)=2p=p+1;k=0;for j1=1:nif B2(j1,6)=2k=k+1;A(p,k)=BI(i,j1);endendendendfor i=1:naif i=naA(i,i)=1./A(i,i);else k=i+1;for j1=k:naA(i,j1)=A(i,j1)./A(i,i);专业文档供参考，如有帮助请下载。endA(i,i)=1./A(i,i);for 。

24、k=i+1:nafor j1=i+1:naA(k,j1)=A(k,j1)-A(k,i)*A(i,j1);endendendendICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;while ICT2=0|ICT3=0ICT2=0;ICT3=0;for i=1:nif i=isbC(i)=0;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k)+BI(i,k)*sin(O(i)-O(k);endDP1(i)=P(i)-V(i)*C(i);DP(i)=DP1(i)./V(i);DET=abs(DP1(i);if DET=prIC。

25、T2=ICT2+1;endendendNp(K)=ICT2;if ICT2=0for i=2:nDP(i)=B(i,i)*DP(i);if i=nIC1=i+1;for k=IC1:nDP(k)=DP(k)-B(k,i)*DP(i);endelsefor LZ=3:iL=i+3-LZ;IC4=L-1;for MZ=2:IC4I=IC4+2-MZ;DP(I)=DP(I)-B(I,L)*DP(L);endend专业文档供参考，如有帮助请下载。endendfor i=2:nO(i)=O(i)-DP(i);endkq=1;L=0;for i=1:nif B2(i,6)=2C(i)=0;L=L+1;fo。

26、r k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k)-BI(i,k)*cos(O(i)-O(k); endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);DET=abs(DQ1(i);if DET =prICT3=ICT3+1;endendendelse kp=0;if kq=0;L=0;for i=1:nif B2(i,6)=2C(i)=0;L=L+1;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k)-BI(i,k)*cos(O(i)-O(k); endDQ1(i)=Q(i)-V(i。

27、)*C(i);DQ(L)=DQ1(i)./V(i);DET=abs(DQ1(i);endendendendNq(K)=ICT3;if ICT3=0L=0;for i=1:naDQ(i)=A(i,i)*DQ(i);if i=na专业文档供参考，如有帮助请下载。for LZ=2:iL=i+2-LZ;IC4=L-1;for MZ=1:IC4I=IC4+1-MZ;DQ(I)=DQ(I)-A(I,L)*DQ(L);endendelseIC1=i+1;for k=IC1:naDQ(k)=DQ(k)-A(k,i)*DQ(i);endendendL=0;for i=1:nif B2(i,6)=2L=L+1;V。

28、(i)=V(i)-DQ(L);endendkp=1;K=K+1;elsekq=0;if kp=0K=K+1;endendfor i=1:nDy(K-1,i)=V(i);endenddisp(迭代次数 );disp(K);disp(每次没有达到精度要求的有功功率个数为);disp(Np);disp(每次没有达到精度要求的无功功率个数为);disp(Nq);for k=1:nE(k)=V(k)*cos(O(k)+V(k)*sin(O(k)*j;O(k)=O(k)*180./pi;end专业文档供参考，如有帮助请下载。disp(各节点的电压标幺值E 为： );disp(E);disp(各节点的电压 。

29、V 大小为： );disp(V);disp(各节点的电压相角O 为： );disp(O);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(各条支路的首端功率Sj 为： );for i=1:nlif B1(i,6)=0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSi(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5。

30、)-conj(E(q)*conj(1./( B1(i,3)*B1(i,5);disp(Si(p,q);enddisp(各条支路的末端功率Sj 为： );for i=1:nlif B1(i,6)=0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSj(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);disp(Sj(q,p);enddisp(各条支路的功率损耗DS 为： );for i=1:nlif B1(i,6)。

31、=0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endDS(i)=Si(p,q)+Sj(q,p);专业文档供参考，如有帮助请下载。disp(DS(i);endfor i=1:KCs(i)=i;for j=1:nDy(K,j)=Dy(K-1,j);endend附录 3进行三相短路容量计算的Matlab 程序代码%短路容量计算程序%-% B1 矩阵： 1、支路首端号； 2、末端号； 3、支路阻抗； 4、支路对地电纳% 5、支路的变比； 6、支路首端处于 K 侧为 1，1 侧为 0% B2 矩阵： 1、该节点发电机功率； 2、该节点负荷功率； 3、节点电。

32、压初始值% 4、PV 节点电压 V 的给定值； 5、节点所接的无功补偿设备的容量% 6、节点分类标号： 1 为平衡节点(应为 1 号节点)； 2 为 PQ 节点； %Yd 为修改后的节点导纳矩阵%-clear all;format long;g1=input(300MW 发电机数： g1=);g2=input(250MW 发电机数： g2=);n=input(请输入节点数 :n=);nl=input( 请输入支路数 :nl=);B1=input( 请输入由各支路参数形成的矩阵:B1=);B2=input( 请输入各节点参数形成的矩阵:B2=);Y=zeros(n);% Y 为修改前节点导纳矩阵。

33、for i=1:nl%支路数if B1(i,6)=0%左节点处于 1 侧p=B1(i,1);q=B1(i,2);else%左节点处于 K 侧p=B1(i,2);q=B1(i,1);end专业文档供参考，如有帮助请下载。Y(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);%对角元 K 侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4);%对角元 1 侧end%-%Y2-Y5 为各 PQ 节点负荷的导纳Y1=0;Y2=conj(B2。

34、(2,2);Y3= conj(B2(3,2);Y4= conj(B2(4,2);Y5= conj(B2(5,2); Xd300=0.51j/0.9504132;XT300=0.033212j/0.9504132;Xd250=0.714j/0.9504132;XT250=0.038747j/0.9504132;Y6=g1/(XT300+Xd300)+g2/(XT250+Xd250);%处理相应的负荷及机组部分的导纳%-C=Y1,Y2,Y3,Y4,Y5,Y6;Yd=Y;for i=1:nYd(i,i)=Yd(i,i)+C(i);%修改各节点自导纳enddisp(Yd);Z = inv(Yd);%求节点阻抗矩阵for j=1:nI(j) = 1/Z(j,j);%电压故障前电压标幺值为1S(j)=abs(I(j);Sn(j)=S(j)*100;%短路电流有名值end%计算完毕 -disp(各节点短路时的短路电流幅值标幺值)disp(abs(I)disp(短路容量有名值Sn=);disp(Sn);专业文档供参考，如有帮助请下载。