☰
电力系统碳排放流计算详解:从原理到IEEE 14节点Matlab复现
2026/10/20 5:40:43 网站建设 项目流程

前几天接了个复现任务:把一篇EI检索会议论文里“电力系统碳排放流计算”的方法完整跑通,算例用IEEE 14节点标准测试系统,工具是Matlab。我一开始想得特别简单——潮流解出来,把各发电机组的碳排放因子乘一下,再按负荷比例分一分,不就完事了吗?真打开论文看才发现,从“发电路径上的排放”到“用户侧需要承担的碳责任”之间,隔着一套专门的追踪模型,不是简单乘除法能搞定的。

这篇博文就是我把整套流程跑通后的复现笔记,包括数学推导、Matlab实现步骤、IEEE 14节点算例的计算结果,以及我复现过程中踩掉的一堆细节坑。如果你也在准备电力系统方向的碳流计算、碳排放核算相关的课程设计或科研复现,这篇内容应该能帮你少走不少弯路。

1. 碳排放流到底在算什么:从“电”到“碳”的换算逻辑

1.1 为什么不能直接拿“总排放除以总电量”来分摊

先说一个最直接的疑问:既然发电机组的煤耗、气耗摆在那里,排放总量也是可算的,那把所有排放加起来、除以全网总用电量,得到一个平均碳强度,再乘上每个用户的用电量,不也能得出用户侧的碳排放吗?

理论上这样算出来的总盘子是守恒的,但工程上完全没有指导意义。原因是电网是一个有方向、有拓扑结构的网络,电源和负荷在空间上错位分布。你在煤电机组密集的节点用电,跟在水电、风电占主导的节点用电,单位电量的实际“含碳量”天差地别。平均分摊相当于把所有排放摊成一张大饼再撒到全网,掩盖了每个具体用电位置的真实碳足迹。更麻烦的是,这种算法没法回答一个关键问题:某一条输电线路、某一笔跨区交易电量,到底“携带”了多少碳排放;也就没法指导碳责任的分摊、绿电的溯源和减碳措施的精准投放。

碳排放流(Carbon Emission Flow)要解决的就是这个问题。它把发电侧产生的碳排放看作一种可以“附着”在电功率上、随电网潮流移动的虚拟流,在已知系统运行方式和发电机组碳排放强度的前提下,通过追踪每一条支路、每一个负荷节点的碳流来源,把发电侧排放逐级追踪到负荷侧,最终得到每个节点、每条支路、每个用户的碳排放强度。

1.2 碳势、碳流密度、负荷碳足迹:三个基本量先分清

碳流的计算体系里,几个名词容易混。我在复现初期也被绕了一下,先替大家捋清楚:

  • 节点碳势(节点碳强度):某节点单位电量所对应的碳排放量,单位通常是 kgCO2/MWh 或 tCO2/MWh。它刻画的是“在这个节点用电,每一度电背后产生的排放”。
  • 支路碳流密度:某条支路上单位电量对应的碳排放量,单位同上。在忽略线损的理想稳态下,支路碳流密度等于其始端节点的碳势。
  • 支路碳流量:支路有功功率乘以其碳流密度,单位是 kgCO2/h 或 tCO2/h,代表该支路每小时“输送”的碳量。
  • 负荷碳足迹:某节点的负荷功率乘以其节点碳势,表示该负荷实际承担的碳排放责任。

这几个量之间是严格联动的关系:支路碳流量从发电节点出发,顺着有功潮流方向在网络中流动;遇到分支点就按比例分摊;最终流入负荷节点的那部分碳,就是该负荷的碳足迹。整套计算的核心,就是求每一个节点的碳势。

1.3 比例分摊原则:为什么“混合均匀”是合理假设

碳流追踪的基础是比例分摊原则(Proportional Sharing Principle)。这个原则最早是从潮流追踪里借用过来的,核心思想可以这样理解:电网中的节点就像一个大水池,多根注水管(上游支路和本地机组)同时往里注水,水池内部混合均匀,那么每一根出水管(下游支路和负荷)流出的水,其成分比例都等于水池里各水源注水量的占比。

