基于MATLAB的潮流计算源程序代码

1、%*电力系统直角坐标系下的牛顿拉夫逊法潮流计算*clearclcload E:dataIEEE014_Node.txtNode=IEEE014_Node;weishu=size(Node);nnum=weishu(1,1); %节点总数 load E:dataIEEE014_Branch.txtbranch=IEEE014_Branch;bwei=size(branch); bnum=bwei(1,1); %支路总数Y=(zeros(nnum);Sj=100;%*节点导纳矩阵*for m=1:bnum; s=branch(m,1); %首节点 e=branch(m,2); %末节点 R=bra

2、nch(m,3); %支路电阻 X=branch(m,4); %支路电抗 B=branch(m,5); %支路对地电纳 k=branch(m,6); if k=0 %无变压器支路情形 Y(s,e)=-1/(R+j*X); %互导纳 Y(e,s)=Y(s,e); end if k=0 %有变压器支路情形 Y(s,e)=-(1/(R+j*X)*k); Y(e,s)=Y(s,e); Y(s,s)=-(1-k)/(R+j*X)*k2); Y(e,e)=-(k-1)/(R+j*X)*k); %对地导纳 end Y(s,s)=Y(s,s)-j*B/2; Y(e,e)=Y(e,e)-j*B/2; %自导纳的

