電力系統(tǒng)潮流計算完整程序及詳細理論說明

1、電力系統(tǒng)潮流計算完整程序及詳細理論說明秦羽風在我剛開始學習潮流程序時，總是找不到一個正確的程序開始模仿學習。后來經(jīng)過多方努力，終于自己寫出了一個結(jié)構(gòu)清晰、完整的潮流程序。此程序是一個通用的程序，只需要修改輸入數(shù)據(jù)的子函數(shù)（PowerFlowsData_K）里面的母線、支路、發(fā)電機、負荷，就能算任意一個網(wǎng)絡結(jié)構(gòu)的交流系統(tǒng)潮流。很適合初學者學習。為了幫助電力系統(tǒng)的同學一起學習，我將我編寫的潮流計算程序分享下來給大家；此程序是在基于牛頓拉夫遜算法的基礎(chǔ)上，編寫的快速解耦算法。每一個子程序我都有備注說明。如果有不對的地方，希望大家指正！下文中呈現(xiàn)的順序為：網(wǎng)絡結(jié)構(gòu)、子程序、主程序、運算結(jié)果、程序設(shè)計理

2、論說明。一、網(wǎng)絡結(jié)構(gòu)：5節(jié)點網(wǎng)絡如下圖。二、子程序（共有9個子程序）子程序1：（其他系統(tǒng)，只需要修改Bus、Branch、Generator、Load，這四個矩陣就行了）function Bus,Branch,Generator,Load=PowerFlowsData_K % 節(jié)點數(shù)據(jù)% 類型：1-平衡節(jié)點；2-發(fā)電機PV節(jié)點；3-負荷PQ節(jié)點；4-發(fā)電機PQ節(jié)點；Bus=%類型電壓相角 1 1.06 0; 2 1 0; 3 1 0; 3 1 0; 3 1 0 ; % 線路數(shù)據(jù)Branch=%發(fā)送接收電阻電感(電導電容)并聯(lián) 1 2 0.02 0.06 0 0.06; 1 3 0.08 0.2

3、4 0 0.05; 2 3 0.06 0.18 0 0.04; 2 4 0.06 0.18 0 0.04; 2 5 0.04 0.12 0 0.03; 3 4 0.01 0.03 0 0.02; 4 5 0.08 0.24 0 0.05; % 發(fā)電機數(shù)據(jù)Generator=%節(jié)點定有功定無功(上限下限)無功 1 0 0 5 -5; 2 0.4 0 3 -3; % 負載數(shù)據(jù)Load=%節(jié)點定有功定無功 2 0.2 0.1; 3 0.45 0.15; 4 0.4 0.05; 5 0.6 0.1;子程序2：% 求解網(wǎng)絡的導納矩陣； function YR,YI=YBus_K(nbb,ntl,tlse

4、nd,tlrec,tlresis,tlreac,tlcond,tlsuscep) YR=zeros(nbb,nbb);YI=zeros(nbb,nbb); for kk=1:ntl ii=tlsend(kk); jj=tlrec(kk); demon=tlresis(kk)2+tlreac(kk)2; YR(ii,ii)=YR(ii,ii)+tlresis(kk)/demon+0.5*tlcond(kk); YI(ii,ii)=YI(ii,ii)-tlreac(kk)/demon+0.5*tlsuscep(kk); YR(ii,jj)=YR(ii,jj)-tlresis(kk)/demon;