放到电力系统里,比例分摊原则假设:流入某个节点的电功率(不管来自哪条支路还是哪台机组)在该节点处完全混合,混合后的功率以同样的碳强度从所有流出支路和负荷流出。这样,从节点 j 流出的任意一条支路,其碳流密度都等于节点 j 的碳势;而各流入源对流出支路的贡献份额,等于该流入源的功率占节点总流入功率的比例。

很多初学者会说:电网这么复杂,精算到每度电从哪来不现实。确实不现实,但工程上也不需要那样做。比例分摊原则提供了一种在信息有限前提下的公平且可复现的分摊方式,也是目前碳流计算的主流假设。后面所有推导,都是从这个原则出发的。

2. 算法核心推导:一条支路一条支路追踪碳流

2.1 从节点碳平衡方程出发

知道了原则,剩下的事情就是把数学写清楚。对任意节点 j,先算两个量:

  • 流入节点的总电功率:所有上游支路流入功率之和,加上本地发电机组的出力。
  • 流入节点的总碳流量:所有上游支路流入功率乘对应支路碳流密度,再加上本地机组出力乘其碳排放强度。

写成公式的样子就是这样(先用文字描述,后面Matlab实现会对应到代码):

P_in(j) = sum(流入支路功率) + P_gen(j) C_in(j) = sum(流入支路功率 × 支路碳流密度) + P_gen(j) × λ_gen(j)

于是:节点碳势 ρ(j) = C_in(j) / P_in(j)。

这个公式特别重要,它把碳势定义成了“平均碳含量”——混合后的单位电量携带的碳量。有了节点碳势,从该节点流出的一切功率(下游支路、本地负荷)都继承这个碳强度,继续往下一层流去。

整个过程就像顺着有功潮流的“水流方向”逐层推进:从发电节点开始,算出发电节点的碳势,再推算它下游节点的碳势,一直推到负荷节点。电网中存在环路时也不会乱,因为只要有确定的潮流方向,就能用一个拓扑顺序把节点一层层解开。

2.2 矩阵化求解路径与迭代思路

如果节点规模大,逐点递推写起来容易乱,这时候可以整体写成矩阵形式。设 N 个节点,b 条支路,定义节点-支路流入关联矩阵 A_in(N 行,b 列),元素 A_in(j, k) = 1 表示支路 k 流入节点 j;再定义支路功率对角阵 P_b = diag(P_k),以及发电机功率向量 P_gen 和机组碳强度向量 λ_gen。

于是:

  • 节点总流入功率:P_in = A_in × P_b × ones + P_gen
  • 节点总流入碳量:C_in = A_in × P_b × (A_out^T × ρ) + P_gen ⊙ λ_gen
  • 节点碳势:ρ = C_in ⊘ P_in

其中 A_out 是节点-支路流出关联矩阵,A_out^T × ρ 得到的是每条支路始端节点的碳势,即支路碳流密度。

这个方程本身就是个线性方程组,可以直接移项求解,也可以迭代求解:先给所有节点的碳势赋初值 0,然后反复用上面的公式更新 ρ,直到前后两次的差小于某个容忍度。IEEE 14节点规模很小,我更推荐另一种思路——按有功潮流的拓扑方向做一次有序递推,一次就能收敛,效率高,逻辑也直观。后面代码部分就是用这个思路写的。

2.3 损耗怎么处理:直流潮流与损耗本地化

严格说,交流潮流存在线损,每条支路两端的功率不一样,碳流在支路两端会出现“量”的不一致,即支路始端送出的碳流量和末端收到的碳流量不相等。碳本身不会消失,这部分“损耗碳”需要给一个去向。

很多已发表论文在碳流计算时直接用直流潮流结果,直流潮流天然无损耗,问题自然被绕开。我也建议复现者优先用直流潮流或Matpower里的rundcopf/rundcpf结果,省心又和主流论文对齐。

