潮流计算matlab牛顿拉夫逊法法.docx
- 文档编号:11519524
- 上传时间:2023-06-01
- 格式:DOCX
- 页数:48
- 大小:125.23KB
潮流计算matlab牛顿拉夫逊法法.docx
《潮流计算matlab牛顿拉夫逊法法.docx》由会员分享,可在线阅读,更多相关《潮流计算matlab牛顿拉夫逊法法.docx(48页珍藏版)》请在冰点文库上搜索。
潮流计算matlab牛顿拉夫逊法法
用matlab潮流计算(牛顿拉夫逊法)
%主程序:
[d]=uigetfile('ieee14.m','SelectDataFile');ifpathname==0
error('youmustselectavaliddatafile')else
l(dfile);
%stripoff.m
eval(d));%打开数据文件end
endend;
bus=[PQ;PV;SW];
newbus=[1:
nb]';
f=bus(:
1);
nodenum=[newbusbus(:
1)];
bus(:
1)=newbus;
fori=1:
nl
forj=1:
2
fork=1:
nb
ifline(i,j)==nodenum(k,2)line(i,j)=nodenum(k,1);break
end
end
end
end
Y=y(bus,line);%形成节点导纳矩阵K=0;%迭代次数初值
Kmax=10;%最大迭代次数
eps1=1.0e-10;
eps2=1.0e-10;
m=nPQ;n=nb;
Um=eye(m,m);
myf=fopen('output1.dat','w');
forK=1:
Kmax
fori=1:
m
forj=1:
m
ifi==j
Um(i,j)=bus(i,2);
end
end
endb=dPQ(Y,bus);C=jac(bus,Y);dX=C\b';dx=dX';
[nx,mx]=size(dx);
fori=1:
n-1%计算相角
bus(i,3)=bus(i,3)-dX(i,1);
end
B=dx(nx,n:
mx)*Um;%计算电压差
bus(1:
m,2)=bus(1:
m,2)-B';%计算电压值
dx(nx,n:
mx)=B;
fprintf(myf,'--第%d次迭代时雅可比矩阵--',K)fprintf(myf,'\n');
fori=1:
(n+m-1)
forj=1:
(n+m-1)
fprintf(myf,'%8.6f',C(i,j));fprintf(myf,'');
end
fprintf(myf,'\n');
end
fprintf(myf,'--第%d次迭代时dPQ的误差--',K)fprintf(myf,'\n');
fori=1:
(n+m-1)
fprintf(myf,'%8.6e',b(1,i));fprintf(myf,'\n');end
fprintf(myf,'\n');
fprintf(myf,'--第%d次迭代时dx(误差)--',K)fprintf(myf,'\n');
fori=1:
(n+m-1)
fprintf(myf,'%8.6e',dX(i,1));fprintf(myf,'\n');
end
fprintf(myf,'\n');
fprintf(myf,'第%d次迭代后节点电压(仅PQ节点儿K)fprintf(myf,'\n');
fori=1:
m
fprintf(myf,'%8.6f',bus(i,2));
fprintf(myf,'\n');
end
fprintf(myf,'\n');
fprintf(myf,'第%d次迭代后相角(角度儿K)fprintf(myf,'\n');
fori=1:
n
fprintf(myf,'%8.6f',bus(i,3)*180/pi);
fprintf(myf,'\n');
endfprintf(myf,'\n');
if(max(abs(dx)) end end %计算功率 P1=0;T=0;%计算平衡节点的有功 forj=1: nT=bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+imag(Y(n,j))*sin(bus(n,3)-bus(j,3)));P1=P1+T; end bus(n,4)=P1; fork=m+1: n%计算PV节点、平衡节点的无功 Q1=0;T=0; forj=1: nT=bus(k,2)*bus(j,2)*(real(Y(k,j))*sin(bus(k,3)-bus(j,3))-imag(Y(k,j))*cos(bus(k,3)-bus(j,3)));Q1=Q1+T; end bus(k,5)=Q1; end bus(: 1)=f;%换回各节点、支路的初始编号 r=zeros(1,mb); fort=1: n forl=t+1: n ifbus(t,1)>bus(l,1)r=bus(t,: );bus(t,: )=bus(l,: );bus(l,: )=r; end end end fori=1: nl forj=1: 2 fork=1: nb ifline(i,j)==nodenum(k,1) line(i,j)=nodenum(k,2); break end end end end fclose(myf); Pf=loss(bus,line);%计算支路潮流及损耗%将节点导纳矩阵、节点潮流计算结果写入文件output2 myf=fopen('output2.dat','w'); fprintf(myf,'--节点导纳矩阵--\n'); fork=1: n forj=1: n fprintf(myf,'%8.6f',real(Y(k,j))); fprintf(myf,'+i*'); fprintf(myf,'%8.6f',imag(Y(k,j))); fprintf(myf,''); end fprintf(myf,'\n');end fprintf(myf,'牛顿-拉夫逊法潮流计算结果\n'); fprintf(myf,'节点计算结果\n'); dS(I,J)--\n'); fork=1: nl forj=1: 5 ifj<=2 fprintf(myf,'%8.1f',Pf(k,j)); fprintf(myf,'');else fprintf(myf,'%8.6f',real(Pf(k,j)));fprintf(myf,'+i*'); fprintf(myf,'%8.6f',imag(Pf(k,j))); fprintf(myf,''); end -line R=X=0 endfprintf(myf,'\n'); endfclose(myf); %根据支路参数建立节点导纳矩阵程序: functionY=y(bus,line)%目的: 根据支路参数建立节点导纳矩阵%输入: 节点参数矩阵--bus;支路参数矩阵%输出: 节点导纳矩阵--Y[nb,mb]=size(bus); [nl,ml]=size(line); Y=zeros(nb,nb); fork=1: nl I=line(k,1); J=line(k,2); Zt=line(k,3)+j*line(k,4); ifZt==0 disp('对地支路'); Yt=inf; else Yt=1/Zt;endYm=line(k,5)+j*line(k,6);K=line(k,7); if(K==0)&(J~=0)%普通线路 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt; Y(J,I)=Y(I,J); endif(K==0)&(J==0)%对地支路K=J=0, Y(I,I)=Y(I,I)+Ym; end ifK>0%K>0时变压器支路 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt/O2; Y(I,J)=Y(I,J)-Yt/K; Y(J,I)=Y(I,J); endifK<0%K<0时变压器支路 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt*KA2; Y(I,J)=Y(I,J)+Yt*K; Y(J,I)=Y(I,J); end end %形成雅克矩阵程序: functionJ=jac(bus,Y) %形成雅可比矩阵 globaln; globalm; fori=1: (n-1)%计算PQ、PV节点的有功 P1=0;T=0; forj=1: nT=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));P1=P1+T; end bus(i,4)=P1; end fori=1: n-1%计算PV、PQ节点的无功 Q1=0;T=0; forj=1: nT=bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));Q1=Q1+T; end bus(i,5)=Q1; end fori=1: n-1%计算H forj=1: n-1 ifi~=jH(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));N(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));K(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));L(i,j)=H(i,j); end end end fori=1: n-1%计算H forj=1: n-1 ifj==i H(i,i)=bus(i,5)+imag(Y(i,i))*bus(i,2)A2; N(i,i)=-bus(i,4)-real(Y(i,i))*bus(i,2)A2; K(i,i)=N(i,i)+2*real(Y(i,i))*bus(i,2)A2; L(i,i)=-H(i,i)+2*imag(Y(i,i))*bus(i,2)A2; end end end N=N(1: n-1,1: m); K=K(1: m,1: n-1); L=L(1: m,1: m); J=[H,N;K,L]; %计算dPQ的程序: functiondPQ=dPQ(Y,bus) %delP--有功偏差量 %delQ--无功偏差量 %Y--节点导纳矩阵 %bus--节点参数(P,Q,U及相角)矩阵 globaln; globalm; delP=zeros(1,n-1); delQ=zeros(1,m); fori=1: (n-1)%形成delP矩阵(PQ、PV节点) P1=0;T=0; forj=1: n T=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3))); P1=P1+T; end delP(1,i)=bus(i,4)-P1; end fori=1: m%形成delQ矩阵(PQ节点) Q1=0;T=0; forj=1: n T=bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3))); Q1=Q1+T; end delQ(1,i)=bus(i,5)-Q1; end dPQ=[delP,delQ]; %计算线路损耗、线路潮流程序: functionPf=loss(bus,line) %计算线路损耗、线路潮流 [nl,ml]=size(line); Pf=zeros(nl,5); fork=1: nl I=line(k,1); J=line(k,2); Zt=line(k,3)+i*line(k,4); ifZt==0 Yt=inf; else Yt=1/Zt; end Ym=line(k,5)+i*line(k,6); K=line(k,7); S(l,J)=bus(l,2)A2*(conj(Yt)+conj(Ym))-bus(l,2)*(cos(bus(l,3))+i*sin(bus(l,3)))*bus(J,2)*(cos(bus(J,3))-i*sin(bus(J,3)))*conj(Yt); S(J,l)=bus(J,2)A2*(conj(Yt)+conj(Ym))-bus(J,2)*(cos(bus(J,3))+i*sin(bus(J,3)))*bus(l,2)*(cos(bus(l,3))-i*sin(bus(l,3)))*conj(Yt); delS(l,J)=S(l,J)+S(J,l); end if(K==0)&(J==0)%对地支路潮流 J=5;S(l,5)=bus(l,2)A2*conj(Ym); end ifK>0%变压器支路k>0时的潮流 S(l,J)=bus(l,2)A2*(conj(Ym+Yt*(1-1/K))+conj(Yt/K))-bus(l,2)*(cos(bus(l,3))+i*sin(bus(l,3)))*bus(J,2)*(cos(bus(J,3))-i*sin(bus(J,3)))*conj(Yt/K); S(J,l)=bus(J,2)A2*(conj(Yt))/KA2-bus(J,2)*(cos(bus(J,3))+i*sin(bus(J,3)))*bus(l,2)*(cos(bus(l,3))-i*sin(bus(l,3)))*conj(Yt/K); delS(l,J)=S(l,J)+S(J,l); end ifK<0%变压器支路k<0时的潮流 S(l,J)=bus(l,2)A2*(conj(Ym+Yt))+bus(l,2)*(cos(bus(l,3))+i*sin(bus(l,3)))*bus(J,2)*(cos(bus(J,3))-i*sin(bus(J,3)))*conj(Yt*K); S(J,l)=bus(J,2)A2*(conj(Yt))*KA2+bus(J,2)*(cos(bus(J,3))+i*sin(bus(J,3)))*bus(l,2)*(cos(bus(l,3))-i*sin(bus(l,3)))*conj(Yt*K); delS(l,J)=S(l,J)+S(J,l); end ifJ==5&Zt==0 Sp=[line(k,1)line(k,2)S(l,5)0S(l,5)]; else Sp=[line(k,1)line(k,2)S(l,J)S(J,l)delS(l,J)]; end Pf(k,: )=Sp; end %输入的参数数据: %datafortestcase %各节点参数: 节点编号,注入有功,注入无功,(Sn=100MVA)电压幅值,电压相位,类型 %类型: 1=PQ节点,2=PV节点,3=平衡节点 %(bus#)(volt)(ang)(p)(q)(bustype) bus=[ 1,1.0,0.0,-0.478,0.039,1; 2,1.0,0.0,-0.076,-0.016,1; 3,1.0,0.0,0.0,0.0,1; 4,1.0,0.0,-0.295,-0.166,1; 5,1.0,0.0,-0.09,-0.058,1; 6,1.0,0.0,-0.035,-0.018,1; 7,1.0,0.0,-0.061,-0.016,1; 8,1.0,0.0,-0.135,-0.058,1; 9,1.0,0.0,-0.149,-0.05,1; 10,1.045,0.0,0.183,0.0,2; 11,1.010,0.0,-0.942,0.0,2; 12,1.70,0.0,-0.112,0.047,2; 13,1.90,0.0,0.0,0.174,2; 14,1.060,0.0,0.0,0.0,3]; %各支路参数: 起点编号,终点编号,电阻,电抗,电导,电纳line=[ 1,2,0.01335,0.04211,0.0,0.0,0; 1,3,0.0,0.20912,0.0,0.0,0; 1,4,0.0,0.55618,0.0,0.0,0; 1,10,0.05811,0.17632,0.0,0.0340,0; 1,11,0.06701,0.17103,0.0,0.0128,0; 2,10,0.05695,0.17388,0.0,0.0346,0;2,12,0.0,0.25202,0.0,0.0,0; 2,14,0.05403,0.22304,0.0,0.0492,0; 3,4,0.0,0.11001,0.0,0.0,0; 3,13,0.0,0.17615,0.0,0.0,0; 4,5,0.03181,0.08450,0.0,0.0,0; 4,9,0.12711,0.27038,0.0,0.0,0; 5,6,0.08205,0.19207,0.0,0.0,0;6,12,0.09498,0.19890,0.0,0.0,0; 7,8,0.22092,0.19988,0.0,0.0,0; 7,12,0.12291,0.25581,0.0,0.0,0; 8,9,0.17093,0.34802,0.0,0.0,0; 8,12,0.06615,0.13027,0.0,0.0,0; 10,11,0.04699,0.19797,0.0,0.0438,0; 10,14,0.01938,0.05917,0.0,0.0528,0]; 输出结果数据1: --第1次迭代时雅可比矩阵-- -38.624033 21.578554 4.781943 1.797979 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 5.346051 5.119505 -0.000000 -0.000000 -10.417258 6.840981 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000-0.000000 21.578554 -38.240787 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 5.427654 -0.000000 6.745496 -0.000000 6.840981 -9.429913 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000-0.000000 4.781943 -0.000000 -24.658288 9.090083 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 10.786262 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000-0.000000 1.797979 -0.000000 9.090083-24.282506 10.365394 -0.000000 -0.000000 -0.000000 3.029050 -0.000000 -0.000000 -0.000000-0.000000 -0.000000 -0.000000 -0.000000 -5.326055 3.902050 -0.000000 -0.000000 -0.0000001.424005 -0.000000 -0.000000 -0.00000010.365394 -14.768338 4.402944 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000-0.000000 -0.000000 -0.000000 -0.000000 3.902050 -5.782934 1.880885 -0.000000 -0.000000-0.000000 -0.000000 -0.000000 -0.000000-0.000000 4.402944 -11.362870 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 6.959926-0.000000 -0.000000 -0.000000 -0.000000 -0.000000 1.880885 -2.467393 -0.000000 -0.000000-0.000000 -0.000000 -0.000000 -0.000000-0.000000 -0.000000 -0.000000 -7.651113 2.251975 -0.000000 -0.000000 -0.000000 5.399139-0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -2.946815 2.489025-0.000000 -0.000000 -0.000000 -0.000000-0.000000 -0.000000 -0.000000 2.251975 -14.941622 2.314963 -0.000000 -0.000000 10.374684-0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 2.489025 -4.5556971.136994 -0.000000
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 潮流 计算 matlab 牛顿 拉夫逊法法