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

电力系统潮流计算完整程序及详细理论说明——秦羽风在我刚开始学习潮流程序时，总是找不到一个正确的程序开始模仿学习。后来经过多方努力，终于自己写出了一个结构清晰、完整的潮流程序。此程序是一个通用的程序，只需要修改输入数据的子函数（PowerFlowsData_K）里面的母线、支路、发电机、负荷，就能算任意一个网络结构的交流系统潮流。很适合初学者学习。为了帮助电力系统的同学一起学习，我将我编写的潮流计算程序分享下来给大家；此程序是在基于牛顿拉夫逊算法的基础上，编写的快速解耦算法。每一个子程序我都有备注说明。如果有不对的地方，希望大家指正！下文中呈现的顺序为：网络结构、子程序、主程序、运算结果、程序设计理论说明。一、网络结构：5节点网络如下图。二、子程序（共有9个子程序）子程序1：（其他系统，只需要修改Bus、Branch、Generator、Load，这四个矩阵就行了）function[Bus,Branch,Generator,Load]=PowerFlowsData_K%%节点数据%类型：1-平衡节点；2-发电机PV节点；3-负荷PQ节点；4-发电机PQ节点；Bus=[%类型电压相角11.060;210;310;310;310];%%线路数据Branch=[%发送接收电阻电感(电导电容)并联120.020.0600.06;130.080.2400.05;230.060.1800.04;240.060.1800.04;250.040.1200.03;340.010.0300.02;450.080.2400.05];%%发电机数据Generator=[%节点定有功定无功(上限下限)无功1005-5;20.403-3];%%负载数据Load=[%节点定有功定无功20.20.1;30.450.15;40.40.05;50.60.1];子程序2：%求解网络的导纳矩阵；function[YR,YI]=YBus_K(nbb,ntl,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep)YR=zeros(nbb,nbb);YI=zeros(nbb,nbb);forkk=1:ntlii=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;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：%求节点的功率净值function[PNET,QNET]=NetPowers_K(nbb,ngn,nld,genbus,PGEN,QGEN,PLOAD,QLOAD,loadbus)PNET=zeros(1,nbb);QNET=zeros(1,nbb);forkk=1:ngnii=genbus(kk);PNET(ii)=PNET(ii)+PGEN(kk);QNET(ii)=QNET(ii)+QGEN(kk);endforkk=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);forii=1:nbbforjj=1:nbbPCAL(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：%检查发电机PV节点是否越限；Function[QNET,bustype]=GeneratorLimit_K(it,ngn,genbus,bustype,QCAL,QMAX,QMIN,nld,loadbus,QLOAD,QNET)flag1=0;ifit2forkk=1:ngnii=genbus(kk);ifbustype(ii)==2ifQCAL(ii)QMAX(kk)QNET(ii)=QMAX(kk);bustype(ii)=3;flag1=1;elseifQCAL(ii)QMIN(kk)QNET(ii)=QMIN(kk);bustype(ii)=3;flag1=1;endifflag1==1forjj=1:nldifloadbus(jj)==iiQNET(ii)=QNET(ii)-QLOAD(jj);endendendendendend子程序6：%计算节点失配功率；Function[DP,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;forii=1:nbbDP(ii)=PNET(ii)-PCAL(ii);DQ(ii)=QNET(ii)-QCAL(ii);ifbustype(ii)==1DP(ii)=0;DQ(ii)=0;elseifbustype(ii)==2DQ(ii)=0;endDPQ(kk)=DP(ii);DPQ(kk+1)=DQ(ii);kk=kk+2;endifabs(DPQ)tolflag=1;end子程序7：%求快速解耦算法中的[B1,B2]function[B1,B2]=FastDecoupled_K(YI,nbb,bustype)B1=-YI;forii=1:nbbifbustype(ii)==1;forjj=1:nbbifii==jjB1(ii,ii)=1;elseB1(ii,jj)=0;B1(jj,ii)=0;endendendendB2=B1;forii=1:nbbifbustype(ii)==2;forjj=1:nbbifii==jjB2(ii,ii)=1;elseB2(ii,jj)=0;B2(jj,ii)=0;endendendend子程序8：%更新状态变量：电压、相角；function[VM,VA]=StateVariablesUpdates_K(nbb,VA,VM,DVA,DVM)forii=1:nbbVA(ii)=VA(ii)+DVA(ii);VM(ii)=VM(ii)+DVM(ii)*VM(ii);end子程序9：%计算支路潮流；function[PQsend,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);forii=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+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);endPQbus=zeros(1,nbb);forii=1:ntlPQbus(tlsend(ii))=PQbus(tlsend(ii))+PQsend(tlsend(ii));PQbus(tlrec(ii))=PQbus(tlrec(ii))+PQrec(tlrec(ii));endforii=1:nldjj=loadbus(ii);forkk=1:ngnll=genbus(kk);ifll==jjPQbus(genbus(kk))=PQbus(genbus(kk))+PLOAD(ii)+1i*QLOAD(ii);endendend二、主程序clearallclc[Bus,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(:,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;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,QLOAD,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=B1\DP';DVM=B2\DQ';[VM,VA]=StateVariablesUpdates_K(nbb,VA,VM,DVA,DVM);it=it+1;end