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

及详细理论说明

——秦羽风

在我刚开始学习潮流程序时，总是找不到一个正确的程序开始模仿学习。后来经过多方努力，终于自己写出了一个结构清晰、完整的潮流程序。此程序是一个通用的程序，只需要修改输入数据的子函数（PowerFlowsData_K）里面的母线、支路、发电机、负荷，就能算任意一个网络结构的交流系统潮流。很适合初学者学习。

为了帮助电力系统的同学一起学习，我将我编写的潮流计算程序分享下来给大家；此程序是在基于牛顿拉夫逊算法的基础上，编写的快速解耦算法。每一个子程序我都有备注说明。如果有不对的地方，希望大家指正！

下文中呈现的顺序为：网络结构、子程序、主程序、运算结果、程序设计理论说明。

一、网络结构：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.24    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,tlsend,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;

    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);

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,YI)

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)));

    end

end

子程序5：

%检查发电机PV节点是否越限；

 Function [QNET,bustype]=GeneratorLimit_K(it,ngn,genbus,bustype,QCAL,QMAX,QMIN,nld,loadbus,QLOAD,QNET)

flag1=0;

if it>2

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)            QNET(ii)=QMIN(kk);

            bustype(ii)=3;

            flag1=1;

        end

        if  flag1==1

            for jj=1:nld

                if  loadbus(jj)==ii

                    QNET(ii)=QNET(ii)-QLOAD(jj);

                end

            end

        end

    end

end

end

子程序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;

for ii=1:nbb

    DP(ii)=PNET(ii)-PCAL(ii);

    DQ(ii)=QNET(ii)-QCAL(ii);

    if  bustype(ii)==1

        DP(ii)=0;

        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)    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

    end

end

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

    end

end

子程序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,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(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,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

    end

end

二、主程序

clear all

clc

[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( it    [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

[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(:,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('    线路发送端序号     线路接收节端序号    发送端发出的功率    接收端接收的功率    线路损失的功率')

disp([Branch(:,1),Branch(:,2),juzhen2,PQloss'])

disp('%线路两端接收的功率理解：')

disp('      当功率为“+”时，表示向线路发送功率；')

disp('      当功率为“-”时，表示从线路吸收功率；')

disp(' *************************************************完成计算****************************************************** ')

三、程序运行结果：

输出的结果：

***********传统潮流计算结果（快速解耦算法）*********

    迭代次数:   27

    系统损失的总功率：S=0.061222-0.10777i

*****************************************************

   节点类型   节点电压  电压相角  接收的有功  接收的无功

    1.0000    1.0600         0    1.7866    1.4799

    2.0000    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.79   -1.06   -0.0966

% 节点接收的功率理解：

      当接收功率为“+”时，表示节点从外部（如发电机）吸收功率；

      当接收功率为“-”时，表示节点向外部（如  负荷）输送功率；

*****************************************************************

    线路发送端序号     线路接收节端序号    发送端发出的功率    接收端接收的功率    线路损失的功率

   1.0000 + 0.0000i   2.0000 + 0.0000i   0.33 + 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.0255i

   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

% 线路两端接收的功率理解：

      当功率为“+”时，表示向线路发送功率；

      当功率为“-”时，表示从线路吸收功率；

****************完成计算****************

四、程序运行结果：

程序一：快速解耦算法

第一步：给定初值；

 包括发电机参数、母线节点数据、初始电压、线路参数等；

第二步：计算导纳矩阵 [YR,YI] ；

第三步：计算节点功率净值；

第四步：计算节点注入功率；

  、QCAL；

第五步：检查发电机PV节点的无功功率QG是否越限；

  根据发电机PV节点： QCAL 是否在 [QGmax、QGmin]内，判断是否越线。

  不越限，不处理；如果越限，如下处理：

1、将发电机节点QG固定在QGmax、QGmin上。

2、改变发电机PV节点类型，改为PQ节点；

3、更新次发电机节点的 无功功率净值：QNET=QG-QL；

第六步：计算所有节点的失配功率：△P、△Q。

  、先统一处理所有节点；

   △P=PNET-PCAL;

   △Q=QNET-QCAL;

  、对于发电机平衡节点（slack）的处理；

平衡节点（slack）的失配功率△P、△Q，不是由上述的计算方法得出的。而是当所有节点功率确定以后，按照功率守恒，最后计算得出的。

因此，平衡节点（slack）的失配功率不应该包含在上述方程里面。所以如下处理：

对于平衡节点：△P=0、△Q=0；

3、对于发电机PV节点的处理；

 对于发电机PV节点，要使节点电压V恒定，使节点吸收或发出无功功率，且次无功功率不能越限。也就是说，此节点的无功功率Q不是按照上述计算方法得到的，而是有节点电压V确定的。

 因此，上述计算过程不应该包含PV节点的失配无功功率△Q。由此分析，发电机PV节点如下处理：△Q=0；

第七步：判断收敛情况；

  判断失配功率是否满足误差精度：

   △P △Q])  ？  tol；

  如果满足，跳出运行，计算结束。得到运行结果：节点电压幅值VM、相角VA。

  如果不满足，继续运行。

第八步：计算雅克比举证 [B1,B2]；

 在快速解耦算法中：

 在其中，B1、B2就相当于牛顿拉夫逊法中的雅克比矩阵JAC。且有，B1与节点导纳矩阵的虚部接近一致；区别在于B1不包含与平衡节点有关的行和列，如下处理：

1、与平衡节点（slack）相关的行列，需要特殊处理；

1）此节点的 失配有功 对 本节点 相角 的偏导数，取1；

2）此节点的 失配无功 对 本节点 电压 的偏导数，取1；

3）其余与此节点失配有功、无功相关的所有行，都取0；

4）其余对节点相角、电压 求偏导数的项，都取0；

B2除了按照上述平衡节点处理后，还需要对发电机PV节点进行如下处理：

3、与PV节点相关的行列，需要特殊处理；

1）此节点的 失配有功 对 本节点 相角 的偏导数，取1；

2）此节点的 失配无功 对 本节点 电压 的偏导数，取1；

3）其余与此节点失配有功、无功相关的所有行，都取0；

4）其余对节点相角、电压 求偏导数的项，都取0；

第九步：更新状态变量VA、VM。

  理论基础：求方程：

  让方程在  处进行泰勒展开，忽略高阶项：

  第i次迭代后得出的公式为：