电网的碳足迹到底该怎么算这是我这几个月一直在折腾的问题。起因是课题组要复现一篇电气工程领域顶级EI期刊上关于电力系统碳排放流分析的文章而且必须在IEEE 14节点系统上跑通。一开始我以为只是把潮流计算跑一遍再把发电机碳排放因子乘上去就完事真正动手才发现碳排放流的计算远不是“功率乘因子”这么简单里面涉及碳流密度的传播、比例分摊原则、平衡节点的处理等一系列容易翻车的细节。这篇就把我完整的复现思路、Matlab代码逻辑以及踩坑记录整理出来希望能帮到正在做双碳方向仿真或者准备复现同类论文的朋友。1. 为什么电力系统需要一本“碳账本”从发电侧核算说起1.1 传统碳核算的盲区电力系统的碳排放传统核算方法是盯住发电厂烧了多少煤、多少天然气乘以对应的排放因子就是发电侧的碳排放总量。这个方法在就地发电、就地消纳的年代问题不大但现在的电网是跨区互联的你在A城市用的电很可能来自几百公里外的B电厂。如果只统计发电侧就会出现一个尴尬局面电厂所在地承担了全部碳排放责任而真正用电的用户却不需要为这些排放负责。这在碳交易、碳考核面前是说不通的。欧盟碳边境调节机制、国内碳市场扩容都在倒逼电网企业把碳排放责任合理分摊到每个负荷节点甚至每个用户头上。可问题是电在网里怎么流动碳就该怎么流动但碳不像有功功率那样有仪表可以直接测。碳排放流分析方法就是在这个背景下被提出来的——它通过潮流计算结果把发电厂的碳排放在电网拓扑中逐级追踪、分摊最后得到每个节点、每条支路、每个负荷的碳流指标。1.2 碳排放流本质上在算什么碳排放流的基本思想很朴素既然电能从发电机经过线路送到负荷那么发电产生的碳排放也应当附着在电能上沿着有功潮流的路径在网络中传播。核心概念有三个节点碳势也叫碳流密度单位是kgCO2/MWh表示该节点上每消纳1MWh电能所对应的碳排放。支路碳流率单位是kgCO2/h表示单位时间内通过某条支路运输的碳排放量。负荷碳流率表示某个负荷节点单位时间内从系统中拿走的碳排放量。这三个量建立了发电排放—网络传输—负荷消耗的完整碳账本让每一个节点的用电都有了可追踪的碳排放强度。搞懂这个概念你再看任何一篇碳排放流相关的EI论文核心模型基本都是这套体系。1.3 为什么选IEEE 14节点作为验证系统IEEE 14节点测试系统是电力系统分析里最经典的算例之一14个节点、20条支路、5台发电机规模不大但麻雀虽小五脏俱全有环网、有变压器支路、有纯负荷节点、也有纯发电机节点节点8只有发电机没有负荷。这种结构特别适合用来验证碳流算法的正确性——结果可以用手算局部验证画图又足够直观。而且绝大多数EI论文和学位论文里的碳排放流算例要么用14节点要么用30节点先跑通14节点后面扩展到30节点、118节点都是顺水推舟的事。2. 碳排放流的核心模型节点碳势、比例分摊与线性方程组推导2.1 基本假设碳跟着有功走碳排放流计算有一个核心假设碳排放流完全跟随有功功率的流动方向与无功功率无关。这条假设的物理依据是电网中的碳排放本质上是发电环节化石燃料燃烧的产物它通过电能传输被携带到负荷侧而有功功率才是真正做功的那部分能量。无功功率主要维持电压水平不承担能量传输任务所以在碳流分析里默认不携带碳。另一条关键假设是比例分摊原则。什么意思假如节点i既有本地发电机注入功率又有从上游线路流入的功率那么从节点i流出的所有支路功率以及节点i上的负荷功率都按照各自功率大小等比例分摊节点i的总碳注入。这个原则保证了碳流计算是线性的、可分解的也是后续矩阵求解的基础。2.2 节点碳势的数学定义节点i的碳势定义为$$e_i \frac{\sum_{g \in G_i} P_g \cdot e_g \sum_{j \in U_i} P_{ji} \cdot e_j}{\sum_{g \in G_i} P_g \sum_{j \in U_i} P_{ji}}$$其中$G_i$是接在节点i上的发电机集合$P_g$是发电机g的有功出力MW$e_g$是发电机g的碳排放强度kgCO2/MWh$U_i$是向节点i注入功率的上游节点集合$P_{ji}$是从节点j流向节点i的支路有功功率MW$e_j$是上游节点j的碳势注意分母是所有注入节点i的有功功率之和不包括负荷。负荷是从节点取走碳不参与碳势的稀释。这一点我在最初实现时搞反过后面会单独说。把上式两边同乘分母并移项得到节点i的碳势方程$$e_i \cdot \left( \sum_{g \in G_i} P_g \sum_{j \in U_i} P_{ji} \right) - \sum_{j \in U_i} P_{ji} \cdot e_j \sum_{g \in G_i} P_g \cdot e_g$$对全网的N个节点都写这个方程就得到一个含N个未知数节点碳势、N个方程的线性方程组。写成矩阵形式$$(D G - A) \cdot \mathbf{e} \mathbf{b}$$其中$D$是对角阵$D_{ii}$等于流入节点i的所有支路功率之和$G$是对角阵$G_{ii}$等于节点i上所有发电机有功出力之和$A$是支路注入矩阵$A_{ij}$表示从节点j流入节点i的支路功率$i \neq j$$\mathbf{b}$是发电碳注入向量$b_i \sum_{g \in G_i} P_g \cdot e_g$这个方程组解出来的$\mathbf{e}$就是每个节点的碳势也是后面所有碳流指标计算的源头。2.3 支路碳流率与负荷碳流率的计算得到节点碳势后一切都顺了。支路碳流率如果支路k的有功潮流从节点i流向节点j且首端有功为$P_{ij}$那么该支路的碳流率为$$R_{ij} P_{ij} \cdot e_i$$这里的物理意义很直接功率从高碳势节点流向低碳势节点时会把碳势带过去反之功率从低碳势节点流向高碳势节点该支路碳流率的方向是相反的。支路两端碳流率之差就是这条支路线损所对应的碳排放这在全网碳平衡校验时会用到。负荷碳流率$$R_i^L P_i^L \cdot e_i$$$P_i^L$是节点i的有功负荷$e_i$是节点i的碳势。这一项最终落到每个用户头上是碳排放责任分摊的核心数据。2.4 完整计算流程整个碳排放流计算可以总结为五步对IEEE 14节点系统做稳态潮流计算得到各节点电压、相角、支路有功潮流、发电机实际出力。根据支路有功潮流的正负判断每条支路的实际功率方向确定上下游关系。组装矩阵$D$、$G$、$A$和向量$\mathbf{b}$。求解线性方程组得到节点碳势。由节点碳势计算支路碳流率、负荷碳流率并做全网碳平衡校验。3. IEEE 14节点系统数据准备Matpower结构与发电机碳排放参数设定3.1 标准数据从哪来复现碳排放流第一步是拿到IEEE 14节点系统的标准潮流数据。我强烈建议直接用Matpower自带的case14.m不要自己手输参数。Matpower是电力系统领域最常用的Matlab开源工具箱内置了标准的IEEE 14节点算例数据包含14个节点bus、20条支路branch、5台发电机gen。用Matpower最省心的地方在于它的数据结构非常规整做潮流计算只要一行runpf(case14)而且计算结果里包含我们需要的所有字段。你可以在Matlab命令行输入open(case14.m)查看具体数据也可以直接运行case14查看系统的节点数、支路数、基准容量等基本参数。注意Matpower中所有功率的单位是MW/Mvar基准容量是100MVA但bus、branch、gen矩阵里的数据已经是标幺值换算后的有名值在case14中通常直接采用100MVA基准这个细节在后续单位换算时要留意。3.2 发电机碳排放强度的设定方案IEEE 14节点系统里有5台发电机分别挂在节点1、2、3、6、8上。但在标准测试数据中节点3、6、8的发电机在基础工况下有功出力为0相当于同步调相机只发无功。如果你直接用原始数据做碳排放流计算会发现全系统只有节点1和节点2两台发电机组在供电碳流分布非常单调压根体现不出多源碳流分摊的复杂性。所以绝大多数EI论文在碳排放流算例中都会对发电机数据进行改造给节点3、6、8的机组分配有功出力并设定不同的燃料类型和碳排放强度。常见的改造方案如下节点发电类型有功出力MW碳排放强度kgCO2/MWh1燃煤机组232.4平衡机由潮流计算确定8502燃煤机组40.08503燃气机组40.04506燃气机组30.04508水电机组30.00这个改造方案不是唯一的。有的文献会让节点8的机组也出力有的文献会给节点1设置更高的排放因子比如900取值差异主要看你要复现的目标论文怎么设定。我的建议是如果目标论文给了明确的排放因子就严格按论文取值如果论文只写了煤电、气电、水电的分类那就用我上面这一组推荐值煤电850、气电450、水电0并在论文里注明假设。需要特别提醒的是节点1是平衡节点它的有功出力不是预先设定的而是潮流计算收敛后自动平衡全系统功率差额得到的结果。所以在碳流计算里节点1的发电碳注入必须用潮流计算后的实际出力而不是你在gen矩阵里初始填的估计值。这个坑我在后面还会展开说。3.3 从Matpower结构到碳流计算的数据映射Matpower的bus、branch、gen矩阵每个字段都有固定的列编号。做碳流计算前我们要手动提取以下几个关键字段bus矩阵的第3列PD节点有功负荷单位MW。branch矩阵的第1列F_BUS、第2列T_BUS支路首端节点、末端节点编号。branch矩阵的第14列PF支路首端有功功率单位MW正值表示功率从首端流向末端。branch矩阵的第15列PT支路末端有功功率正值表示功率从末端流出。gen矩阵的第1列GEN_BUS发电机所在节点编号。gen矩阵的第2列PG发电机有功出力单位MW。在碳流计算中最核心的是PF列首端有功。通过判断PF的正负就能确定支路潮流的实际方向PF为正功率从节点F流向节点TPF为负功率从节点T流向节点F。这里千万不能用PT列直接判断方向我最初就是在这里栽了跟头下面会有专门一节讲。4. Matlab代码实现矩阵装配、求解与碳流指标计算4.1 代码总体框架整个Matlab实现分为四个模块潮流计算、数据提取、碳势求解、碳流指标输出。下面是我在复现过程中最终使用的代码总体结构%% 碳排放流计算主程序 - IEEE 14节点系统 clear; clc; close all; % 1. 载入IEEE 14节点系统并运行潮流计算 mpc loadcase(case14); mpopt mpoption(verbose, 0, out.all, 0); result runpf(mpc, mpopt); % 2. 提取潮流数据 bus result.bus; branch result.branch; gen result.gen; % 3. 设定发电机碳排放强度kgCO2/MWh % 节点1、2为煤电节点3、6为气电节点8为水电 genEmission [850; 850; 450; 450; 0]; % 4. 计算节点碳势 e solveCarbonPotential(bus, branch, gen, genEmission); % 5. 计算支路碳流率和负荷碳流率 branchCarbon calcBranchCarbonFlow(branch, e); loadCarbon bus(:, 3) .* e; % 6. 输出结果 disp(节点碳势(kgCO2/MWh):); disp(e);4.2 核心函数节点碳势求解节点碳势求解是整个碳流计算的心脏。我把它封装成一个独立函数方便后续扩展到其他节点系统时复用function e solveCarbonPotential(bus, branch, gen, genEmission) nb size(bus, 1); % 节点数 ng size(gen, 1); % 发电机数 % 提取支路首端有功和首末端节点编号 PF branch(:, 14); F branch(:, 1); T branch(:, 2); % 初始化矩阵 D zeros(nb, nb); % 支路注入对角阵 A zeros(nb, nb); % 上游注入矩阵 G zeros(nb, nb); % 发电机注入对角阵 b zeros(nb, 1); % 发电碳注入向量 % 遍历支路根据潮流方向组装D和A for k 1:length(PF) if PF(k) 0 % 功率从F流向T D(T(k), T(k)) D(T(k), T(k)) PF(k); A(T(k), F(k)) A(T(k), F(k)) PF(k); else % 功率从T流向F p -PF(k); D(F(k), F(k)) D(F(k), F(k)) p; A(F(k), T(k)) A(F(k), T(k)) p; end end % 遍历发电机组装G和b for g 1:ng i gen(g, 1); % 发电机所在节点 pg gen(g, 2); % 有功出力 if pg 0 G(i, i) G(i, i) pg; b(i) b(i) pg * genEmission(g); end % pg为0或负值的机组不注入碳也不参与碳势计算 end % 组装系统矩阵并求解 M D G - A; e M \ b; end这段代码的核心逻辑就三步第一遍历所有支路根据PF正负判断方向把支路功率灌进对应节点的对角元和非对角元。注意方向处理的技巧我统一用PF判断PF为负时取反这样不管潮流方向如何变化代码都不需要分叉处理。第二遍历所有发电机把有功出力乘碳排放强度作为碳注入向量。这里有个细节如果某台发电机有功出力为0比如原版case14中节点3、6、8的机组它不产生碳注入如果出现负出力比如抽水蓄能电站充电工况也不应该计入碳注入所以用pg 0做判断是稳妥的。第三解线性方程组。这里直接用Matlab的左除运算符\求解。对于14节点这种小系统直接LU分解完全够用不需要迭代法。实际测试中这个方程组的求解耗时在毫秒级几乎可以忽略。4.3 支路碳流率计算节点碳势算出来之后支路碳流率就很简单了function branchCarbon calcBranchCarbonFlow(branch, e) nbr size(branch, 1); branchCarbon zeros(nbr, 1); PF branch(:, 14); F branch(:, 1); T branch(:, 2); for k 1:nbr if PF(k) 0 % 功率从F流向T碳流率 F节点碳势 * 有功功率 branchCarbon(k) PF(k) * e(F(k)); else % 功率从T流向F碳流率 T节点碳势 * 反向功率 branchCarbon(k) (-PF(k)) * e(T(k)); end end end这里要注意支路碳流率的正负号含义是方向不是大小。实际做图或分析时建议单独用一个方向标志数组记录每条支路的碳流方向不要只靠正负号判断。4.4 结果验证的三种方法代码跑通只是开始验证结果正确性才是完美复现的关键。我总结了三层验证方法每一层都能筛出不同的bug第一层全网碳平衡校验。全系统发电碳注入总量应该等于所有负荷碳流率之和加上全网线损碳流率。用代码表示为totalGenCarbon sum(sum(gen(:, 2) .* genEmission)); totalLoadCarbon sum(loadCarbon); % 线损碳流 全网支路首端碳流率之和 - 支路末端碳流率之和 branchLossCarbon 0; for k 1:size(branch, 1) pt branch(k, 15); % 末端有功 if branch(k, 14) 0 branchLossCarbon branchLossCarbon branch(k, 14)*e(branch(k,1)) - pt*e(branch(k,2)); else branchLossCarbon branchLossCarbon (-branch(k,14))*e(branch(k,2)) - pt*e(branch(k,1)); end end理论上totalGenCarbon应该等于totalLoadCarbon加上lineLossCarbon。如果误差超过1%大概率是方向判断写错了。第二层与目标论文逐表对比。把节点碳势、支路碳流率做成表格和目标EI论文里的结果对比。这一步最花时间但也最有效。对比时注意不同论文对发电机排放因子的假设不同碳势数值会有整体偏移这时候要看趋势是否一致——哪些节点碳势高、哪些低高碳潮流从哪个区域流向哪个区域这些结构特征应该是一样的。第三层敏感性分析。把某台发电机的碳排放强度人为调高10%观察全网节点碳势的变化。理论上离该发电机电气距离越近的节点碳势变化越明显电气距离越远的节点变化越小。如果出现该发电机附近节点碳势不变远处节点反而大变的情况说明你的网络关联矩阵组装有误。5. 结果分析怎么做碳势分布、碳流方向与高碳通道识别5.1 结果可视化的两种有效方式碳排放流的计算结果光看数字表格很难形成直观印象。我实践中最好用的可视化方案有两种。第一种是节点碳势热力着色图。把IEEE 14节点的拓扑按标准坐标画出来网上能找到很多现成的坐标文件然后用节点颜色的深浅表示碳势高低。碳势越高颜色越红碳势越低颜色越绿。这张图放在论文里审稿人扫一眼就能看出哪个区域是高碳区。第二种是支路碳流方向箭头图。在每条支路中间画一个箭头箭头方向代表碳流方向箭头粗细代表碳流率大小。这张图能直观展示碳排放的主通道尤其是从煤电大机组向负荷中心输送碳流的那几条关键线路。Matlab里画拓扑图可以用graph对象加plot命令也可以直接用quiver画箭头再叠加text标节点号。我建议先画碳势着色图再在它的基础上叠加碳流箭头一张图把两个维度的信息都表达清楚。5.2 从结果中读出哪些关键信息完成计算和可视化后分析结果时我通常会按以下路径读信息第一高碳节点识别。找出碳势最高的节点通常它就是离大型煤电机组电气距离最近的节点。在改造后的IEEE 14算例中节点1和节点2是煤电挂接点碳势一般在800以上是整个系统的碳势峰值区。第二碳流传输路径。顺着高碳势区域往外看找出碳流率最大的几条支路这些就是系统碳排放的主干道。如果某条支路连接着90MW的负荷节点4或节点5而这条支路又直接承接节点1送出的功率那这条支路的碳流率通常很大值得在分析中单独点名。第三负荷碳流率排名。把14个节点的负荷碳流率从大到小排序找出碳流率最高的负荷节点。这个排序直接反映了谁的用电造成了最多的碳排放是碳责任分摊的核心结论。第四水电的稀释作用。节点8的水电机组虽然容量不大30MW但它的存在会显著降低附近节点的碳势。计算时你会发现节点8附近的节点碳势被稀释到很低这就是零碳电源在碳流分析中的直观体现。论文里讨论清洁能源消纳的碳减排效益时这个计算结果就是最好的论据。5.3 输出图表与论文配图的衔接复现代码最终要服务于论文写作。我的建议是主程序末尾加上自动导出图表的功能%% 绘制节点碳势图 figure; % 这里用你准备好的14节点坐标画拓扑图 plotNodeTopology(coord); hold on; scatter(coord(:, 1), coord(:, 2), 300, e, filled); colorbar; title(IEEE 14节点系统节点碳势分布 (kgCO2/MWh)); %% 绘制支路碳流方向图 figure; plotNodeTopology(coord); hold on; % 用quiver画碳流方向箭头 for k 1:size(branch, 1) f branch(k, 1); t branch(k, 2); mid (coord(f, :) coord(t, :)) / 2; % 根据PF方向决定箭头指向 if branch(k, 14) 0 dirVec coord(t, :) - coord(f, :); else dirVec coord(f, :) - coord(t, :); end dirVec dirVec / norm(dirVec); quiver(mid(1), mid(2), dirVec(1), dirVec(2), 0.1, LineWidth, 2); end导出的图用PDF或SVG格式保存放到论文里做单栏图效果很好。6. 复现过程中躲不开的坑平衡节点、潮流方向与排放因子取值6.1 平衡节点的出力必须用潮流计算后的值这是我最想强调的一个坑。IEEE 14节点系统中节点1是平衡节点它的有功出力不是人为设定的而是潮流计算收敛后自动填补全网功率缺额。很多人在写碳流代码时直接拿case14里gen矩阵的初始出力值比如232.4MW算碳注入这在某些工况下会出错。正确的做法是先跑完runpf然后从result.gen里提取节点1的实际出力。因为潮流计算可能会对平衡节点出力做修正尤其是你手动调整了负荷水平或其它机组出力后如果你用的是修改前的初始值碳平衡校验一定会出问题。细节是魔鬼但魔鬼也怕你把代码流程写规范。6.2 支路潮流方向PF列 vs PT列Matpower的branch矩阵里PF列是首端有功PT列是末端有功。判断支路潮流方向理论上看PF和PT都能判断PF 0或PT 0都表示功率从首端流向末端但在碳流计算中我建议统一使用PF列理由有两点第一PF列是潮流计算收敛后直接得到的支路首端注入功率物理意义最明确。PT列要通过首末端功率平衡才和PF关联如果遇到线损很大的重载支路PT的正负号在临界工况下可能和PF不一致功率从首端流入但在末端由于损耗变为负值这种情况在理论上是可能的虽然正常潮流不会出现。第二统一用PF可以避免代码里出现多种判断逻辑减少出错概率。我之前某版代码里PF和PT混用结果在几条损耗较大的支路上方向判断反了碳平衡校验死活对不上查了半天才发现是方向问题。6.3 多发电机节点与零出力机组的处理IEEE 14节点中每个发电机节点目前只挂一台发电机但实际系统中一个节点挂多台机组很常见。碳流计算时如果节点i有多台发电机要把所有发电机的有功出力求和作为G_ii把所有发电机的碳注入求和作为b_i。另外case14原版数据里节点3、6、8的发电机有功出力为0同步调相机模式。如果你在碳流计算前不手动调整出力这些机组在模型里就是隐形的——不注入功率也不注入碳排放但占用了一个发电机位置。我在代码里用pg 0作为判断条件就是为了忽略这些零出力机组避免矩阵G出现零行导致方程奇异。6.4 碳排放因子取值对结果的影响节点碳势对发电机碳排放因子的敏感度非常高。如果把节点1的煤电排放因子从850改成900不仅节点1的碳势会上升整个下游区域的碳势都会跟着往上走。这意味着如果你要复现某篇论文但没有拿到论文的完整参数表只照着论文里煤电、气电、水电的文字描述去猜排放因子结果大概率对不上。我的建议是复现时优先找目标论文的附录或数据公开部分看有没有给出详细排放因子如果找不到就退而求其次先按行业内通用值跑一版然后在论文或代码注释里明确写出你假设的排放因子。审稿人关注的是方法逻辑只要参数合理并有说明一般不会因为排放因子取值不同而卡你。6.5 完美复现的理性边界最后聊一句顶级EI完美复现这个说法。我复现过几篇EI论文后最大的感受是除非作者完全公开了输入数据和代码否则完美复现基本是一个理论目标。不同论文对发电机排放因子的假设、对网损碳流的处理方式、对平衡节点的处理细节都存在微小差异导致最终数值不可能完全一致。真正的做法是复现出论文的核心结构和趋势——哪些节点碳势高、哪些支路碳流大、负荷碳流率排名如何——这些结构特征一致复现就算成功。不必纠结于小数点后三位是否完全一致那没有意义。6.6 一个隐藏的细节负荷单位与碳流率单位的一致性碳流率计算中有个隐蔽的单位陷阱。如果IEEE 14节点系统的负荷PD单位是MW碳势单位是kgCO2/MWh那负荷碳流率的单位就是kgCO2/h。这个单位本身没问题但当你把负荷碳流率累加起来和发电碳注入对比时一定要统一单位。发电机碳注入的计算是PGMW乘以排放因子kgCO2/MWh得到的结果单位同样是kgCO2/h。两边单位一致碳平衡校验才有效。如果有人习惯把排放因子改成kgCO2/kWh那所有公式都要注意数量级的变化千万不能混用。从14节点到更大系统的扩展思路与最终体会代码在IEEE 14节点系统上跑通之后我顺手把它扩展到了IEEE 30节点和IEEE 118节点系统。核心算法完全不用改只需要换掉loadcase的数据文件再把发电机碳排放强度向量改成对应的长度。14节点系统跑一次的时间不到0.1秒118节点系统也基本是眨眼工夫就出结果这说明碳排放流计算的复杂度主要由节点数和支路数决定对于小规模系统Matlab的矩阵左除完全够用。我在实际做这个课题时最大的体会是碳排放流计算本身的技术门槛不高真正的门槛在于对电力系统潮流的理解——支路方向怎么判断、平衡节点怎么处理、线损碳流怎么算这些细节全部都能追溯到对潮流计算本质的理解。如果你正在准备电力系统方向的论文或者导师给了一个双碳相关的仿真任务我建议你从这套IEEE 14节点的碳排放流代码入手先把手推计算和代码输出对齐一遍把每个指标背后的物理含义弄清楚再考虑扩展到更大的系统。最后分享一个小技巧算完碳流后不妨把每个节点的碳势和这个节点的边际电价做一张散点图。在统一碳价机制下两者通常有正相关关系。这个发现经常能成为论文里一个有意思的讨论点也是碳流分析研究里一个值得深挖的交叉方向。