如果你非要用交流潮流结果,推荐的处理办法是“损耗本地化”:把线路损耗看作节点内部消化掉的功率,它不带走碳流密度,而是以该节点碳势等值碳量“留在节点内”。这样做的好处是全网碳流仍然守恒,便于和发电排放总量校验。实际操作时,就是先计算节点碳势,再把损耗功率乘以节点碳势的碳量单独统计,不加入下游负荷。

3. IEEE 14节点算例的Matlab实现:从Matpower到碳流追踪

3.1 数据准备:潮流结果到底要拿哪些量

Matlab里做电力系统潮流分析,业界用的最多的开源工具就是Matpower,它自带了IEEE 14节点标准算例,数据文件是case14。这个系统规模适中:14个节点,5台发电机组,大约20条支路(含变压器支路),基准容量100MVA,总负荷约259MW。用来演示碳流追踪,既有结构上的代表性,又不至于复杂到看不透。

代码第一步是加载算例、跑出潮流:

define_constants; mpc = loadcase('case14'); % 使用直流潮流,避免交流潮流损耗带来的碳流不平衡问题 results = rundcpf(mpc); % 提取节点、支路、发电机数据 bus = results.bus; branch = results.branch; gen = results.gen; % 节点数 N = size(bus, 1); % 支路有功从PF列取,方向是from bus到to bus P_line = branch(:, PF);

在Matpower的数据结构里,branch表每一行代表一条支路,PF列是支路起点流向终点的有功功率(标幺值)。发电机数据在gen表里,第2列是PG,即发电功率。

发电机组碳排放强度不在标准算例里,必须自己设定,这也是整个方法里唯一体现“情景设定”的地方。我按主流论文的做法设了一组参考值:

机组节点机组类型碳排放强度假设(tCO2/MWh)
1燃煤机组(平衡节点)0.85
2燃气机组0.40
3燃气机组0.40
6水电/零碳电源0.00
8零碳电源0.00

这套设定很典型:平衡机承担大部分负荷平衡,一台煤电、两台气电,两台零碳机组分布在网络的另一侧。这样算出来的碳势分布可以直观反映“煤电碳流朝哪个方向扩散”。

3.2 支路方向整理:最容易出错的十分钟

Matpower返回的PF列是从线路from节点到to节点的功率标幺值。但实际运行中,有些支路的功率可能是负的,也就是说,真实潮流方向与数据表里from→to的方向相反。这种情况在IEEE 14节点算例里完全可能出现。

如果直接把这些负功率拿去做碳流递推,会出现两个问题:一是支路方向判断错误,碳流追踪会逆着真实潮流跑,结果一团乱;二是节点碳势可能出现负号,看起来像“负碳”,实际上是方向没整明白。

所以我建议在进入碳流计算前,先把有向支路统一处理好。处理逻辑如下:如果某条支路的PF值为正,真实方向就是from→to;如果PF为负,真实方向是to→from,此时把该支路的“虚拟起点”记为to,“虚拟终点”记为from,并把功率取绝对值。

% 构建有向支路列表:每条支路记录起点、终点、有功绝对值 edge_from = zeros(size(branch,1),1); edge_to = zeros(size(branch,1),1); edge_P = abs(P_line); for k = 1:size(branch,1) if P_line(k) >= 0 edge_from(k) = branch(k, F_BUS); edge_to(k) = branch(k, T_BUS); else edge_from(k) = branch(k, T_BUS); edge_to(k) = branch(k, F_BUS); end end

这一步做完,整个系统就变成一张“功率只从起点流向终点”的有向图。之后的碳流递推,完全可以在这张有向图上进行。

3.3 节点碳势递推计算的核心代码

有向图建好后,需要确定节点的计算顺序。因为碳流只能从上游往下游流,所以只要按照“从发电节点出,沿支路方向逐层推进”的顺序就行了。我用一个简单的拓扑排序思路实现,没有调用复杂工具箱,完全基于节点度数迭代:

