现代电力系统分析.doc

(完整word)现代电力系统分析 选用牛顿拉佛森方法，利用matlab软件计算基于PQ节点情况下的潮流计算。一所用公式 其中 二、程序流程图开始形成节点导纳矩阵输入原始数据设节点电压，i=1,2,n,is置迭代次数置节点号i=1 计算雅克比矩阵元素 计算节点的，节点的，求解修正方程式，得，雅克比矩阵是否已全部形成？计算平衡节点及PV节点功率求，迭代次数 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.800000 2 1。077916 17.853530 -2。000000 + j -1。000000 3 1.036411 4。281930 -3.700000 + j -1。300000 4 1.050000 21.843319 5。000000 + j 1.813084 5 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 + j 0。263480 1 3 -0.133819 + j 0.390924 0.156788 + j 0。471315 0。022969 + j 0。080391 2 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。324917 4 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。7460i 0 0 +63.4921i 0 0 66.6667i 0 0 0 0 +31。7460i 0 0 -33.3333i 图2。 Matlab运行结果 结果分析：此程序的运行结果和试验程序给出的结果是一致的。说明程序无误，但在精确度上有微小差异，这主要是和导纳矩阵的精确度以及显示精度有关。 心得：本程序分模块进行,先是排序，再是求导纳阵，然后求雅阁比，再进行迭代运算，程序本身很简洁明了，运行的时候只需要在matlab里输入main就行了，然后打开BUS和line所在的。m文件，结果就会自动存在result文件中了,通过编写牛顿拉夫逊法matlab潮流计算程序复习了潮流计算的知识，也实现了计算机算法附录：实验源程序:Main函数:cleardfile,pathname=uigetfile（*.m，Select Data File)；if pathname = 0 error（ you must select a valid data file)else lfile =length（dfile）； strip off .m eval（dfile(1：lfile2）；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 long EPS = 1.0e10； % 设定误差精度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 = JdP；dQ;dAng = dAngU(1:nb-1,1）； 计算相角修正量dU = UD（dAngU(nb:nb+nPQ1，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）break end % 判断是否满足精度误差，如满足则跳出，否则返回继续迭代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.6fn，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 + j10。6fn,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 K0 变压器线路： 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 K0 变压器线路： Zt和Ym为折算到K侧的值，K在i侧Y（I，I）=Y(I，I)+Yt+Ym；Y(J,J)=Y（J，J)+KKYt；Y（I,J)=Y（I,J)+KYt；Y(J,I）=Y(I，J);endEnddPQ.m 作用为计算功率偏差function dP，dQ =dPQ（Y,bus,nPQ,nPV） % nPQ、nPV为相应节点个数n = nPQ + nPV +1; 总节点个数dP = bus(1:n1,4);dQ = bus(1:nPQ,5）； 对dP和dQ赋初值 PV节点不需计算dQ 平衡节点不参与计算for i = 1:n1 for j = 1：n dP（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 inPQ+1 dQ（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); end endend 利用循环计算求取dP和dQ”Jac。m” 作用为计算雅克比矩阵 function J = Jac(bus,Y，nPQ）nb,mb=size（bus);H = zeros（nb1,nb-1)；N = zeros(nb-1,nPQ);K = zeros（nPQ,nb1）；L = zeros（nPQ,nPQ）; 将雅克比矩阵分块，即:J = H N； K L，并初始化Qi = zeros(nb-1,1)；Pi = zeros（nb-1，1）;for i = 1:nb-1 for j = 1：nb Pi（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和Pi for i = 1：nb-1 for j = 1:nb1 if i=j H（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）)）； else H(i，j)=Qi(i，1)+imag（Y（i,j）)*(（bus(i，2）)2）； end % 分别计算H矩阵的对角及非对角元素 if j nPQ+1 if i=j 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)； else N(i,j)=-Pi（i，1)-real(Y(i,j)*(（bus(i,2)2）； end end % 分别计算N矩阵的对角及非对角元素 if i nPQ+1 if i=j 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)； else K(i,j）=-Pi（i,1）+real（Y（i，j）(bus（i，2）)2）； end % 分别计算K矩阵的对角及非对角元素 if j nPQ+1 if i=j L(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)）; else L（i，j)=Qi(i,1)+imag（Y(i,j）)(bus（i，2)）2）; end end 分别计算L矩阵的对角及非对角元素 end endendJ = H N; K L; 生成雅克比矩阵PQ。m” 作用为计算每个节点的功率注入function bus = PQ(bus，Y,nPQ,nPV）n = nPQ+nPV+1; n为总节点数for i = nPQ+1:n1 for j = 1：n bus（i，5）=bus(i,5)+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 利用公式计算PV节点的无功注入for j =1：n bus(n,4）=bus(n，4)+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）); bus(n，5）=bus（n,5)+bus(n,2）bus(j,2）（real(Y（n，j）)*sin（bus(n,3)bus（j，3)-imag(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 ：nb for j = 1 : nb if bus(j,1） = k bus_temp（k，：） = bus（j,：)； k = k + 1; end endend 利用bus矩阵的首列编号重新对bus矩阵排序并存入bus_temp矩阵中bus = bus_temp; % 重新赋值回bus，完成bus矩阵的重新编号for i=1:nl for j=1：2 for k=1:nb if line(i,j）=nodenum(k,1） line(i，j)=nodenum（k，2）; break end end endend % 恢复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：nl I=line(k，1); J=line(k,2)； Zt=line(k，3）+j*line(k，4）; if Zt=0 Yt=1/Zt； end Ym=line（k，5）+j*line(k，6）； K=line（k,7）； if (K=0)&(J=0) 普通线路： K=0 YtYm(k，3) = Yt; YtYm（k,4) = Ym; YtYm(k,5） = Ym; end if （K=0）&(J=0） % 对地支路： K=0,J=0，R=X=0 YtYm（k,4） = Ym； end if K0 变压器线路： Zt和Ym为折算到i侧的值,K在j侧 YtYm（k,3） = Yt/K; YtYm(k，4） = Ym+Yt*(K-1）/K; YtYm（k，5) = Yt（1K)/K/K； end if K0 变压器线路: Zt和Ym为折算到K侧的值，K在i侧 YtYm(k,3) = YtK; YtYm(k，4） = Ym+Yt（1+K）； YtYm(k,5） = Yt（K2+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：nl I = 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)*(cos(bus(I，3）)+j*sin(bus(I，3）)（conj(cos（bus（J,3）)+jsin（bus(J，3)）*conj(YtYm(k,3）)； S_res(k，4)=bus(J,2)2(conj(YtYm（k,3）+conj(YtYm(k,5)）bus(I，2）bus(J，2）*（conj（cos（bus（I，3)）+j*sin（bus（I，3)）（cos(bus(J，3）)+j*sin（bus（J，3）)*conj（YtYm（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)； else S_res(k,4)=bus(J，2)2conj（YtYm(k，5）); S_res(k,5）=S_res(k,3)+S_res（k,4）； % 利用公式计算接地支路的潮流 end endendend