☰
电力系统碳排放流计算:原理、Matlab实现与工程实践
2026/10/14 7:02:15 网站建设 项目流程

前阵子整理代码库,翻出了一套用Matlab写的电力系统碳排放流计算脚本。最初是因为某项目要核算区域电网里不同负荷节点到底该为多少碳排放“买单”,全网统一用平均排放因子那种做法,算出来的账总觉得不太对劲——同样是用一度电,火电密集区附近取的电和新能源富集区附近取的电,隐含碳怎么可能一样。碳排放流计算就是专门解决这个问题的:把电网里的功率流当作一根根携带碳的“管道”,从发电机侧一路追踪到负荷侧,把碳排放责任按实际电气路径分摊下去。

这篇文章我会把这套计算的完整思路、核心公式、Matlab实现步骤和踩坑记录一次性聊透。代码从潮流输入、碳流追踪矩阵构建到结果可视化都验证过,适合正在做双碳核算、绿电溯源、电网低碳规划的研究生和工程师参考,也适合刚接触这个方向、想理解碳流理论到底在算什么的同学读。

1. 为什么要把碳排放流单独拎出来算

1.1 平均排放因子算不出的那笔“区域差异账”

现在很多碳排放核算还是用省级或区域级的平均排放因子:总发电碳排放除以总发电量,得到单位电量的平均碳排,然后按用电量一乘,完事。这个方法最大的问题在于它假设全网的每一度电都是同质的,不管电从哪来、经过哪些线路、送到哪里去,碳成本完全相同。

但在实际电网里,电源结构差异非常大。一个节点附近可能既有大型煤电基地,又有大规模风电光伏基地,而另一个节点只能靠远距离送电。如果都用同一个平均值,新能源比例高的省份实际上被“平均”拖累了,火电密集区超排的部分却变得不明显。更进一步,在电力市场或绿电交易的环境下,用户想知道自己买的那一部分电量到底对应多少碳排放,平均因子根本给不出这种颗粒度的答案。

碳排放流的核心价值,就是利用现有潮流计算结果,把发电侧的碳排放按支路功率的实际走向分配到每个节点和每条支路上,从而得到空间差异化的碳排放分布。它不是另起炉灶做生命周期评价,而是对电网潮流结果做一次“碳视角”的后处理。

1.2 碳排放流到底在“流”什么

要理解碳排放流,可以先做个类比。电网就像一套水力系统:发电机是水源,水质(碳排放浓度)取决于源头属性;功率从水源出发沿着线路流向负荷,就像水沿着管道流到用户。沿途不同水源的水在某处汇合,混合后的水再分流到不同支路,这时候各支路水的“含碳浓度”应该用混合后的平均浓度来算。

映射到电力系统里,功率从一个节点流向另一个节点时,它把上游发电机产生的碳排放也“携带”过去了。我们把单位功率所携带的碳排放称为碳势,节点碳势由所有注入该节点的功率碳势加权平均得到。然后,流出节点的每一条支路,其碳流量就等于支路有功功率乘以源端节点的碳势。这样追踪下来,每个负荷节点的碳流就是它消费功率乘以所在节点碳势。所有负荷碳流加起来,在忽略网损的情况下恰好等于所有发电机碳排放之和,排放总量守恒。

这套追踪过程本质上是对节点功率平衡方程做了一次碳排放维度的扩展。因为功率平衡方程对每个节点都成立,那么配套的碳平衡方程只要守恒条件满足,也全节点成立,整网计算是可闭合验证的。

1.3 算出一套碳流结果之后,能干什么

数值推出来不是终点,关键是能回答业务问题。

第一类应用是节点碳势核算。对电网里各节点的用电碳排放做差异化定价或者考核,理论上可以把全网的碳排精确到节点级,比宏观平均因子细了不止一个量级。

第二类是支路碳流溯源。如果某条输电线路向高耗能产业园区送电,碳流结果能告诉我们这部分电的碳是从哪个电源基地来的、混合比例是多少。这对绿电凭证发放很有价值,可以识别出“证书与物理电力实现同步”的路径。

第三类是低碳规划。改扩建线路、新增电源、调整机组出力,都会影响全网的碳势分布。碳流计算能给出灵敏度方向,比如在哪个节点接入新能源对整网碳势的稀释效果最明显。很多论文研究低碳经济调度、碳费电价机制,底层都要先过碳流这一关。