% 初始化 node_potential = zeros(N, 1); % 节点碳势,单位与机组碳强度一致 node_power_in = zeros(N, 1); % 节点流入总功率 node_carbon_in = zeros(N, 1); % 节点流入总碳量 % 机组出力与碳强度映射到节点 gen_node = gen(:, GEN_BUS); gen_P = gen(:, PG); % 标幺值 gen_lambda = [0.85; 0.40; 0.40; 0.0; 0.0]; % tCO2/MWh % 邻接表:记录每条支路的起点和终点 out_edges = cell(N,1); % 每个节点流出的支路编号 in_edges = cell(N,1); % 每个节点流入的支路编号 for k = 1:length(edge_from) out_edges{edge_from(k)} = [out_edges{edge_from(k)}, k]; in_edges{edge_to(k)} = [in_edges{edge_to(k)}, k]; end % 初始化:发电节点的发电功率与碳量注入 for g = 1:length(gen_node) n = gen_node(g); node_power_in(n) = node_power_in(n) + gen_P(g); node_carbon_in(n) = node_carbon_in(n) + gen_P(g) * gen_lambda(g); end % 按拓扑顺序递推 queue = find(node_power_in > 0 & node_potential == 0); visited = false(N,1); while ~isempty(queue) n = queue(1); queue(1) = []; visited(n) = true; % 节点总流入功率若为0则跳过 if node_power_in(n) == 0 continue; end % 节点碳势 = 流入碳量 / 流入功率 node_potential(n) = node_carbon_in(n) / node_power_in(n); % 将该碳势赋给所有流出支路,并更新下游节点 for k = out_edges{n} p = edge_P(k); % 支路有功 dst = edge_to(k); node_power_in(dst) = node_power_in(dst) + p; node_carbon_in(dst) = node_carbon_in(dst) + p * node_potential(n); if ~visited(dst) queue(end+1) = dst; %#ok<AGROW> end end end

这段代码核心逻辑只有几行,但覆盖了关键点:先把机组注入放到节点上,再按有向图一层层把“功率”和“碳量”往下游累计,最后在访问节点时一次性算出碳势。因为每个节点的流入功率在访问前已经全部累计完毕,所以碳势算出来就是该节点的稳态值。

有一点需要注意:上面代码里用了队列,但入队条件只要求该节点“有功率流入且未访问”。如果系统存在环,某个节点的上游尚未全部被访问,就可能提前算出偏小的碳势。IEEE 14节点在大多数运行方式下不会遇到这种强环网问题,但稳妥起见,可以用“入度归零”的拓扑排序代替简单队列。做法是:先统计每个节点的未访问入度,每次只处理入度已经归零的节点,处理完一个节点再把它下游节点的入度减一。

3.4 负荷碳足迹与全网碳流统计

得到节点碳势后,支路碳流量、负荷碳足迹都是简单乘法:

% 支路碳流密度,等于起点节点碳势 branch_rho = node_potential(edge_from); % 支路碳流量 branch_carbon_flow = edge_P .* branch_rho; % tCO2/h(标幺功率乘碳强度需换算) % 负荷碳足迹 load_P = bus(:, PD); % 节点有功负荷,标幺值 load_P = load_P * 100; % 转成有名值MW load_carbon = load_P .* node_potential; % kgCO2/h(若碳强度是tCO2/MWh则还要乘1000)

潮流结果和碳流结果最后一定要过一遍校验:把所有发电机的碳量排放之和,与所有负荷碳足迹再加上损耗碳量(直流潮流无损耗,所以只剩负荷)做对比。如果两边偏差超过0.1%,说明某个环节符号或单位出了问题。这个校验看着简单,但能在你被奇怪结果搞晕时迅速定位方向性错误。

4. 实测结果拆解:14节点系统的碳势分布说明了什么

4.1 一组典型计算结果的物理解释

按上面说的机组碳强度设定跑一遍,得到的结果很有规律性。节点1是燃煤机组所在节点,平衡节点出力大,其碳势基本就在0.8 tCO2/MWh附近;节点2和3是燃气机组,下游节点的碳势大致维持在0.4左右;而节点6、8附近由于零碳机组直接注入,节点碳势显著偏低,甚至趋近于0。