3、计算情形endfor t=1:nnum; Y(t,t)=-sum(Y(t,:)+Node(t,12)+j*Node(t,13); %求支路自导纳end G=real(Y); %电导 B=imag(Y); %电纳%*节点分类*pq=0; pv=0; blancenode=0;pqnode=zeros(1,nnum);pvnode=zeros(1,nnum);for m=1:nnum; if Node(m,2)=3 blancenode=m; %平衡节点编号 else if Node(m,2)=0 pq=pq+1; pqnode(1,pq)=m; %PQ节点编号 else if Node(m,2)

4、=2 pv=pv+1; pvnode(1,pv)=m; %PV节点编号 end end endend%*设置电压初值*Uoriginal=zeros(1,nnum); %对各节点电压矩阵初始化for n=1:nnum Uoriginal(1,n)=Node(n,9); %对各点电压赋初值 if Node(n,9)=0; Uoriginal(1,n)=1; %该节点为非PV节点时，将电压值赋为1 endendPresion=input(请输入误差精度要求：Presion=);disp(该电力系统节点数:);disp(nnum);xiumax=0.1;counter=0;while xiumaxP

5、resion %*计算不平衡量*e=real(Uoriginal); %取初始电压的实部f=imag(Uoriginal); %取初始电压的虚部deta=zeros(2*pq+2*pv,1); %构造储存功率变化量的列矩阵n=1; for m=1:pq; Pi=0;Qi=0; for t=1:nnum; Pi=Pi+e(1,pqnode(1,m)*(G(pqnode(1,m),t)*e(1,t)-B(pqnode(1,m),t)*f(1,t)+f(1,pqnode(1,m)*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t)*e(1,t);%计算该PQ节点的负荷

6、有功 Qi=Qi+f(1,pqnode(1,m)*(G(pqnode(1,m),t)*e(1,t)-B(pqnode(1,m),t)*f(1,t)-e(1,pqnode(1,m)*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t)*e(1,t);%计算该PQ节点的负荷无功 end S1(1,pqnode(1,m)=Pi+j*Qi; P=(Node(pqnode(1,m),7)-Node(pqnode(1,m),5)/Sj-Pi;%计算该PQ节点的实际有功功率 deta(n,1)=P;%在该列向量中储存有功功率 n=n+1; Q=(Node(pqnode(1,m)

7、,8)-Node(pqnode(1,m),6)/Sj-Qi;%计算该PQ节点的实际无功功率 deta(n,1)=Q;%在该列向量中储存无功功率 n=n+1; endfor m=1:pv; Pv=0; Qv=0; for t=1:nnum; Pv=Pv+e(1,pvnode(1,m)*(G(pvnode(1,m),t)*e(1,t)-B(pvnode(1,m),t)*f(1,t)+f(1,pvnode(1,m)*(G(pvnode(1,m),t)*f(1,t)+B(pvnode(1,m),t)*e(1,t);%计算该PV节点的负荷有功 Ui=e(1,pvnode(1,m)2+f(1,pvnode

8、(1,m)2;%计算该节点的负荷电压值 Qv=Qv+f(1,pqnode(1,m)*(G(pqnode(1,m),t)*e(1,t)-B(pqnode(1,m),t)*f(1,t)-e(1,pqnode(1,m)*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t)*e(1,t); end S1(1,pvnode(1,m)=Pv+j*Qv; P=(Node(pvnode(1,m),7)-Node(pvnode(1,m),5)/Sj-Pv; %计算该节点的实际有功功率 deta(n,1)=P; %储存该有功功率 n=n+1; U=Node(pvnode(1,m),3

9、)2-Ui; %计算电压变化量 deta(n,1)=U; %储存该电压变化量 n=n+1;enddeta;cerate=zeros(pq+pv,1); for k=1:pq cerate(k,1)=pqnode(1,k);end for v=1:pv cerate(pq+v,1)=pvnode(1,v);end%*雅克比矩阵* Jacob=ones(2*nnum-2); L=0;J=0;H=0;N=0; R=0;S=0;n=1;k=1; for m=1:pq; %m表示雅克比矩阵中pq节点的行数 for u=1:pq+pv; %u表示雅克比矩阵中pq节点的列数 t=cerate(u,1); %

10、t为中间变量，用来标记雅克比矩阵中指定元素的个数 if pqnode(1,m)=t %非对角元素的情况 H=G(pqnode(1,m),t)*f(1,pqnode(1,m)-B(pqnode(1,m),t)*e(1,pqnode(1,m); N=G(pqnode(1,m),t)*e(1,pqnode(1,m)+B(pqnode(1,m),t)*f(1,pqnode(1,m); L=H; J=-N; else if pqnode(1,m)=t %对角线元素时的情况 I=0; for g=1:nnum I=Y(t,g)*Uoriginal(1,g)+I; %计算节点的注入电流 end aii=re

11、al(I); bii=imag(I); H=-B(t,t)*e(1,pqnode(1,m)+G(t,t)*f(1,pqnode(1,m)+bii; N=G(t,t)*e(1,pqnode(1,m)+B(t,t)*f(1,pqnode(1,m)+aii; L=-B(t,t)*e(1,pqnode(1,m)+G(t,t)*f(1,pqnode(1,m)-bii; J=-G(t,t)*e(1,pqnode(1,m)-B(t,t)*f(1,pqnode(1,m)+aii; end end Jacob(n,k)=H; k=k+1; Jacob(n,k)=N; k=k-1;n=n+1; Jacob(n,k

12、)=J; k=k+1; Jacob(n,k)=L; n=n-1; k=k+1; %按照雅克比矩阵的排列规则排列pq节点的雅克比元素 end k=1; n=2*m+1; %将光标定位于下一个待排列PQ节点元素的第一个位置endn=2*pq+1; k=1; %定位于PV节点的第一个位置处 for m=1:pv; for u=1:pq+pv; t=cerate(u,1); %t为中间变量，用来标记雅克比矩阵中指定元素的位置 if pvnode(1,m)=t %非对角线元素情况 H=G(pvnode(1,m),t)*f(1,pvnode(1,m)-B(pvnode(1,m),t)*e(1,pvnode

13、(1,m); N=G(pvnode(1,m),t)*e(1,pvnode(1,m)+B(pvnode(1,m),t)*f(1,pvnode(1,m); R=0; S=0; end if pvnode(1,m)=t %对角线元素情况 I=0; for g=1:nnum I=Y(t,g)*Uoriginal(1,g)+I; %计算PV节点的注入电流 end aii=real(I); bii=imag(I); H=-B(t,t)*e(1,pvnode(1,m)+G(t,t)*f(1,pvnode(1,m)+bii; N=G(t,t)*e(1,pvnode(1,m)+B(t,t)*f(1,pvnode

14、(1,m)+aii; R=2*f(1,pvnode(1,m); S=2*e(1,pvnode(1,m); end Jacob(n,k)=H; k=k+1; Jacob(n,k)=N; k=k-1;n=n+1; Jacob(n,k)=R; k=k+1; Jacob(n,k)=S; n=n-1;k=k+1; %按照雅克比矩阵的排列规则排列PV节点的雅克比元素 end k=1; n=n+2; %定位于下一个待排列PV节点的雅克比元素第一个位置 end %*电压变化量化的计算与存储*Detau=inv(Jacob)*deta; %构建电压的变化量的列向量f=zeros(1,nnum); %给电压实部赋

15、初值0 e=zeros(1,nnum); %给电压虚部赋初值0 for p=1:(pq+pv); f(1,cerate(p,1)=j*Detau(2*p-1,1);%将电压变量的奇数行赋值给f e(1,cerate(p,1)=Detau(2*p,1); %将电压变量的偶数行赋值给eendt=e+f;xiumax=abs(Detau(1,1); %将电压变化量的第一个元素赋值给最大允许误差for n=2:2*nnum-2; if abs(Detau(n,1)xiumax xiumax=abs(Detau(n,1); %找出最大的电压误差 endendUoriginal=Uoriginal+t;

16、%迭代修正后的电压值counter=1+counter; %统计迭代次数enddisp(迭代次数counter:);disp(counter);%*平衡节点功率及显示*m=blancenode; t=0; for n=1:nnum; t=t+(G(m,n)-j*B(m,n)*(real(Uoriginal(1,n)-j*imag(Uoriginal(1,n); end S1(1,m)=Uoriginal(1,m)*t;%*直角坐标下各节点电压及显示*U=zeros(1,nnum);for n=1:nnum Ui(n,1)=Node(n,1); U(1,n)=real(Uoriginal(1,n

17、)+i*imag(Uoriginal(1,n); %将电压值由极坐标转化为直角坐标形式 Ui(n,2)=U(1,n); Ui(n,3)=S1(1,n); end disp(各节点电压直角坐标形式及节点注入功率:);disp( 节点号 节点电压值 节点注入功率);disp(Ui);disp(修正电压的最大误差:)disp(xiumax);%*功率损耗* for m=1:bnum %支路功率及损耗 startnode=branch(m,1); endnode=branch(m,2); %终止节点 y=sum(Y,2); Sij=Uoriginal(1,startnode)*(conj(Uoriginal(1,startnode)*conj(y(startnode,1)+(conj(Uoriginal(1,startnode)-conj(Uoriginal(1,endnode)*conj(-Y(startnode,endnode); Sji=Uoriginal(1,endnode)*(conj(Uoriginal(1,endnode)*conj(y(endnode