Matlab-牛拉法代码

牛顿-拉夫逊法潮流计算的Matlab实现(2009-09-2613:23:26)转载标签：电力系统潮流计算牛顿-拉夫逊牛顿拉夫逊matlabnewton-raphson写了近一个礼拜的程序，藏着可惜，拿出来给需要的人转载请注明出处有不少朋友反映程序报错：???Outputargumentbus(andmaybeothers)notassignedduringcalltoE:\MAT\best\OpDF_.mOpDF_.是因为输入数据的格式问题输入数据中是包含bus,line的，并非只有两个矩阵，如：bus=[11.000.00-0.30-0.181;21.000.00-0.55-0.131;31.100.000.500.002;41.050.000.000.003];line=[120.100.400.00.015280;420.080.400.00.014130;140.120.500.00.01920;310.000.30.00.0-1.1];所有蓝色文字都为数据文件中的内容---------------------分---------------------割---------------------线---------------------【注意】一、此程序只适用于求解节点电压以极坐标形式表示的潮流方程，没有考虑节点优化编号二、程序在Matlab6.5上测试通过，应该适用于目前其所有后续版本三、程序主要通过文件方式输入输出（同时程序也返回结果向量），对输入文件格式有严格要求，具体如下：1.输入文件可以直接在Matlab中新建m文件编写，也内容可以以文本方式编写，但最后必须改后缀名成为“.m”文件。文件名第一个字符必须是字母，后面可以跟字母、数字和下划线的任何组合，但不能和已有文件和函数冲突，不能含中文。2.输入文件内容格式：内容中包括bus（节点数据）和line（线路数据）两个矩阵其中bus的格式为：bus=[节点编号节点电压节点相角（弧度制）有功注入无功注入节点类型];line的格式为：line=[节点i节点j线路电阻电抗电导电纳变压器变比（普通线路为零）];四、程序分主程序和各子程序，需全部放入Matlab当前目录下调用，计算开始后会在当前目录下生成文件Result.m保存运算结果。---------------------分---------------------割---------------------线---------------------以下为程序代码（红色部分，“%”后为注释）：主程序PowerFlow_NR.mfunction[bus_res,S_res]=PowerFlow_NR_2%牛顿-拉夫逊法解潮流方程的主程序%%%%%%%%%%%%%%%%%%%%%bylongdinhohe%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[bus,line]=OpDF_;%打开数据文件的子程序，返回bus（节点数据）和line（线路数据）回主程序[nb,mb]=size(bus);[nl,ml]=size(line);%计算bus和line矩阵的行数和列数[bus,line,nPQ,nPV,nodenum]=Num_(bus,line);%对节点重新排序的子程序Y=y_(bus,line);%计算节点导纳矩阵的子程序myf=fopen('Result.m','w');fprintf(myf,'---------------bylongdinhohe------------------------------\n\n');fclose(myf);%在当前目录下生成“Result.m”文件，写入节点导纳矩阵formatlongEPS=1.0e-10;%设定误差精度fort=1:100%开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出[dP,dQ]=dPQ_(Y,bus,nPQ,nPV);%计算功率偏差dP和dQ的子程序J=Jac_(bus,Y,nPQ);%计算雅克比矩阵的子程序UD=zeros(nPQ,nPQ);fori=1:nPQUD(i,i)=bus(i,2);%生成电压对角矩阵enddAngU=J\[dP;dQ];dAng=dAngU(1:nb-1,1);%计算相角修正量dU=UD*(dAngU(nb:nb+nPQ-1,1));%计算电压修正量bus(1:nPQ,2)=bus(1:nPQ,2)-dU;%修正电压bus(1:nb-1,3)=bus(1:nb-1,3)-dAng;%修正相角if(max(abs(dU))EPS)&(max(abs(dAng))EPS)breakend%判断是否满足精度误差，如满足则跳出，否则返回继续迭代endbus=PQ_(bus,Y,nPQ,nPV);%计算每个节点的有功和无功注入的子程序[bus,line]=ReNum_(bus,line,nodenum);%对节点恢复编号的子程序YtYm=YtYm_(line);%计算线路的等效Yt和Ym的子程序，以计算线路潮流bus_res=bus_res_(bus);%计算节点数据结果的子程序S_res=S_res_(bus,line,YtYm);%计算线路潮流的子程序myf=fopen('Result.m','a');fprintf(myf,'---------------牛顿－拉夫逊法潮流计算结果----------\n\n节点计算结果：\n节点节点电压节点相角（角度）节点注入功率\n');fori=1:nbfprintf(myf,'%2.0f',bus_res(i,1));fprintf(myf,'%10.6f',bus_res(i,2));fprintf(myf,'%10.6f',bus_res(i,3));fprintf(myf,'%10.6f+j%10.6f\n',real(bus_res(i,4)),imag(bus_res(i,4)));endfprintf(myf,'\n线路计算结果：\n节点I节点J线路功率S(I,J)线路功率S(J,I)线路损耗dS(I,J)\n');fori=1:nlfprintf(myf,'%2.0f',S_res(i,1));fprintf(myf,'%2.0f',S_res(i,2));fprintf(myf,'%10.6f+j%10.6f',real(S_res(i,3)),imag(S_res(i,3)));fprintf(myf,'%10.6f+j%10.6f',real(S_res(i,4)),imag(S_res(i,4)));fprintf(myf,'%10.6f+j%10.6f\n',real(S_res(i,5)),imag(S_res(i,5)));endfclose(myf);%迭代结束后继续在“Result.m”写入节点计算结果和线路计算结果程序结束子程序1OpDF_.m作用为打开数据文件function[bus,line]=OpDF_[d]=uigetfile('*.m','SelectDataFile');%数据文件类型为m文件，窗口标题为“SelectDataFile”ifpathname==0error('youmustselectavaliddatafile')%如果没有选择有效文件，则出现错误提示elsel(dfile);%去除文件后缀(d));%打开文件end子程序2Num.m作用为对节点重排序，并修改相应的线路数据function[bus,line,nPQ,nPV,nodenum]=Num_(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);nSW=0;%nSW为平衡节点个数nPV=0;%nPV为PV节点个数nPQ=0;%nPQ为PQ节点个数fori=1:nb,%nb为总节点数type=bus(i,6);iftype==3,nSW=nSW+1;SW(nSW,:)=bus(i,:);%计算并储存平衡节点elseiftype==2,nPV=nPV+1;PV(nPV,:)=bus(i,:);%计算并储存PV节点elsenPQ=nPQ+1;PQ(nPQ,:)=bus(i,:);%计算并储存PQ节点endendbus=[PQ;PV;SW];%对bus矩阵按PQ、PV、平衡节点的顺序重新排序nodenum=[[1:nb]'bus(:,1)];%生成新旧节点对照表fori=1:nlforj=1:2fork=1:nbifline(i,j)==nodenum(k,2)line(i,j)=nodenum(k,1);breakendendendend%按排序以后的节点顺序对line矩阵重新编号子程序3y_.m作用为计算节点导纳矩阵functionY=y_(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);%对导纳矩阵赋初值0fork=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+j*line(k,4);%读入线路参数ifZt~=0Yt=1/Zt;%接地支路不计算YtendYm=line(k,5)+j*line(k,6);K=line(k,7);if(K==0)&(J~=0)%普通线路:K=0Y(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=0,J=0,R=X=0Y(I,I)=Y(I,I)+Ym;endifK0%变压器线路:Zt和Ym为折算到i侧的值,K在j侧Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt/K/K;Y(I,J)=Y(I,J)-Yt/K;Y(J,I)=Y(I,J);endifK0%变压器线路:Zt和Ym为折算到K侧的值,K在i侧Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+K*K*Yt;Y(I,J)=Y(I,J)+K*Yt;Y(J,I)=Y(I,J);endend子程序4dPQ_.m作用为计算功率偏差function[dP,dQ]=dPQ_(Y,bus,nPQ,nPV)%nPQ、nPV为相应节点个数n=nPQ+nPV+1;%总节点个数dP=bus(1:n-1,4);dQ=bus(1:nPQ,5);%对dP和dQ赋初值PV节点不需计算dQ平衡节点不参与计算fori=1:n-1forj=1:ndP(i,1)=dP(i,1)-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)));ifinPQ+1dQ(i,1)=dQ(i,1)-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)));endendend%利用循环计算求取dP和dQ子程序5Jac_.m作用为计算雅克比矩阵functionJ=Jac_(bus,Y,nPQ)[nb,mb]=size(bus);H=zeros(nb-1,nb-1);N=zeros(nb-1,nPQ);K=zeros(nPQ,nb-1);L=zeros(nPQ,nPQ);%将雅克比矩阵分块，即：J=[HN;KL]，并初始化Qi=zeros(nb-1,1);Pi=zeros(nb-1,1);fori=1:nb-1forj=1:nbPi(i,1)=Pi(i,1)+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)));Qi(i,1)=Qi(i,1)+