2. 整体实现思路与方案选型

2.1 算法路线:先有潮流,再追碳流

碳排放流计算不是独立于潮流计算的新方程体系,而是在已收敛潮流结果基础上追加一个线性计算过程。

完整流程分四步:

  1. 输入电网拓扑与运行数据:节点、支路参数、发电机出力与碳排放强度、负荷。
  2. 执行常规潮流计算,得到各支路有功功率的大小和方向、各节点功率注入。
  3. 根据支路功率方向构建“碳流追踪矩阵”,把发电碳排放作为方程右端项。
  4. 求解线性方程组得到节点碳势,再向后计算出支路碳流与负荷碳流。

这个顺序很重要。碳流计算的一切依据都是潮流结果,如果潮流本身不收敛或者数据方向符号处理错了,碳流结果就会完全失真。所以我在架构上把潮流结果单独保存成结构体,碳流模块只看这个结构体,不直接碰原始数据。这样以后换其他潮流程序算完,把格式统一一下,碳流模块可以原封不动地复用。

2.2 追踪的核心假设:比例共享原则

碳流追踪必须回答一个问题:当多个上游电源的功率汇入同一个节点,再分裂成多个下游支路流出时,下游每一条支路拿到的碳是按什么比例来的?

目前被广泛接受的工程假设叫“比例共享原则”。意思是节点是一个均匀混合池,流出功率的碳势就等于节点碳势。节点碳势由全部注入功率按照功率大小加权平均来的,不会区分这条支路的功率具体来自哪个发电机组。

这个假设物理上并不严格可证,因为交流电网中电子走哪条路径并没有唯一性,这是一个病态的物理路径辨识问题。但工程上它是最公平合理的折中方案:每个节点都按混合平均来处理,既简单又满足常识预期。类似的“追踪”思想还用在输配电成本分摊、网损分摊等很多领域,不是碳流独创的。

在实现时我用矩阵形式,把比例共享原则写成线性方程组。这是整个代码的数学骨架。

2.3 为什么用Matlab实现

选Matlab做这个方向,有几个很实际的理由。

一是建模和调试快。碳流计算核心就是稀疏矩阵的构建和一次线性方程组求解,Matlab对这两种操作支持非常成熟。相比C++,省掉了大量内存管理和矩阵索引的繁琐工作,可以快速迭代算法逻辑。

二是电力系统领域生态好。很多电网算例和潮流工具都有Matlab版本,尤其大家熟悉的IEEE标准测试系统,在Matlab里能直接加载,不用手工敲数据。数据格式统一,对比别人的结果非常方便。

三是可视化顺手。碳势柱状图、支路碳流分布图、时序碳势曲线,Matlab的绘图命令直接出图,对解释计算结果很友好。后期写报告或者论文插图都不用再导给其他画图工具。

当然它的短板也很明显,大规模系统性能不如编译型语言。但我后面会讲怎么用稀疏矩阵和预分解来弥补,几千节点的系统在普通电脑上跑起来也没压力。

2.4 程序模块怎么设计才不混乱

代码如果全部堆在一个脚本里,第一眼能看,改几版之后就想扔了。我最终把它拆成六个相对独立的部分,每个文件职责明确:

  • 主入口脚本:控制整个流程,负责调用下面的模块。
  • 数据预处理函数:从标准算例或外部文件读取并整理节点、支路、发电数据。
  • 潮流计算封装:调用潮流求解器,统一提取支路有功大小与方向。
  • 碳流核心计算函数:输入潮流结果,输出节点碳势、支路碳流、负荷碳流。
  • 正确性校验函数:做碳排放守恒测试和同质化短路测试。
  • 可视化函数:绘制节点碳势与碳流分布图。

这样的好处是,以后想替换潮流求解器、想接入新的电源碳强度数据,只需要改一小部分函数,核心计算逻辑完全不用动。

3. 核心公式与Matlab实现细节

3.1 节点碳势方程的推导