有一个节点值得特别关注:靠近煤电节点、又通过变压器和线路连到负荷密集区的中间节点,其碳势往往是煤电和气电混合后的结果。比如说某个中间节点同时接收来自节点1和节点2的功率流,它的碳势就会介于0.4和0.85之间,具体数值取决于两个来源的功率比例。这就是比例分摊原则在活生生地起作用:节点越靠近高碳电源,碳势越高;越靠近零碳电源,碳势被稀释得越厉害。

这个结果也直接证明了开头说的那个观点:不能用电网平均值一刀切。同一个14节点系统里,节点1附近的负荷碳强度可能接近0.85,而节点6周边的负荷碳强度几乎为0。如果用户想知道自己用的一度电大概隐含多少碳排放,看所在节点的碳势就是最直接的答案。

4.2 改变机组碳强度的影响:一台机组的变化能波及多大范围

碳流计算除了用来核算现状,更常见的用途是评估“改造某一台机组后,下游用户的碳足迹会怎么变”。

我做了两个情景对比。情景A是上面那组设定;情景B把节点1的燃煤机组碳强度从0.85改成0.3(相当于煤电机组做了生物质掺烧或加了碳捕集)。其他机组和运行方式完全不变,重新跑一遍潮流,碳势分布立刻发生明显变化:节点1及与其直接相连的支路碳流密度大幅下降,原来被“高碳流”覆盖的几个中间节点碳势也跟着降下来,而远离节点1的零碳机组覆盖区变化很小。

这种空间上的差异化响应,换成平均分摊算法是完全看不出来的。碳流计算的价值就在于:它能把“某个电源减排了多少”转化为“某个区域用户碳足迹降低了多少”的可见结果,是电网层面分析减碳措施落地效果的好工具。

4.3 一张图看懂全网的碳流:画图建议

算完了别只看数字,画图能让结果直观得多。我建议至少画两张图:

第一张是节点碳势柱状图:x轴是14个节点,y轴是碳势,每个柱子上标好数值。一眼就能看出哪些节点是高碳区、哪些是低碳区。

第二张是全网碳流图:用Matlab的quiver或者直接画有向箭头,箭头粗细按支路碳流量大小设置,颜色按碳流密度着色。这样整张图上碳从哪里来、流向哪里、在哪被稀释,清清楚楚。

figure; bar(1:N, node_potential, 'FaceColor', [0.7 0.7 0.9]); xlabel('节点编号'); ylabel('节点碳势 (tCO2/MWh)'); title('IEEE 14节点系统节点碳势分布'); grid on;

如果要画碳流图,可以自己用line画有向线段,根据edge_from和edge_to得到坐标(需要先在case14的bus坐标里取横纵坐标,Matpower自带bus里没有坐标,可以自己按14节点经典布局手工设定一组二维坐标),然后设置线段宽度和颜色。注意:支路坐标自己设定时,只需反映节点相对位置,不需要严格对应地理地貌。

5. 复现过程中踩过的四个典型坑

5.1 支路功率负号没处理,碳势出现负值

这是我最开始犯的错误。当时拿到潮流结果,直接遍历每个节点的下游支路做递推,把所有PF值当成正向功率用。结果出来一看,好几个节点的碳势是负数,支路碳流量也出现了“负碳流”,从物理上完全解释不通。

后来定位到原因:有一条支路在算例里的from/to方向和实际潮流方向相反,PF是负值。我没做方向翻转,等于让碳流沿着错误方向跑,本来该接收上游碳流的节点反而被当成“上游”输出碳流,自然算出负值。解决方式就是上面代码里的那一步:PF为负时交换from/to。这是整个复现里最简单也最致命的细节,强烈建议一进数据处理就先把方向统一。

5.2 交流潮流的损耗让碳流不平衡

