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

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; 3 1 0; 3 1 0 ; % 线路数据Branch=%发送接收电阻电感(电导电容)并联 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; % 发电机数据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：% 求解网络的导纳矩阵； 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：%求节点的功率净值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：%求网络节点的注入功率 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：%检查发电机PV节点是否越限； 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：% 计算节点失配功率； 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：% 更新状态变量：电压、相角； 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;% 网络参数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); % 通用参数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( *传统潮流计算结果（快速解耦算法）* )disp( 迭代次数: ,num2str(it)disp( 系统损失的总功率：S=,num2str(sum(PQloss)disp(*)disp( 节点类型 节点电压 电压相角 接收的有功 接收的无功)disp(juzhen1)disp(%节点接收的功率理解：)disp( 当接收功率为“+”时，表示节点从外部（如发电机）吸收功率；)disp( 当接收功率为“-”时，表示节点向外部（如 负荷）输送功率；)disp(*)disp( 线路

20、发送端序号 线路接收节端序号 发送端发出的功率 接收端接收的功率 线路损失的功率)disp(Branch(:,1),Branch(:,2),juzhen2,PQloss)disp(%线路两端接收的功率理解：)disp( 当功率为“+”时，表示向线路发送功率；)disp( 当功率为“-”时，表示从线路吸收功率；)disp( *完成计算* )三、程序运行结果：输出的结果：*传统潮流计算结果（快速解耦算法）* 迭代次数: 27 系统损失的总功率：* 节点类型 节点电压 电压相角 接收的有功 接收的无功 1.0000 1.0600 0 1.7866 1.4799 2.0000 1.0000 -2.06

21、12 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% 节点接收的功率理解： 当接收功率为“+”时，表示节点从外部（如发电机）吸收功率； 当接收功率为“-”时，表示节点向外部（如 负荷）输送功率；* 线路发送端序号 线路接收节端序号 发送端发出的功率 接收端接收的功率 线路损失的功率 1.0000 + 0.0000i 2.0000 + 0.0000i 0.8933 + 0.7400i -0

22、.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.0255i 2.0000 +

23、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% 线路两端接收的功率理解： 当功率为“+”时，表示向线路发送功率； 当功率为“-”时，表示从线路吸收功

24、率；*完成计算*四、程序运行结果：程序一：快速解耦算法第一步：给定初值；包括发电机参数、母线节点数据、初始电压、线路参数等；第二步：计算导纳矩阵 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、Q。1、先统一处理所有节点；P=PNET-PCAL;Q=QNET-QCAL;2、对于发电机平衡节点（slack）的处理；平衡节点（slack）的失配功率P、Q，不是由上述的计算方法得出的。而是当所有节点功率确定以后，按照功率守恒，最后计算得出的。因此，平衡节点（slack）的失配功率不应该包含在上述方程里面。所以如下