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

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

2、論說明。一、網(wǎng)絡(luò)結(jié)構(gòu)：5節(jié)點網(wǎng)絡(luò)如下圖。二、子程序（共有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;4 1 0;5 1 0;%線路數(shù)據(jù)Branch=% 發(fā)送接收電阻電感（電導(dǎo) 電容）并聯(lián)120.020.0600.06;130.080.2400.05;23

3、0.060.1800.04;240.060.1800.04;250.040.1200.03;340.010.0300.02;450.080.2400.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)絡(luò)的導(dǎo)納矩陣；functionYR,YI=YBus_K(nbb,ntl,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep)YR=zeros(n

4、bb,nbb);YI=zeros(nbb,nbb);for kk=1:ntlii=tlsend(kk);jj=tlrec(kk);demon=tlresis(kk)A2+tlreac(kk)A2;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;YI(ii,jj)=YI(ii,jj)+tlreac(kk)/demon;YR(jj,ii)=YR(jj,ii)-

5、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é)點的功率凈值functionPNET,QNET=NetPowers_K(nbb,ngn,nld,genbus,PGEN,QGEN,PLOAD,QLOAD,10ad bus)PNET=zeros(1,nbb);QNET=zeros(1,nbb);for

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

7、PCAL(ii)+VM(ii)*VM(jj)*(YR(ii,jj)*cos(VA(ii)-VA(jj)+YI(ii,j j)*sin(VA(ii)-VA(jj);QCAL(ii)=QCAL(ii)+VM(ii)*VM(jj)*(YR(ii,jj)*sin(VA(ii)-VA(jj)-YI(ii,j j)*cos(VA(ii)-VA(jj);end end子程序5:%僉查發(fā)電機pv節(jié)點是否越限；FunctionQNET,bustype=GeneratorLimit_K(it,ngn,genbus,bustype,QCAL,QMAX,QMIN, nld,loadbus,QLOAD,QNET)fla

8、g1=0;if it>2for kk=1:ngnii=genbus(kk);if bustype(ii)=2if 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;endif flag1=1for jj=1:nldif loadbus(jj)=ii QNET(ii)=QNET(ii)-QLOAD(jj);endendendendendend子程序6:%計算節(jié)點失配功率；FunctionD

9、P,DQ,flag=MismatchPowers_K(nbb,PNET,QNET,PCAL,QCAL,bustype,tol,fla g)DPQ=zeros(1,2*nbb);DP=zeros(1,nbb);DQ=zeros(1,nbb);kk=1;for ii=1:nbbDP(ii)=PNET(ii)-PCAL(ii);DQ(ii)=QNET(ii)-QCAL(ii);if bustype(ii)=1DP(ii)=0;DQ(ii)=0;elseif bustype(ii)=2DQ(ii)=0;endDPQ(kk)=DP(ii);DPQ(kk+1)=DQ(ii);kk=kk+2;endif a

10、bs(DPQ)<tolflag=1;end子程序7:%求快速解耦算法中的B1,B2function B1,B2=FastDecoupled_K(YI,nbb,bustype)B1=-YI;for ii=1:nbbif bustype(ii)=1;for jj=1:nbbif ii=jjB1(ii,ii)=1;elseB1(ii,jj)=0;B1(jj,ii)=0;endend end endB2=B1;for ii=1:nbbif bustype(ii)=2; for jj=1:nbb if ii=jj B2(ii,ii)=1;else B2(ii,jj)=0; B2(jj,ii)=0;

11、endendendend子程序8:%更新狀態(tài)變量：電壓、相角；function VM,VA=StateVariablesUpdates_K(nbb,VA,VM,DVA,DVM)for ii=1:nbbVA(ii)=VA(ii)+DVA(ii);VM(ii)=VM(ii)+DVM(ii)*VM(ii);end子程序9:%計算支路潮流；functionPQsend,PQrec,PQloss,PQbus=PQFlows_K(ntl,VM,VA,tlsend,tlrec,tlresis, tlreac,tlcond,tlsuscep,nbb,nld,ngn,genbus,loadbus,PLOAD,Q

12、LOAD) PQsend=zeros(1,ntl);PQrec=zeros(1,ntl);PQloss=zeros(1,ntl);for ii=1:ntlVsend=VM(tlsend(ii)*(cos(VA(tlsend(ii)+1i*sin(VA(tlsend(ii);Vrec=VM(tlrec(ii)*(cos(VA(tlrec(ii)+1i*sin(VA(tlrec(ii);serice_RL=tlresis(ii)+1i*tlreac(ii);shunt_RC=tlcond(ii)+1i*tlsuscep(ii);current=(Vsend-Vrec)/serice_RL+Vsen

13、d*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);endPQbus=zeros(1,nbb);for ii=1:ntlPQbus(tlsend(ii)=PQbus(tlsend(ii)+PQsend(tlsend(ii);PQbus(tlrec(ii)=PQbus(tlrec(ii)+PQrec(tlrec(ii);endfor ii=

14、1:nldjj=loadbus(ii);for kk=1:ngnll=genbus(kk);if ll=jjPQbus(genbus(kk)=PQbus(genbus(kk)+PLOAD(ii)+1i*QLOAD(ii);endendend二、主程序clear allclcBus,Branch,Generator,Load=PowerFlowsData_K;%網(wǎng)絡(luò)參數(shù)nbb=length(Bus(:,1);bustype=Bus(:,1);VM=Bus(:,2);VA=Bus(:,3);ntl=length(Branch(:,1);tlsend=Branch(:,1);tlrec=Branch

15、(:,2);tlresis=Branch(:,3);tlreac=Branch(:,4);tlcond=Branch(:,5);tlsuscep=Branch(:,6);ngn=size(Generator(:,1);genbus=Generator(:,1);PGEN=Generator(:,2);QGEN=Generator(:,3);QMAX=Generator(:,4);QMIN=Generator(:,5);nld=length(Load(:,1);loadbus=Load(:,1);PLOAD=Load(:,2);QLOAD=Load(:,3);%通用參數(shù)itmax=100;nma

16、x=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,QLOAD,10ad bus);while ( it<itmax && flag=0)PCAL,QCAL=NetPowers_in_K(nbb,VM,VA,YR,YI);QNET,bustype=GeneratorLimit_K(it,ngn,genbus,bustype

17、,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=B2DQ'VM,VA=StateVariablesUpdates_K(nbb,VA,VM,DVA,DVM);it=it+1;endPQsend,PQrec,PQloss,PQbus=PQFlows_K(ntl,VM,VA,tlsend,tlrec,tlresi

18、s, 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(:,5)=imag(PQbus);juzhen2(:,1)=PQsend;juzhen2(:,2)=PQrec;disp（耦算法）* I傳統(tǒng)潮流計算結(jié)果（快速解）disp( disp( disp('迭代次數(shù)：'，num2str(it)系統(tǒng)損

19、失的總功率： S=' ,num2str(sum(PQloss)I*disp(節(jié)點類型節(jié)點電壓電壓相角）接收的有功接收的無功）disp(juzhen1)disp(disp(disp(disp(%節(jié)點接收的功率理解：當接收功率為“當接收功率為“'）+ ”時，表示節(jié)點從外部（如發(fā)電機）吸收功率;”時，表示節(jié)點向外部（如負荷）輸送功率;*disp（' 線路發(fā)送端序號功率 線路損失的功率）線路接收節(jié)端序號發(fā)送端發(fā)出的功率')')接收端接收的disp(Branch(:,1),Branch(:,2),juzhen2,PQloss')disp( disp( di

20、sp( disp(%線路兩端接收的功率理解：）當功率為“ + 當功率為“-時，表示向線路發(fā)送功率;時，表示從線路吸收功率;* I')')完成計算）三、程序運行結(jié)果：輸出的結(jié)果：*傳統(tǒng)潮流計算結(jié)果（快速解耦算法）*迭代次數(shù)：27系統(tǒng)損失的總功率：S=0.061222-0.10777i*節(jié)點類型節(jié)點電壓電壓相角接收的有功接收的無功1.00001.060001.78661.47992.00001.0000-2.06121.05100.42953.00000.9872-4.6367-0.2375-0.03223.00000.9841-4.9570-0.2679-0.03393.0000

21、0.9717-5.7649-1.0689-0.0966%節(jié)點接收的功率理解：當接收功率為“ +”時，表示節(jié)點從外部當接收功率為“-”時，表示節(jié)點向外部*線路接收節(jié)端序號線路發(fā)送端序號 線路損失的功率1.0000 + 0.0000i1.0000 + 0.0000i2.0000 + 0.0000i2.0000 + 0.0000i2.0000 + 0.0000i3.0000 + 0.0000i4.0000 + 0.0000i2.0000 + 0.0000i3.0000 + 0.0000i3.0000 + 0.0000i4.0000 + 0.0000i5.0000 + 0.0000i4.0000 +

22、0.0000i5.0000 + 0.0000i（如發(fā)電機）吸收功率；（如 負荷）輸送功率;發(fā)送端發(fā)出的功率0.8933 + 0.7400i0.4179 + 0.1682i0.2447 - 0.0252i0.2771 - 0.0172i0.5466 + 0.0556i0.1939 + 0.0286i0.0660 + 0.0052i-0.8685 - 0.7291i-0.4027 - 0.1751i-0.2411 - 0.0035i-0.2725 - 0.0083i-0.5344 - 0.0483i-0.1935 - 0.0469i-0.0656 - 0.0517i接收端接收的功率0.0249 -

23、 0.0109i0.0152 + 0.0069i0.0036 + 0.0287i0.0046 + 0.0255i0.0122 - 0.0073i0.0004 + 0.0182i0.0004 + 0.0465i%線路兩端接收的功率理解：當功率為“ +”時，表示向線路發(fā)送功率;當功率為“-”時，表示從線路吸收功率；* 完成計算 *四、程序運行結(jié)果：程序一：快速解耦算法第一步：給定初值；包括發(fā)電機參數(shù)、母線節(jié)點數(shù)據(jù)、初始電壓、線路參數(shù)等；第二步：計算導(dǎo)納矩陣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> AQo1、先統(tǒng)一處理所有節(jié)點； P=PNET-PCAL; Q=QNET-QCAL;2、對于發(fā)電機平衡節(jié)點（slack）的處理；平衡節(jié)點（slack