(文章复现)电力系统碳排放流的计算方法初探(含matlab代码)

参考文献：周天睿,康重庆,徐乾耀,陈启鑫.电力系统碳排放流的计算方法初探[J].电力系统自动化,2012,36(11):44-49.

        双碳目标最近在电力系统领域还是挺火的，碳排放流的计算也是比较热门的话题，这篇论文可以说是最基础的模型，被引次数也是很高的，我就写个博客简单复现一下这篇论文，朋友们有需要的可以参考。

一、基本原理

        下面以IEEE14节点系统为例，说明一下电力系统碳排放流分析的基本原理和步骤。假设系统共有N个节点，K个节点存在发电机功率注入，M个节点存在负荷。

1.1潮流计算

                碳排放流和电力系统的潮流是密切相关的但又有所不同，整体来说，潮流计算是碳排放流计算的基础。因此在计算碳排放流之前首先要对系统进行潮流计算。潮流计算的方法有很多，像最常用的牛拉法、PQ分解法，或者在matlab里调用matpower就能直接出结果，这里不再过多介绍，具体步骤可以参考牛顿法和P-Q分解法IEEE14系统潮流计算（附matlab代码）

        得到潮流计算的结果后，需要用到支路潮流、发电机注入功率、节点注入功率和负荷大小等数据用于碳排放流的计算。

1.2相关矩阵介绍

（1）支路潮流分布矩阵

        支路潮流分布矩阵为N阶方阵,用表示。定义该矩阵的目的是为了描述电力系统的有功潮流分布，计算公式如下：

式中，表示流经支路ij的有功功率。

%% 支路潮流分布矩阵PB
for k=1:length(branch(:,1))
    branch(k,14)=PB(branch(k,1),branch(k,2));
    branch(k,15)=PB(branch(k,2),branch(k,1));
end
for k=1:14
    for kk=1:14
        if PB(k,kk)<0
           PB(k,kk)=0;
        end
    end
end

（2）机组注入分布矩阵

        机组注入分布矩阵为K X N阶矩阵，用表示。定义该矩阵的目的是为了描述所有发电机组向系统中注入的有功功率，计算公式如下：

式中，表示发电机的集合。这个式子的含义就是若第k台发电机组接入节点j，则等于发电机注入功率，否则为零。

%% 机组注入分布矩阵
PG=zeros(K,N);
for k=1:length(gen(:,1))
    PG(k,gen(k,1))=gen(k,2);
end

（3）负荷分布矩阵

        负荷分布矩阵为MXN阶矩阵，用表示。定义该矩阵的目的是描述所有用电负荷与电力系统的连接关系以及有功负荷量，计算公式如下：

式中，表示负荷集合。这个式子的含义就是若第m个负荷接入节点j，则等于有功负荷大小，否则为零。

PL=zeros(M,N);
bus_copy=bus;
for k=1:M
    for kk=1:N
        if bus_copy(kk,3)>0
            PL(k,kk)=bus_copy(kk,3);
            bus_copy(kk,3)=0;
            break
        end
    end
end

（4）节点通量矩阵

        节点有功通量矩阵为N阶对角阵,用表示，计算公式为：

式中，为引入的辅助变量，是一个元素全为1的行向量。

%% 辅助变量PZ
xigama=ones(1,K+N);
PZ=[PB;PG];
%% 节点有功通量矩阵
PN=diag(xigama*PZ);

1.3碳排放流相关计算

（1）发电机组碳排放强度向量

        发电机组碳排放强度向量是一个K维的列向量，用表示，一般都作为已知条件。

（2）节点碳势向量

        节点碳势向量是一个N维的列向量，用表示，计算公式如下：

%% 节点碳势分布向量
EN=(PN-PB')^(-1)*PG'*EG;

（2）支路碳流率分布矩阵

        计算得到节点碳势向量后，和潮流计算类似，可进一步得到系统各个支路的碳流率。由此定义支路碳流率分布矩阵为N阶方阵，用表示，计算公式为：

%% 支路碳硫率分布矩阵
RB=PB*diag(EN);

（3）负荷碳流率向量

计算得到节点碳势向量后，节点负荷的用电碳排放强度与该节点碳势相等。结合负荷分布矩阵，可得所有负荷对应的碳流率,物理意义为发电侧为供应节点负荷每单位时间产生的碳排放量。则负荷碳流率向量用表示，计算公式为：

%% 负荷碳流率向量
RL=PL*EN;

二、运行结果和原文献对比

        由于原文献中未考虑支路损耗，采用直流潮流方法求系统潮流分布，我在代码里用的是极坐标下牛拉法，所以结果会有些误差，但基本一致。

2.1节点有功通量和节点碳势 

（1）原文结果

（2）计算结果

 2.2负荷碳流率

（1）原文结果

（2）计算结果

 PS：完整代码戳这个链接可以看到：

《电力系统碳排放流的计算方法初探》matlab代码复现