电力系统潮流计算完整程序及详细理论说明

1、电力系统潮流计算完整程序及详细理论说明秦羽风在我刚开始学习潮流程序时，总是找不到一个正确的程序开始模仿学习。后来经过多方努力，终于自己写出了一个结构清晰、完整的潮流程序。此程序是一个通用的程序，只需要修改输入数据的子函数(PowerFlowsData_K里面的母线、支路、发电机、负荷，就能算任意一个网络结构的交流系统潮流。很适合初学者学习。为了帮助电力系统的同学一起学习，我将我编写的潮流计算程序分享下来给大家；此程序是在基于牛顿拉夫逊算法的基础上，编写的快速解耦算法。每一个子程序我都有备注说明。如果有不对的地方，希望大家指正！下文中呈现的顺序为：网络结构、子程序、主程序、运算结果、 程序设计理

2、论说明。一、网络结构：5节点网络如下图。二、子程序（共有9个子程序）子程序1 ：（其他系统，只需要修改Bus、Branch> Generator> Load,这四个矩阵就行了）function Bus,Branch,Generator,Load=PowerFlowsData_K%节点数据%类型：1-平衡节点；2-发电机PV节点；3-负荷PQ节点；4-发电机PQ节点;Bus=% 类型电压相角1 1.06 0;2 1 0;3 1 0;4 1 0;5 1 0;%线路数据Branch=% 发送接收电阻电感（电导 电容）并联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;%发电机数据Generator=% 节点定有功定无功(上限 下限)无功1 0 0 5 -5;2 0,4 0 3 -3；%负载数据Load=% 节点定有功 定无功2 0.2 0.1;3 0.45 0.15;4 0.4 0.05;5 0.6 0.1;子程序2:%求解网络的导纳矩阵；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:%求节点的功率净值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:%求网络节点的注入功率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:%佥查发电机pv节点是否越限；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:%计算节点失配功率；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:%更新状态变量：电压、相角；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;%网络参数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);%通用参数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传统潮流计算结果（快速解）disp( disp( disp('迭代次数：'，num2str(it)系统损

19、失的总功率： S=' ,num2str(sum(PQloss)I*disp(节点类型节点电压电压相角）接收的有功接收的无功）disp(juzhen1)disp(disp(disp(disp(%节点接收的功率理解：当接收功率为“当接收功率为“'）+ ”时，表示节点从外部（如发电机）吸收功率;”时，表示节点向外部（如负荷）输送功率;*disp（' 线路发送端序号功率 线路损失的功率）线路接收节端序号发送端发出的功率')')接收端接收的disp(Branch(:,1),Branch(:,2),juzhen2,PQloss')disp( disp( di

20、sp( disp(%线路两端接收的功率理解：）当功率为“ + 当功率为“-时，表示向线路发送功率;时，表示从线路吸收功率;* I')')完成计算）三、程序运行结果：输出的结果：*传统潮流计算结果（快速解耦算法）*迭代次数：27系统损失的总功率：S=0.061222-0.10777i*节点类型节点电压电压相角接收的有功接收的无功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%节点接收的功率理解：当接收功率为“ +”时，表示节点从外部当接收功率为“-”时，表示节点向外部*线路接收节端序号线路发送端序号 线路损失的功率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（如发电机）吸收功率；（如 负荷）输送功率;发送端发出的功率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%线路两端接收的功率理解：当功率为“ +”时，表示向线路发送功率;当功率为“-”时，表示从线路吸收功率；* 完成计算 *四、程序运行结果：程序一：快速解耦算法第一步：给定初值；包括发电机参数、母线节点数据、初始电压、线路参数等；第二步：计算导纳矩阵YR,YI;第三步：计算节点功率净值；PNET=PG-PL;QNET=QG-QL;第四步：计算节点注入功率；PCAL QCAL第五步：检查发电机 PV节点的无功功率 QG是否越限；根据发电机PV节点：QCAL是否在QGmax、QGmin内，判断是否越线。不越限，不处理；如果越限，如下处理：1、将发电机节点 QG固定在 QGmax、QGmin上。2、改变发电机PV节点类型，改为 PQ节点；3、更新次发电机节点的 无功功率净值：QNET=QG-QL ;第六步：计算所有节点的失配功率： P> AQo1、先统一处理所有节点； P=PNET-PCAL; Q=QNET-QCAL;2、对于发电机平衡节点（slack）的处理；平衡节点（slack