我第一次跑的是交流潮流,最开始用的是rundcpf,后来为了“看起来更接近实际”改成runpf,结果全网碳流总和对不上:所有负荷碳足迹加起来,比发电机总碳量少了大约百分之三到五——刚好对应网络损耗的比例。机理就是损耗的碳量没地方去。

我当时用“损耗本地化”的办法处理,也就是把每条支路两端有功差(损耗)乘上该支路起点碳势,单独计入一个“损耗碳”汇总变量,并在校验时把它加进负荷侧。这么做物理上说得通,但后来对比了一下,很多论文根本不做交流潮流,直接用直流潮流结果,整个损耗问题就不存在。所以我的最终建议是:如果你只是为了复现碳流算法本身,用直流潮流;如果你想深究损耗碳流分配,再做交流潮流加损耗本地化的扩展。

5.3 参考平衡节点对结果的影响

IEEE 14节点的平衡节点是节点1,它承担的不只是自己的发电计划,还包括全网不平衡功率的调整。在无损耗情况下,平衡节点出力等于总负荷减去其他机组出力;在有损耗情况下,还得加上损耗。

复现时容易忽略的是:潮流结果里平衡节点的出力是求解出来的,不预先设定。如果平衡节点是煤电,系统的边际碳强度就在相当程度上由它决定。我在情景分析时试过把节点1的碳强度设成0,结果由于平衡节点出力很大,整个系统的碳势分布发生剧烈变化,很多节点的碳势直接跌到0.2以下。这提醒我们:平衡节点的碳强度设定要格外慎重,不能只按机组类型拍脑袋,要看论文里怎么约定的。

5.4 单位混用导致量级离谱

碳流计算涉及的功率单位有两种:Matpower内部是标幺值(基准容量100MVA),发电机组的碳排放强度常用单位是kgCO2/kWh或tCO2/MWh。两者直接相乘,出来的碳流量单位完全是乱的。

我一开始在算负荷碳足迹时,直接拿标幺值功率乘以tCO2/MWh的碳强度,得到的结果小了100倍,反应是“这系统碳排放也太环保了吧”。后来把所有功率转成有名值MW,碳强度统一成tCO2/MWh,负荷碳足迹才落在一个合理区间。建议在代码一开始就定义一个baseMVA=100,所有功率计算用标幺值做内部运算,但最终输出前统一转成MW再乘以碳强度,避免中途单位混乱。

6. 这套方法能往哪里扩展

6.1 多区域碳流与跨区绿电追踪

IEEE 14节点只是起点。实际省级电网动辄上千节点,碳流计算的价值在这种大规模系统里更明显——因为跨区输电、省间交易都会导致碳责任跨区域转移。把碳流算法封装成函数,输入是潮流结果和机组碳强度向量,输出是节点碳势和支路碳流密度,这个接口天然适合嵌入更大规模系统。很多研究里提到的“绿电追踪”“绿证溯源”,本质上也是在碳流框架下做类似的功率溯源,只是把碳强度换成绿电占比。

6.2 碳流信息在用户侧碳核算中的应用

另一个很实际的扩展方向是把它和用户侧碳排放核算结合。比如一个大型工业园区从电网购电,它的电力间接排放到底按什么因子算?如果直接采用省级电网平均排放因子,可能和实际相差很大。而如果所在节点靠近清洁能源基地,按节点碳势算出来的排放量会明显低于平均值,这也为后续绿电交易、碳减排认证提供了技术依据。我最近在做的项目里,就是把碳流计算输出作为前端展示和后端核算的中间层,效果还不错。

最后再分享一个小经验:复现这类算法,最重要的不是半小时跑通代码,而是把每一步的物理含义搞清楚。碳势为什么是这个值、支路碳流为什么往这个方向跑、换了机组碳强度后哪些节点会跟着变——只有这些都能解释清楚,才算真的“完美复现”。如果你按这篇文章的思路把IEEE 14节点的碳流算例跑通了,建议再自己改一改机组碳强度,或者换个更大一点的算例试试,很多细节会在对比里看得更透彻。

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

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

立即咨询