设网络有 n 个节点,K 条支路。用 \rho_i 表示节点 i 的碳势,单位是吨二氧化碳每兆瓦时(tCO2/MWh)。设节点 i 上的机组集合为 G_i,第 g 台机组有功出力为 P_{g},碳排放强度为 e_{g}。同时,设从节点 j 流向节点 i 的支路有功功率为 P_{ji},如果某条支路实际流向与编号方向相反,则 P_{ji} 取负值并归到对应拓扑关系中。

节点 i 的碳平衡方程写成:

[ \left(\sum_{g \in G_i} P_g + \sum_{j \ne i} P_{ji}\right) \rho_i = \sum_{g \in G_i} e_g P_g + \sum_{j \ne i} \rho_j P_{ji} ]

左边括号内是节点 i 的总注入有功功率,它乘以节点碳势表示节点上所有“功率流”携带的碳总量;右边第一项是机组直接注入的碳排放,第二项是上游支路功率携带进来的碳流。

整理一下,把所有 \rho 项移到同一侧:

[ \left(\sum_{g \in G_i} P_g + \sum_{j \ne i} P_{ji}\right) \rho_i

  • \sum_{j \ne i} \rho_j P_{ji} = \sum_{g \in G_i} e_g P_g ]

这就是全系统的线性方程组:A \rho = b。其中 b 向量每个元素是对应节点的机组碳排放注入总量。

把所有节点的方程组合起来,系数矩阵 A 的对角线是节点总注入功率,非对角元素就是“从其他节点流入本节点功率的相反数”。只要存在至少一个平衡节点且网络连通,方程组通常可解。

3.2 方向判断与功率流入矩阵构建

这是代码里最容易出错的一步。潮流结果里每条支路给出的首端有功功率有正有负,正负号取决于你在数据中定义的支路起点和终点。如果不先统一方向判断,后面的矩阵符号就乱套了。

我的做法是:先用潮流结果取出每条支路首端有功功率 Pf,单位标幺值,再乘以基准容量转成兆瓦。然后按正负号判断真实流向:

  • Pf 为正:功率从支路起点 f 流向终点 t。
  • Pf 为负:功率实际从终点 t 流向起点 f,大小是 abs(Pf)。

用一个 n 行 n 列的稀疏矩阵 Win 来记录真实流向,Win(from, to) 存从 from 流向 to 的功率。遍历所有支路,按上述规则填充。构建完成后,对每个节点 i,sum(Win(:, i)) 就是所有线路流入该节点的总功率,sum(Win(i, :)) 是流出该节点的总功率。两者与节点注入、负荷之间应该满足功率平衡:

[ \sum_{j} \text{Win}(i,j) - \sum_{j} \text{Win}(j,i) + P_{D,i} - P_{G,i} \approx 0 ]

我一般把这条校验写进代码,如果偏差超过 1 兆瓦就报警。这会救你很多次。

核心代码片段如下:

function [rho, branch_carbon, load_carbon] = calc_carbon(...) % 输入: branch为nb*2的支路起点终点编号数组, Pf为nb*1的首端有功功率(MW) % P_G_bus为n*1发电出力, C_G_bus为n*1机组碳排放注入, P_load为n*1负荷 n = size(P_G_bus, 1); Win = sparse(n, n); % 用稀疏矩阵省内存 for k = 1:length(Pf) f = branch(k, 1); t = branch(k, 2); p = Pf(k); if p > 1e-6 Win(f, t) = Win(f, t) + p; elseif p < -1e-6 Win(t, f) = Win(t, f) - p; end end % 节点总流入功率 = 机组出力 + 线路流入 E_inflow = P_G_bus + sum(Win, 1)'; % 构造矩阵: A = diag(E_inflow) - Win' A_mat = diag(E_inflow) - Win'; b_vec = C_G_bus; % 求解节点碳势 rho = A_mat \ b_vec;

注意这里 A_mat 用了 Win 的转置,不是 Win 本身。因为 \rho_i 方程中要减去的是“流入 i 的功率乘以上游节点碳势”,流入 i 的功率在 Win 里是 Win(j, i),转置后正好落在 A_mat(i, j) 位置。我第一次写代码时在这里被绕晕,删掉转置直接跑,结果碳势一大堆负值,怎么调都不对。

3.3 支路碳流与负荷碳流计算

解出节点碳势后,后续计算就非常直接了。每条支路的碳流等于支路功率乘以其送端节点的碳势。这里送端节点是功率实际流出的节点,同样需要用方向判断结果来确定。

