选用牛顿-拉佛森方法,利用matlab 软件计算基于PQ 节点情况下的潮流计算。
一.所用公式112222[()()][()()]()j ni i i ij j ij j i ij j ij j j j n i i i ij j ij j i ij j ij j j i i i iP P e G e B f f G f B e Q Q f G e B f e G f B e U U e f ====⎧∆=--++⎪⎪⎪⎪∆=---+⎨⎪⎪∆=-+⎪⎪⎩∑∑i j ≠2200i ij ij i ij iii ij ij i ij ii i ijij i ij i ij ii ij ij i ij i iji i iji i ij i P H B e G f f P N G e B f e Q J G e B f N f Q L B e G f H e U R f U S e ∂⎧==-+⎪∂⎪⎪∂==+⎪∂⎪⎪∂==--=-⎪∂⎪⎨∂⎪==-+=∂⎪⎪∂⎪==⎪∂⎪∂⎪==⎪∂⎩i j=2222i ii ij i ij i iiiiii ij i ij i iii iii ij i ij i ii ii ii ij i ij i ii i i iii i i ii ii P H B e G f b f P N G e B f a e Q J G e B f a f Q L B e G f b e U R f f U S e e ∂⎧==-++⎪∂⎪⎪∂==++⎪∂⎪⎪∂==--+⎪∂⎪⎨∂⎪==-+-∂⎪⎪∂⎪==⎪∂⎪∂⎪==⎪∂⎩其中11()()j nii ii i ii i ij j ij j j j i j niiii i ii i ij j ij jj j i a G e B f G e B f b G f B e G f B e ==≠==≠⎧=-+-⎪⎪⎪⎨⎪=+++⎪⎪⎩∑∑二、程序流程图开始形成节点导纳矩阵输入原始数据设节点电压(0)(0)i ie f,i=1,2…,n,i≠s置迭代次数0k=置节点号i=1计算雅克比矩阵元素计算PQ节点的()kiP∆,()kiQ∆,PV节点的()kiP∆,()2kiU∆求解修正方程式,得()kie∆,()kif∆雅克比矩阵是否已全部形迭代次数k=k+1i=i+1计算各节点电压的新值:牛顿-拉佛森例题中的简单模型系统一、系统单线图图1 简单模型系统二、系统参数节点矩阵%(bus#)(volt)(ang)(p)(q)(bus type)bus=[1 1.00 0.00 -1.60 -0.80 1;2 1.00 0.00 -2.00 -1.00 1;3 1.00 0.00 -3.70 -1.30 1;4 1.05 0.00 5.00 0.00 2;5 1.05 0.00 0.00 0.00 3];线路矩阵%b#1 b#2 (R)(X)(G)(B)(K)line=[1 2 0.04 0.25 0 0.25 0;1 3 0.10 0.35 0 0.00 0;2 3 0.08 0.30 0 0.25 0;5 3 0.00 0.03 0 0.00 1.05;4 2 0.00 0.015 0 0.00 1.05] ;三、计算结果:牛顿-拉夫逊法潮流计算结果节点计算结果:n节点节点电压节点相角(角度)节点注入功率1 0.862150 -4.778511 -1.600000 + j -0.8000002 1.077916 17.853530 -2.000000 + j -1.0000003 1.036411 -4.281930 -3.700000 + j -1.3000004 1.050000 21.843319 5.000000 + j 1.8130845 1.050000 0.000000 2.579427 + j 2.299402n线路计算结果:n节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)1 2 -1.466181 + j -0.409076 1.584546 + j 0.672556 0.118365 + j0.2634801 3 -0.133819 + j -0.390924 0.156788 + j 0.471315 0.022969 + j 0.0803912 3 1.415454 + j -0.244333 -1.277360 + j 0.203170 0.138093 + j -0.041163 5 3 2.579427 + j 2.299402 -2.579427 + j -1.974485 0.000000 + j 0.3249174 2 5.000000 + j 1.813084 -5.000000 + j -1.428223 0.000000 + j 0.384861导纳矩阵:Y =[ 1.3787 - 6.2917i -0.6240 + 3.9002i -0.7547 + 2.6415i 0 0-0.6240 + 3.9002i 1.4539 -66.9808i -0.8299 + 3.1120i 0 +63.4921i 0-0.7547 + 2.6415i -0.8299 + 3.1120i 1.5846 -35.7379i 0 0 +31.7460i0 0 +63.4921i 0 0 -66.6667i 00 0 0 +31.7460i 0 0 -33.3333i ]图2. Matlab运行结果结果分析:此程序的运行结果和试验程序给出的结果是一致的。
说明程序无误,但在精确度上有微小差异,这主要是和导纳矩阵的精确度以及显示精度有关。
心得:本程序分模块进行,先是排序,再是求导纳阵,然后求雅阁比,再进行迭代运算,程序本身很简洁明了,运行的时候只需要在matlab里输入main就行了,然后打开BUS和line所在的.m文件,结果就会自动存在result文件中了,通过编写牛顿拉夫逊法matlab潮流计算程序复习了潮流计算的知识,也实现了计算机算法附录:实验源程序:Main函数:clear[dfile,pathname]=uigetfile('*.m','Select Data File');if pathname == 0error(' you must select a valid data file')elselfile =length(dfile);% strip off .meval(dfile(1:lfile-2));end;global bus;global 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,'计算结果');fclose(myf); % 在当前目录下生成“Result.m”文件,写入节点导纳矩阵format longEPS = 1.0e-10; % 设定误差精度for t = 1:100 % 开始迭代计算,设定最大迭代次数为100,以便不收敛情况下及时跳出[dP,dQ] = dPQ(Y,bus,nPQ,nPV); % 计算功率偏差dP和dQ的子程序J = Jac(bus,Y,nPQ); % 计算雅克比矩阵的子程序UD = zeros(nPQ,nPQ);for i = 1:nPQUD(i,i) = bus(i,2); % 生成电压对角矩阵 endenddAngU = 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');for i = 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');for i = 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”写入节点计算结果和线路计算结果程序结束"Num.m" 作用为对节点重排序,并修改相应的线路数据function [bus,line,nPQ,nPV,nodenum] = Num(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);nSW = 0; % number of swing bus counternPV = 0; % number of PV bus counternPQ = 0; % number of PQ bus counterfor i = 1:nb, % nb为总节点数type= bus(i,6);if type == 3,nSW = nSW + 1; % increment swing bus counterSW(nSW,:)=bus(i,:);elseif type == 2,nPV = nPV +1; % increment PV bus counterPV(nPV,:)=bus(i,:);elsenPQ = nPQ + 1; % increment PQ bus counterPQ(nPQ,:)=bus(i,:);endend;bus=[PQ;PV;SW];newbus=[1:nb]';nodenum=[newbus bus(:,1)];bus(:,1)=newbus;for i=1:nlfor j=1:2for k=1:nbif line(i,j)==nodenum(k,2)line(i,j)=nodenum(k,1);breakendendendend"y.m" 作用为计算节点导纳矩阵function Y = y(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1); %读入线路参数J=line(k,2);Zt=line(k,3)+j*line(k,4);Yt=1/Zt;Ym=line(k,5)+j*line(k,6);K=line(k,7);if (K==0)&(J~=0) % 普通线路: K=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=0,J=0,R=X=0;Y(I,I)=Y(I,I)+Ym;endif K>0 % 变压器线路: 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);endif K<0 % 变压器线路: 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"dPQ.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 平衡节点不参与计算for i = 1:n-1for j = 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)));if i<nPQ+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"Jac.m" 作用为计算雅克比矩阵function J = 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 = [H N; K L],并初始化Qi = zeros(nb-1,1);Pi = zeros(nb-1,1);for i = 1:nb-1for j = 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)+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)));endend % 初始化并计算Qi和Pifor i = 1:nb-1for j = 1:nb-1if i~=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)));elseH(i,j)=Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);end % 分别计算H矩阵的对角及非对角元素if j < nPQ+1if i~=jN(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)));elseN(i,j)=-Pi(i,1)-real(Y(i,j))*((bus(i,2))^2);endend % 分别计算N矩阵的对角及非对角元素if i < nPQ+1if i~=jK(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*s in(bus(i,3)-bus(j,3)));elseK(i,j)=-Pi(i,1)+real(Y(i,j))*((bus(i,2))^2);end % 分别计算K矩阵的对角及非对角元素if j < nPQ+1if i~=jL(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)));elseL(i,j)=-Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2); endend % 分别计算L矩阵的对角及非对角元素endendendJ = [H N; K L]; % 生成雅克比矩阵"PQ.m" 作用为计算每个节点的功率注入function bus = PQ(bus,Y,nPQ,nPV)n = nPQ+nPV+1; % n为总节点数for i = nPQ+1:n-1for j = 1:nbus(i,5)=bus(i,5)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-ima g(Y(i,j))*cos(bus(i,3)-bus(j,3)));endend % 利用公式计算PV节点的无功注入for j =1:nbus(n,4)=bus(n,4)+bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+ima g(Y(n,j))*sin(bus(n,3)-bus(j,3)));bus(n,5)=bus(n,5)+bus(n,2)*bus(j,2)*(real(Y(n,j))*sin(bus(n,3)-bus(j,3))-ima g(Y(n,j))*cos(bus(n,3)-bus(j,3)));end % 计算平衡节点的无功及有功注入"ReNum.m" 作用为对节点和线路数据恢复编号function [bus,line] = ReNum(bus,line,nodenum)[nb,mb]=size(bus);[nl,ml]=size(line);bus_temp = zeros(nb,mb); % bus_temp矩阵用于临时存放bus矩阵的数据k = 1;for i = 1 :nbfor j = 1 : nbif bus(j,1) == kbus_temp(k,:) = bus(j,:);k = k + 1;endendend % 利用bus矩阵的首列编号重新对bus矩阵排序并存入bus_temp矩阵中bus = bus_temp; % 重新赋值回bus,完成bus矩阵的重新编号for i=1:nlfor j=1:2for k=1:nbif line(i,j)==nodenum(k,1)line(i,j)=nodenum(k,2);breakendendendend % 恢复line的编号"YtYm_.m" 作用为计算线路的等效Yt和Ym,以计算线路潮流function YtYm = YtYm_(line)[nl,ml]=size(line);YtYm = zeros(nl,5); % 对YtYm矩阵赋初值0YtYm(:,1:2) = line(:,1:2); % 矩阵前两列为线路两段节点编号,后三列分别为线路等效Yt,i侧的等效Ym,j侧的等效Ymj = sqrt(-1);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+j*line(k,4);if Zt~=0Yt=1/Zt;endYm=line(k,5)+j*line(k,6);K=line(k,7);if (K==0)&(J~=0) % 普通线路: K=0YtYm(k,3) = Yt;YtYm(k,4) = Ym;YtYm(k,5) = Ym;endif (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0YtYm(k,4) = Ym;endif K>0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧YtYm(k,3) = Yt/K;YtYm(k,4) = Ym+Yt*(K-1)/K;YtYm(k,5) = Yt*(1-K)/K/K;endif K<0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧YtYm(k,3) = -Yt*K;YtYm(k,4) = Ym+Yt*(1+K);YtYm(k,5) = Yt*(K^2+K);endend"bus_res_.m" 计算并返回节点数据结果function bus_res = bus_res_(bus)[nb,mb]=size(bus);bus_res = zeros(nb,4); % bus_res矩阵储存着节点计算结果bus_res(:,1:2) = bus(:,1:2);bus_res(:,3) = bus(:,3) *180 / pi; % 相角采用角度制bus_res(:,4) = bus(:,4) + (sqrt(-1))*bus(:,5); % 注入功率"S_res_.m" 计算并返回线路潮流function S_res = S_res_(bus,line,YtYm)[nl,ml]=size(line);S_res = zeros(nl,5); % S_res矩阵储存着线路潮流计算结果S_res(:,1:2) = line(:,1:2); % 前两列为节点编号for k=1:nlI = S_res(k,1);J = S_res(k,2);if (J~=0)&(I~=0)S_res(k,3)=bus(I,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4)))-bus(I,2)*bus(J,2)*(c os(bus(I,3))+j*sin(bus(I,3)))*(conj(cos(bus(J,3))+j*sin(bus(J,3))))*conj(YtY m(k,3));S_res(k,4)=bus(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))-bus(I,2)*bus(J,2)*(c onj(cos(bus(I,3))+j*sin(bus(I,3))))*(cos(bus(J,3))+j*sin(bus(J,3)))*conj(YtY m(k,3));S_res(k,5)=S_res(k,3) + S_res(k,4); % 利用公式计算非接地支路的潮流 else if(J==0 )S_res(k,3)=bus(I,2)^2*conj(YtYm(k,4));S_res(k,5)=S_res(k,3)+S_res(k,4);elseS_res(k,4)=bus(J,2)^2*conj(YtYm(k,5));S_res(k,5)=S_res(k,3)+S_res(k,4); % 利用公式计算接地支路的潮流endendendend。