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

%*************************电力系统直角坐标系下的牛顿拉夫逊法潮流计算**********clearclcloadE:\data\IEEE014_Node.txtNode=IEEE014_Node;weishu=size(Node);nnum=weishu(1,1);%节点总数loadE:\data\IEEE014_Branch.txtbranch=IEEE014_Branch;bwei=size(branch);bnum=bwei(1,1);%支路总数Y=(zeros(nnum));Sj=100;%********************************节点导纳矩阵*******************************form=1:bnum;s=branch(m,1);%首节点e=branch(m,2);%末节点R=branch(m,3);%支路电阻X=branch(m,4);%支路电抗B=branch(m,5);%支路对地电纳k=branch(m,6);ifk==0%无变压器支路情形Y(s,e)=-1/(R+j*X);%互导纳Y(e,s)=Y(s,e);endifk~=0%有变压器支路情形Y(s,e)=-(1/((R+j*X)*k));Y(e,s)=Y(s,e);Y(s,s)=-(1-k)/((R+j*X)*k^2);Y(e,e)=-(k-1)/((R+j*X)*k);%对地导纳endY(s,s)=Y(s,s)-j*B/2;Y(e,e)=Y(e,e)-j*B/2;%自导纳的计算情形endfort=1:nnum;Y(t,t)=-sum(Y(t,:))+Node(t,12)+j*Node(t,13);%求支路自导纳endG=real(Y);%电导B=imag(Y);%电纳%******************节点分类**************************************pq=0;pv=0;blancenode=0;pqnode=zeros(1,nnum);pvnode=zeros(1,nnum);form=1:nnum;ifNode(m,2)==3blancenode=m;%平衡节点编号elseifNode(m,2)==0pq=pq+1;pqnode(1,pq)=m;%PQ节点编号elseifNode(m,2)==2pv=pv+1;pvnode(1,pv)=m;%PV节点编号endendendend%*****************************设置电压初值**********************************Uoriginal=zeros(1,nnum);%对各节点电压矩阵初始化forn=1:nnumUoriginal(1,n)=Node(n,9);%对各点电压赋初值ifNode(n,9)==0;Uoriginal(1,n)=1;%该节点为非PV节点时，将电压值赋为1endendPresion=input('请输入误差精度要求：Presion=');disp('该电力系统节点数:');disp(nnum);xiumax=0.1;counter=0;whilexiumaxPresion%****************************计算不平衡量***********************************e=real(Uoriginal);%取初始电压的实部f=imag(Uoriginal);%取初始电压的虚部deta=zeros(2*pq+2*pv,1);%构造储存功率变化量的列矩阵n=1;form=1:pq;Pi=0;Qi=0;fort=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节点的负荷有功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节点的负荷无功endS1(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),8)-Node(pqnode(1,m),6))/Sj-Qi;%计算该PQ节点的实际无功功率deta(n,1)=Q;%在该列向量中储存无功功率n=n+1;endform=1:pv;Pv=0;Qv=0;fort=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(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));endS1(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)^2-Ui;%计算电压变化量deta(n,1)=U;%储存该电压变化量n=n+1;enddeta;cerate=zeros(pq+pv,1);fork=1:pqcerate(k,1)=pqnode(1,k);endforv=1:pvcerate(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;form=1:pq;%m表示雅克比矩阵中pq节点的行数foru=1:pq+pv;%u表示雅克比矩阵中pq节点的列数t=cerate(u,1);%t为中间变量，用来标记雅克比矩阵中指定元素的个数ifpqnode(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;elseifpqnode(1,m)==t%对角线元素时的情况I=0;forg=1:nnumI=Y(t,g)*Uoriginal(1,g)+I;%计算节点的注入电流endaii=real(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;endendJacob(n,k)=H;k=k+1;Jacob(n,k)=N;k=k-1;n=n+1;Jacob(n,k)=J;k=k+1;Jacob(n,k)=L;n=n-1;k=k+1;%按照雅克比矩阵的排列规则排列pq节点的雅克比元素endk=1;n=2*m+1;%将光标定位于下一个待排列PQ节点元素的第一个位置endn=2*pq+1;k=1;%定位于PV节点的第一个位置处form=1:pv;foru=1:pq+pv;t=cerate(u,1);%t为中间变量，用来标记雅克比矩阵中指定元素的位置ifpvnode(1,m)~=t%非对角线元素情况H=G(pvnode(1,m),t)*f(1,pvnode(1,m))-B(pvnode(1,m),t)*e(1,pvnode(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;endifpvnode(1,m)==t%对角线元素情况I=0;forg=1:nnumI=Y(t,g)*Uoriginal(1,g)+I;%计算PV节点的注入电流endaii=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(1,m))+aii;R=2*f(1,pvnode(1,m));S=2*e(1,pvnode(1,m));endJacob(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节点的雅克比元素endk=1;n=n+2;%定位于下一个待排列PV节点的雅克比元素第一个位置end%*************************电压变化量化的计算与存储************************************Detau=inv(Jacob)*deta;%构建电压的变化量的列向量f=zeros(1,nnum);%给电压实部赋初值0e=zeros(1,nnum);%给电压虚部赋初值0forp=1:(pq+pv);f(1,cerate(p,1))=j*Detau(2*p-1,1);%将电压变量的奇数行赋值给fe(1,cerate(p,1))=Detau(2*p,1);%将电压变量的偶数行赋值给eendt=e+f;xiumax=abs(Detau(1,1));%将电压变化量的第一个元素赋值给最大允许误差forn=2:2*nnum-2;ifabs(Detau(n,1))xiumaxxiumax=abs(Detau(n,1));%找出最大的电压误差endendUoriginal=Uoriginal+t;%迭代修正后的电压值counter=1+counter;%统计迭代次数enddisp('迭代次数counter:');disp(counter);%**************************平衡节点功率及显示**********************************m=blancenode;t=0;forn=1:nnum;t=t+(G(m,n)-j*B(m,n))*(real(Uoriginal(1,n))-j*imag(Uoriginal(1,n)));endS1(1,m)=Uoriginal(1,m)*t;%