nb = size(branch, 1); branch_carbon = zeros(nb, 1); for k = 1:nb f = branch(k, 1); t = branch(k, 2); p = Pf(k); if p >= 0 branch_carbon(k) = p * rho(f); else branch_carbon(k) = (-p) * rho(t); end end % 负荷碳流:节点负荷乘以该节点碳势 load_carbon = P_load .* rho;

这样算出来的负荷碳流,是节点负荷对应的碳排放责任。如果节点上有多个负荷,可按功率比例再细分。

这里需要说明一个处理细节:网损。实际线路有损耗,功率在线路两端并不守恒,碳流到下游支路时实际上有一点功变成了热消耗,这部分碳到底算谁的,理论上是有争议的。工程上最常用的处理是忽略网损,因为线路损耗通常占传输功率的几个百分点以内,对节点碳势的数值影响非常小。如果实在要求严谨,可以在全网守恒校验时把线路损耗作为一项单独显示,不做分摊。

3.4 正确性校验:两个快速测试

代码写完不是直接交差的,必须做校验。我每次改完代码必跑两个测试。

第一个叫同质机组测试。把所有发电机的碳排放强度都设成同一个值,比如 0.5 tCO2/MWh。此时无论网络多复杂,所有节点的碳势必须严格等于 0.5。如果算出来不是,基本可以断定矩阵构造或者方向判断有问题。这个测试一票否决,比什么都好使。

第二个是碳排放守恒测试。在忽略网损的前提下,全网所有机组的碳排放总和应该约等于所有负荷碳流之和。写成代码就是 sum(C_G_bus) 与 sum(load_carbon) 做差,偏差控制在1%以内。如果对不上,大概率是 Win 矩阵漏了支路或者方向判断阈值设得太大。

我帖一下校验代码的写法:

total_emission = sum(C_G_bus); total_load_carbon = sum(load_carbon); dev = abs(total_emission - total_load_carbon) / total_emission; if dev > 0.01 warning('碳守恒偏差过大: %.4f%%', dev * 100); else fprintf('碳守恒校验通过, 偏差: %.4f%%\n', dev * 100); end

这两个测试加在一起,能覆盖绝大多数低级错误,强烈建议保留在流程里。

4. 实操演示:在14节点标准算例上跑通全流程

4.1 输入数据准备与字段约定

碳流计算的输入数据分成四块:节点数据、支路数据、发电机数据、碳排放强度。以大家常用的14节点标准算例为例,它的数据格式是这样的:

  • 节点数据 bus:每行一个节点,必须包含节点编号、有功负荷、无功负荷等基础信息。
  • 支路数据 branch:每行一条支路,必须包含起点节点编号、终点节点编号、电阻、电抗、对地电纳等,支路编号不用节点编号替代,单独按顺序排列。
  • 发电机数据 gen:每行一台机组,包含所在节点编号、有功出力、无功出力、运行上下限等。
  • 碳排放强度表:每台机组对应一个数值,单位 tCO2/MWh。这个数据外部输入,不存在标准算例里。

如果你用Matlab电力系统领域的公共工具箱加载算例,那么节点、支路、发电机数据都能自动读进来,唯一需要额外准备的就是每台机组对应的碳排放强度。可以做成 Excel 表,第二列填机组对应的节点编号,第三列填排放强度,用的时候读进去并跟 gen 数据对齐。

我给个小建议:所有输入数据统一在入口处做一次完整性检查。检查节点编号是否越界、支路两端编号是否存在、机组所在节点是否有效、每条支路阻抗是否为零(会出现除以零的数值问题)。这些细节看起来笨拙,却能直接把问题拦在计算之前,省得后面排查矩阵异常时一头雾水。

4.2 一个可直接运行的主流程代码

下面这段是去掉详细注释后仍然可直接跑的全流程代码结构,把前面几节的内容串在一起。假设潮流结果已经由某个计算函数得到,或者你能从公开标准算例的结果文件里提取到支路首端功率。