5、YI(ii,jj)=YI(ii,jj)+tlreac(kk)/demon; YR(jj,ii)=YR(jj,ii)-tlresis(kk)/demon; YI(jj,ii)=YI(jj,ii)+tlreac(kk)/demon; YR(jj,jj)=YR(jj,jj)+tlresis(kk)/demon+0.5*tlcond(kk); YI(jj,jj)=YI(jj,jj)-tlreac(kk)/demon+0.5*tlsuscep(kk); end子程序3：%求節(jié)點的功率凈值function PNET,QNET=NetPowers_K(nbb,ngn,nld,genbus,PGEN,QGEN

6、,PLOAD,QLOAD,loadbus)PNET=zeros(1,nbb);QNET=zeros(1,nbb); for kk=1:ngn ii=genbus(kk); PNET(ii)=PNET(ii)+PGEN(kk); QNET(ii)=QNET(ii)+QGEN(kk); end for kk=1:nld ii=loadbus(kk); PNET(ii)=PNET(ii)-PLOAD(kk); QNET(ii)=QNET(ii)-QLOAD(kk);end子程序4：%求網(wǎng)絡節(jié)點的注入功率 function PCAL,QCAL=NetPowers_in_K(nbb,VM,VA,YR,Y

7、I)PCAL=zeros(1,nbb);QCAL=zeros(1,nbb); for ii=1:nbb for jj=1:nbb PCAL(ii)=PCAL(ii)+VM(ii)*VM(jj)*(YR(ii,jj)*cos(VA(ii)-VA(jj)+YI(ii,jj)*sin(VA(ii)-VA(jj); QCAL(ii)=QCAL(ii)+VM(ii)*VM(jj)*(YR(ii,jj)*sin(VA(ii)-VA(jj)-YI(ii,jj)*cos(VA(ii)-VA(jj); endend子程序5：%檢查發(fā)電機PV節(jié)點是否越限； Function QNET,bustype=Genera

8、torLimit_K(it,ngn,genbus,bustype,QCAL,QMAX,QMIN,nld,loadbus,QLOAD,QNET) flag1=0; if it2 for kk=1:ngn ii=genbus(kk); if bustype(ii)=2 if QCAL(ii)QMAX(kk) QNET(ii)=QMAX(kk); bustype(ii)=3; flag1=1; elseif QCAL(ii)QMIN(kk) QNET(ii)=QMIN(kk); bustype(ii)=3; flag1=1; end if flag1=1 for jj=1:nld if loadbu

9、s(jj)=ii QNET(ii)=QNET(ii)-QLOAD(jj); end end end endendend子程序6：% 計算節(jié)點失配功率； FunctionDP,DQ,flag=MismatchPowers_K(nbb,PNET,QNET,PCAL,QCAL,bustype,tol,flag)DPQ=zeros(1,2*nbb);DP=zeros(1,nbb);DQ=zeros(1,nbb);kk=1; for ii=1:nbb DP(ii)=PNET(ii)-PCAL(ii); DQ(ii)=QNET(ii)-QCAL(ii); if bustype(ii)=1 DP(ii)=0

10、; DQ(ii)=0; elseif bustype(ii)=2 DQ(ii)=0; end DPQ(kk)=DP(ii); DPQ(kk+1)=DQ(ii); kk=kk+2;end if abs(DPQ)tol flag=1;end子程序7：%求快速解耦算法中的 B1,B2 function B1,B2=FastDecoupled_K(YI,nbb,bustype) B1=-YI;for ii=1:nbb if bustype(ii)=1; for jj=1:nbb if ii=jj B1(ii,ii)=1; else B1(ii,jj)=0; B1(jj,ii)=0; end end e

11、ndend B2=B1;for ii=1:nbb if bustype(ii)=2; for jj=1:nbb if ii=jj B2(ii,ii)=1; else B2(ii,jj)=0; B2(jj,ii)=0; end end endend子程序8：% 更新狀態(tài)變量：電壓、相角； function VM,VA=StateVariablesUpdates_K(nbb,VA,VM,DVA,DVM) for ii=1:nbb VA(ii)=VA(ii)+DVA(ii); VM(ii)=VM(ii)+DVM(ii)*VM(ii);end子程序9：% 計算支路潮流； function PQsend

12、,PQrec,PQloss,PQbus=PQFlows_K(ntl,VM,VA,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,nbb,nld,ngn,genbus,loadbus,PLOAD,QLOAD)PQsend=zeros(1,ntl);PQrec=zeros(1,ntl);PQloss=zeros(1,ntl); for ii=1:ntl Vsend=VM(tlsend(ii)*(cos(VA(tlsend(ii)+1i*sin(VA(tlsend(ii); Vrec=VM(tlrec(ii)*(cos(VA(tlrec(ii)+1i*sin(

13、VA(tlrec(ii); serice_RL=tlresis(ii)+1i*tlreac(ii); shunt_RC=tlcond(ii)+1i*tlsuscep(ii); current=(Vsend-Vrec)/serice_RL+Vsend*shunt_RC*0.5; PQsend(ii)=Vsend*conj(current); current=(Vrec-Vsend)/serice_RL+Vrec*shunt_RC*0.5; PQrec(ii)=Vrec*conj(current); PQloss(ii)=PQsend(ii)+PQrec(ii);end PQbus=zeros(1

14、,nbb);for ii=1:ntl PQbus(tlsend(ii)=PQbus(tlsend(ii)+PQsend(tlsend(ii); PQbus(tlrec(ii)=PQbus(tlrec(ii)+PQrec(tlrec(ii);end for ii=1:nld jj=loadbus(ii); for kk=1:ngn ll=genbus(kk); if ll=jj PQbus(genbus(kk)=PQbus(genbus(kk)+PLOAD(ii)+1i*QLOAD(ii); end endend二、主程序clear allclcBus,Branch,Generator,Load

15、=PowerFlowsData_K;% 網(wǎng)絡參數(shù)nbb=length(Bus(:,1);bustype=Bus(:,1);VM=Bus(:,2);VA=Bus(:,3); ntl=length(Branch(:,1);tlsend=Branch(:,1);tlrec=Branch(:,2);tlresis=Branch(:,3);tlreac=Branch(:,4);tlcond=Branch(:,5);tlsuscep=Branch(:,6); ngn=size(Generator(:,1);genbus=Generator(:,1);PGEN=Generator(:,2);QGEN=Gen

16、erator(:,3);QMAX=Generator(:,4);QMIN=Generator(:,5); nld=length(Load(:,1);loadbus=Load(:,1);PLOAD=Load(:,2);QLOAD=Load(:,3); % 通用參數(shù)itmax=100;nmax=2*nbb;tol=1e-12;it=1;flag=0; YR,YI=YBus_K(nbb,ntl,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep);PNET,QNET=NetPowers_K(nbb,ngn,nld,genbus,PGEN,QGEN,PLOAD,QL

17、OAD,loadbus); while( ititmax & flag=0) PCAL,QCAL=NetPowers_in_K(nbb,VM,VA,YR,YI); QNET,bustype=GeneratorLimit_K(it,ngn,genbus,bustype,QCAL,QMAX,QMIN,nld,loadbus,QLOAD,QNET); DP,DQ,flag=MismatchPowers_K(nbb,PNET,QNET,PCAL,QCAL,bustype,tol,flag); B1,B2=FastDecoupled_K(YI,nbb,bustype); DVA=B1DP; DVM=B2

18、DQ; VM,VA=StateVariablesUpdates_K(nbb,VA,VM,DVA,DVM); it=it+1;end PQsend,PQrec,PQloss,PQbus=PQFlows_K(ntl,VM,VA,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,nbb,nld,ngn,genbus,loadbus,PLOAD,QLOAD);VA=(VA*180/pi); juzhen1(:,1)=Bus(:,1);juzhen1(:,2)=VM;juzhen1(:,3)=VA;juzhen1(:,4)=real(PQbus);juzhen1(:

19、,5)=imag(PQbus); juzhen2(:,1)=PQsend;juzhen2(:,2)=PQrec;disp( *傳統(tǒng)潮流計算結(jié)果（快速解耦算法）* )disp( 迭代次數(shù): ,num2str(it)disp( 系統(tǒng)損失的總功率：S=,num2str(sum(PQloss)disp(*)disp( 節(jié)點類型 節(jié)點電壓 電壓相角 接收的有功 接收的無功)disp(juzhen1)disp(%節(jié)點接收的功率理解：)disp( 當接收功率為“+”時，表示節(jié)點從外部（如發(fā)電機）吸收功率；)disp( 當接收功率為“-”時，表示節(jié)點向外部（如 負荷）輸送功率；)disp(*)disp( 線路

20、發(fā)送端序號 線路接收節(jié)端序號 發(fā)送端發(fā)出的功率 接收端接收的功率 線路損失的功率)disp(Branch(:,1),Branch(:,2),juzhen2,PQloss)disp(%線路兩端接收的功率理解：)disp( 當功率為“+”時，表示向線路發(fā)送功率；)disp( 當功率為“-”時，表示從線路吸收功率；)disp( *完成計算* )三、程序運行結(jié)果：輸出的結(jié)果：*傳統(tǒng)潮流計算結(jié)果（快速解耦算法）* 迭代次數(shù): 27 系統(tǒng)損失的總功率：S=0.-0.10777i* 節(jié)點類型 節(jié)點電壓 電壓相角 接收的有功 接收的無功 1.0000 1.0600 0 1.7866 1.4799 2.0000

21、 1.0000 -2.0612 1.0510 0.4295 3.0000 0.9872 -4.6367 -0.2375 -0.0322 3.0000 0.9841 -4.9570 -0.2679 -0.0339 3.0000 0.9717 -5.7649 -1.0689 -0.0966% 節(jié)點接收的功率理解： 當接收功率為“+”時，表示節(jié)點從外部（如發(fā)電機）吸收功率； 當接收功率為“-”時，表示節(jié)點向外部（如 負荷）輸送功率；* 線路發(fā)送端序號 線路接收節(jié)端序號 發(fā)送端發(fā)出的功率 接收端接收的功率 線路損失的功率 1.0000 + 0.0000i 2.0000 + 0.0000i 0.8933

22、 + 0.7400i -0.8685 - 0.7291i 0.0249 - 0.0109i 1.0000 + 0.0000i 3.0000 + 0.0000i 0.4179 + 0.1682i -0.4027 - 0.1751i 0.0152 + 0.0069i 2.0000 + 0.0000i 3.0000 + 0.0000i 0.2447 - 0.0252i -0.2411 - 0.0035i 0.0036 + 0.0287i 2.0000 + 0.0000i 4.0000 + 0.0000i 0.2771 - 0.0172i -0.2725 - 0.0083i 0.0046 + 0.02

23、55i 2.0000 + 0.0000i 5.0000 + 0.0000i 0.5466 + 0.0556i -0.5344 - 0.0483i 0.0122 - 0.0073i 3.0000 + 0.0000i 4.0000 + 0.0000i 0.1939 + 0.0286i -0.1935 - 0.0469i 0.0004 + 0.0182i 4.0000 + 0.0000i 5.0000 + 0.0000i 0.0660 + 0.0052i -0.0656 - 0.0517i 0.0004 + 0.0465i% 線路兩端接收的功率理解： 當功率為“+”時，表示向線路發(fā)送功率； 當功率為

24、“-”時，表示從線路吸收功率；*完成計算*四、程序運行結(jié)果：程序一：快速解耦算法第一步：給定初值；包括發(fā)電機參數(shù)、母線節(jié)點數(shù)據(jù)、初始電壓、線路參數(shù)等；第二步：計算導納矩陣 YR,YI ；第三步：計算節(jié)點功率凈值；PNET=PG-PL;QNET=QG-QL;第四步：計算節(jié)點注入功率；PCAL、QCAL；第五步：檢查發(fā)電機PV節(jié)點的無功功率QG是否越限；根據(jù)發(fā)電機PV節(jié)點： QCAL 是否在 QGmax、QGmin內(nèi)，判斷是否越線。不越限，不處理；如果越限，如下處理：1、 將發(fā)電機節(jié)點QG固定在QGmax、QGmin上。2、 改變發(fā)電機PV節(jié)點類型，改為PQ節(jié)點；3、 更新次發(fā)電機節(jié)點的 無功功率凈值：QNET=QG-QL；第六步：計算所有節(jié)點的失配功率：P、Q。1、先統(tǒng)一處理所有節(jié)點；P=PNET-PCAL;Q=QNET-QCAL;2、對于發(fā)電機平衡節(jié)點（slack）的處理；平衡節(jié)點（slack）的失配功率P、Q，不是由上述的計算方法得出的。而是當所有節(jié)點功率確定以后，按照功率守恒，最后計算得出的。因此，平衡節(jié)點（slack）的失配功率不應該包含在上述方程里面。所以