%% 碳流计算主流程 clc; clear; close all; % 1. 读入数据(以公共算例工具为例) mpc = standard_case_14(); % 返回bus, branch, gen等结构 bus_data = mpc.bus; branch_data = mpc.branch; gen_data = mpc.gen; baseMVA = mpc.baseMVA; % 2. 读取外部碳强度表:每行[节点编号, 排放强度tCO2/MWh] emission_table = load_emission_table('emission.csv'); n = size(bus_data, 1); % 3. 计算潮流并提取支路首端有功功率 result = run_power_flow(mpc); % 返回含branch功率的结果 Sf = result.branch(:, 14); % 支路首端复功率标幺值 Pf = real(Sf) * baseMVA; % 转成兆瓦 % 4. 构造节点发电出力与碳排放注入 P_G_bus = zeros(n, 1); C_G_bus = zeros(n, 1); for k = 1:size(gen_data, 1) bid = gen_data(k, 1); Pg = gen_data(k, 2); % 发电机有功出力, 单位MW e_g = emission_table(find(emission_table(:,1)==bid), 2); P_G_bus(bid) = P_G_bus(bid) + Pg; C_G_bus(bid) = C_G_bus(bid) + Pg * e_g; end P_load = bus_data(:, 3); % 有功负荷 MW % 5. 调用碳流核心计算函数(见3.2与3.3) [rho, branch_carbon, load_carbon] = calc_carbon(... branch_data(:, [1 2]), Pf, P_G_bus, C_G_bus, P_load); % 6. 输出关键结果到控制台 for i = 1:n fprintf('节点 %3d:碳势 %.4f tCO2/MWh,负荷碳流 %.2f tCO2/h\n', ... bus_data(i,1), rho(i), load_carbon(i)); end

如果你没有潮流工具,也可以用网上的直流潮流计算替代。直流潮流只算有功功率,对碳流计算来说通常足够,因为碳流关心的就是有功潮流方向与大小。只要把直流潮流算出来的支路有功功率代入 Pf,后面的碳流计算流程完全不用改。

4.3 计算结果怎么解读

以14节点标准算例为例,我拿一套参考数据跑了一次,得到的结果非常有规律性。发电机接入节点的碳势直接等于该节点机组的碳排放强度,因为初始注入是没有被上游“污染”的。随着功率往外输送,节点碳势会逐渐向全网机组的加权平均水平靠拢,但因为不同节点从不同机组取电的比例不同,节点碳势仍然有明显差异。

靠近高排放煤电节点的母线碳势明显偏高,可能达到 0.8 tCO2/MWh 上下;靠近水电或新能源节点的母线碳势就低一些,可能只有 0.4 左右。负荷端的碳流受碳势和负荷大小两个因素影响,一个大功率高碳势节点可能贡献全系统接近两成的负荷碳排放。

读结果的时候要注意一点:节点碳势反映的是“该节点功率的边际平均碳含量”,不代表某个具体用户在物理上真的用了某台特定机组的电。碳流是核算方法,不是物理测量,它的价值是提供一种合理且可复现的责任分摊方式。

可视化我做了三张图:

  • 节点碳势柱状图:一眼看出全网哪些节点低碳、哪些高碳。
  • 负荷碳流占比饼图:体现主要用电节点的碳排放责任。
  • 系统总排放与总负荷碳流对比表:做守恒性展示。

可视化代码比较简单,用 bar 和 pie 就能完成。重点是把节点编号标注清楚,结果图直接贴在报告里别人能看懂。

5. 常见问题、坑点排查与性能优化

5.1 矩阵奇异、潮流不收敛与碳势NaN

碳流矩阵 A 奇异最常见的原因是网络中存在“孤立节点”或“零注入节点”。比如某个节点与全网其他节点没有任何支路相连,或者连接支路停运了,它的 E_inflow 为零。此时对应方程变成 0 = 0,矩阵必然奇异,Matlab 中会出现一条警告并解出 NaN。

遇到这种情况,先不要急着怀疑算法,用这几步排查:

  1. 把潮流结果里每一条支路的状态列拿出来,确认没有把停运支路纳入碳流计算。
  2. 检查 Win 矩阵每行每列的总和,找出总流入或总流出都为零的节点。
  3. 对全网络做一次连通性检查。用图论工具箱里的连通分量函数,看是不是有节点脱离主网。

还有一种常见情况是潮流本身没收敛,但代码忽略了这个错误继续算。如果你用的潮流函数返回值里含收敛标志位,一定要判断它。潮流不收敛时支路功率可能是噪声,碳流计算结果自然毫无意义,而且经常表现为负数。

if ~converged error('潮流未收敛,请先检查输入数据或网络参数'); end

5.2 碳势出现负值,突破口在方向判断

节点碳势出现负值,理论上完全不合理。我调试初期遇到最多的问题就出在这里,而且几乎每次都是方向判断出错。

最典型的坑有两个。第一个是把潮流结果的首端功率当作流动方向时,忘记处理负号。一条支路如果首端功率为负数,说明它实际从末端节点流向首端节点。如果在构建 Win 矩阵时还是往 f 到 t 方向填负值,矩阵里就会出现负的流入功率,算出的碳势自然乱套。第二个坑是把阈值设得太大。我见过有人把方向判断阈值设为 0.1 兆瓦,结果正常运行功率小的支路被当成零功率丢弃,节点碳势出现极大跳变。

我的经验是阈值取 1e-4 左右即可,既能过滤数值噪声,又不会漏掉实际小额功率。方向判断完以后,一定在调试阶段打印每个节点的流入功率明细,肉眼扫一遍,比所有自动排查都直观。

5.3 处理网损时怎么取舍

前面提到我默认忽略网损,但很多初学者会在这一步反复纠结。网损对碳势的数值影响通常在千分之一级别,对负荷碳流的影响也有限,但在做总量守恒校验时会把“不平衡量”全归到网损上。

如果项目要求必须考虑网损,有一个相对简单的做法:先把各支路损耗功率算出来,然后将损耗功率按比例分摊到该支路两端节点的负荷碳流上,或者引入一个虚拟负荷节点来吸收损耗功率。更精细的方法是把网损当作支路终点侧的功率消耗来迭代计算节点碳势,但因为碳势本身跟着潮流路径走,迭代计算量不大,收敛性也不错。

不过在大多数场景下,经典论文和工程实务都是直接忽略网损的,不用觉得这是重大缺陷。

5.4 让计算跑得更快的小技巧

碳流计算本身规模很小,但如果你要做时序碳流分析,比如8760小时的全年逐小时碳势曲线,计算速度就很重要了。

我常用的优化手段有三个:

  1. 用 sparse 矩阵而不是全矩阵。几百节点时全矩阵可能还看不出差距,上千节点时全矩阵的内存开销会造成很大压力。
  2. 如果节点矩阵结构不随小时变化,可预先把 A_mat 做 LU 分解,然后每一小时只需要回代一次,速度能提升很多倍。
  3. 避免在循环里重复构建 Win 矩阵。时序计算中如果网络拓扑不变,可以提前把每条支路是否参与Win构建的索引算好,循环里只更新功率数值。

以下是一个预分解加速的示意:

[L, U] = lu(A_mat); % 仅分解一次 % 在每小时循环内 rho_h = U \ (L \ b_vec_h); % 每次只需两次回代

这套做法我实际用在一个包含几百节点的区域电网碳排放时序分析上,全年8760个时段的碳势计算大概几十秒跑完,性能完全够用。

5.5 后续可以怎么扩展

碳流计算的框架搭好以后,扩展方向非常多。

一是把潮流结果换成时序断面,做全年的节点碳势曲线。这能反映新能源出力波动对碳排放的空间影响,比如白天光伏大发时节点碳势压低,夜间负荷高峰时煤电顶上来碳势升高。

二是把碳排放强度和发电成本一起纳入优化模型。做低碳经济调度的时候,目标函数里同时考虑煤耗成本和碳排成本,节点碳势作为碳价格信号的替代指标,可以指导机组出力和负荷响应策略。

三是扩展到区域间碳流追踪。把电网边界等效为虚拟支路,跨区传输功率携带的碳流就能在区域层面做量化核算,这个在国家层面的碳市场、省间交易中很有价值。

我个人的体会是,碳排放流计算真正的价值不在那几行矩阵求解,而在于它逼你把“电从哪里来”这件事想清楚。平均因子的时代迟早会过去,空间差异化的碳核算会越来越成为标准配置。代码本身不难,难的是把每个节点的功率方向、每条支路的碳流责任都梳理得明明白白。你把这套原理和代码跑通一次,以后再碰到碳核算相关的问题,思路会清晰很多。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询