做电力系统低碳化研究的朋友，十有八九都绕不开碳排放流这个工具。我在读文献时第一次看到“碳排放流”四个字，以为又是某篇论文里玄乎的新概念，直到自己动手在IEEE 14节点系统上把计算流程完整走了一遍，才发现它其实就是一套把“发电侧碳排放”按电网实际潮流“追责”到每个负荷节点的核算方法，原理不复杂，但真要写成Matlab代码跑出可靠结果，还是有相当多细节容易踩坑。
这篇博文就把我这套已在IEEE 14节点系统上完整验证过的碳排放流计算方法拆开来讲：从理论公式、数据准备、Matlab实现到结果分析一条线串下来。代码框架可以直接改造成你自己的算例，适合正在做电力系统低碳规划、碳追踪、碳责任分摊相关课题的研究生，也适合刚接触碳流计算、想把理论落到仿真上的工程师。
1. 碳排放流到底在算什么
1.1 电网是一池混合的“水”，碳是里面的“色”
理解碳排放流，最形象的方式是把它想成一根水管网络，发电机组从不同的水源往管网里注入带有不同颜色的水，负荷从管网末端取水。每条管道里流出的水是什么颜色，取决于它上游所有水源的混合比例——但你没法说清某一滴水的颜色具体来自哪个水源，只能按比例估算。
电网里的“颜色”就是碳排放强度，单位是 kgCO2/MWh（或者 tCO2/MWh）。煤电厂的碳强度很高，可能到0.9左右；燃气机组低一些，可能在0.4上下；水电、风电、光伏基本是0。电网是一个电气上连通的网络，功率从发电机流向负荷，天然会把不同机组的电能混在一起，因此在任何一个节点上取电，都要按该节点流入功率的“构成比例”来分摊碳排放责任，这就是碳排放流的核心逻辑。
这套思路在文献中称为比例分担原则（Proportional Sharing Principle），也是目前碳排放流计算最主流的基础假设。
1.2 三个核心指标：节点碳势、支路碳流密度、碳流率
碳排放流理论里有三个指标，几乎所有研究都绕不开它们。
第一个是节点碳势，可以理解为“在这个节点上取1MWh电所对应的碳排放量”，单位是 tCO2/MWh 或 kgCO2/MWh。节点碳势是一个状态量，只和该节点所有注入功率的碳强度加权平均有关，和负荷大小无关。
第二个是支路碳流密度，概念和电流密度类似，指的是单位有功功率流过支路时所携带的碳流量。在忽略网损或按送端节点处理时，支路碳流密度等于该支路送端节点的碳势。
第三个是碳流率，单位是 tCO2/h，表示单位时间内流经支路或被负荷消耗的碳流量。支路碳流率等于支路有功功率乘以该支路的碳流密度，负荷碳流率等于负荷功率乘以对应节点的碳势。
如果用一句话概括：节点碳势回答“这里用电有多‘脏’”，碳流率回答“这里用电产生了多少碳排放”，一虚一实，构成了碳流计算的基本输出。
1.3 为什么选IEEE 14节点系统做演示
IEEE 14节点是电力系统分析里的经典算例，规模适中，数据公开，全网只有14个节点、20条支路、5台发电/调相机组、11个负荷点。相比3节点等简单系统，它有足够的拓扑复杂度，能够体现环网中碳流按潮流方向分配的过程；相比IEEE 39节点或118节点，它又不至于让数据准备和调试图形过分繁琐。
更重要的是，14节点系统里既有常规发电机组（节点1、2、3），又有不带净有功出力的同步调相机（节点6、8），碳流计算时必须把这种情况单独处理，否则容易出错。把这个系统跑通了，后面换更大系统只是数据规模问题，算法逻辑不用变。
2. 系统数据准备：把“原料”备齐再开火
2.1 IEEE 14节点系统的拓扑构成说明
做任何仿真，先搞清楚系统长什么样。
IEEE 14节点系统的基准容量取100 MVA。三个常规发电机组分别接在节点1、2、3上，其中节点1是平衡节点，节点2和3是PV节点；节点6和8接有同步调相机，它们向系统提供无功支撑，但净有功出力视为0。负载主要分布在节点2、3、4、5、6、9、10、11、12、13、14，其中节点4、5、9等是较重的负荷节点。
支路部分包括变压器支路和输电线路支路，尤其要注意节点5到6之间、节点4到7等位置包含变压器，支路参数里存在非标准变比，潮流计算时会直接影响支路有功流向，而碳流计算完全依赖有功潮流结果，所以支路数据不能填错。
一般来说，IEEE 14节点原始数据里已经给出了所有节点、支路、发电机的完整参数，直接从公开数据源抓下来整理成Matlab可读的矩阵即可。
2.2 数据清单与存储格式
我在复现时按照三个矩阵来组织数据，尽量贴近Matlab的索引习惯：
bus = [ 1 1 0 0 0 0 100 1 0 0 0 230 1 1 2 2 21.7 12.7 0 0 100 1 0 0 0 230 1 1 % ... ];每一行是一个节点的编号、类型、有功负荷、无功负荷、并联电导电纳等。发电机数据单独放一个矩阵，每一行包含所在节点编号、有功出力、无功出力、电压幅值设定值。支路数据每一行是首端节点、末端节点、电阻、电抗、对地电纳、变比等。
这里有一个很重要的习惯：发电机数据矩阵里的“节点编号”不要直接当数组索引用，而是单独建一个索引向量，把发电机所在节点映射到发电机序号。否则后面的碳流矩阵拼接会乱套。
2.3 碳排放强度参数的设定逻辑
计算机组注入碳流率前，必须先给每台发电机组设定碳排放强度。实际工程中，碳强度可以来自实测、机组类型缺省值或碳配额数据；学术演示时一般按机组类型直接给典型值。
本算例我按以下方式设定：
| 节点 | 机组类型 | 有功出力/MW | 碳排放强度/(tCO2/MWh) |
|---|---|---|---|
| 1 | 燃煤机组 | 约232.4 | 0.9 |
| 2 | 燃气机组 | 40.0 | 0.4 |
| 3 | 燃气机组 | 40.0 | 0.45 |
| 6 | 同步调相机 | 0 | 0 |
| 8 | 同步调相机 | 0 | 0 |
注意潮流计算所得节点1平衡机出力会随负荷水平变化，碳排放强度取0.9这一典型煤电值。调相机没有净有功出力，其碳排放强度设为0不影响结果。这里也可以把节点2和节点3的强度设成不同的值，从而更明显地在结果里体现“不同机组上网电量混在一起之后的碳势差异”。
3. Matlab代码实现：从潮流结果到碳流矩阵
3.1 程序整体架构
我建议把整个计算分成三个模块：潮流计算模块、碳流计算主模块、结果可视化模块。这样换系统时只改数据文件和柱子即可，不需要动核心算法。
- 潮流计算模块：读取网络数据，用牛顿-拉夫逊法或者调用MATPOWER的runpf得到全系统精确潮流解。本次复现我以MATPOWER作为潮流求解器，因为它处理14节点这种小系统非常稳定，还能直接输出支路有功矩阵。
- 碳流计算主模块：输入潮流结果、发电机有功、碳排放强度，计算节点碳势、支路碳流密度和碳流率。
- 结果可视化模块：把碳势和碳流率用图形方式展示出来。
3.2 支路有功潮流的提取与矩阵化
碳流计算最关键的输入是全网各支路的有功潮流量值和方向。潮流计算后，支路有功功率是一个 from-to 二维矩阵，其中 F(i,j) 表示从节点 i 流向节点 j 的有功功率。提取出来后要按节点序号建立完整的 n×n 矩阵。
% MATPOWER的branch结果: [from, to, pf, qt, ...] % F矩阵初始化 F = zeros(nb, nb); for k = 1:size(branch, 1) from_bus = branch(k, 1); to_bus = branch(k, 2); pf_flow = branch(k, 14); % MATPOWER里pf列对应的支路有功 F(from_bus, to_bus) = pf_flow; F(to_bus, from_bus) = -pf_flow; % 反向为负 end这里要特别小心：负功率并不代表“倒流”，而是说明在当前潮流解下，实际有功功率方向与该支路参数定义的首末端方向相反。碳流计算必须以实际流向为准，后续构造流入矩阵时要取正值的部分。
3.3 节点碳势的线性方程组构造
根据比例分担原则，节点 k 的碳势 E_N(k) 等于流入该节点的总碳流率除以流入该节点的总有功功率，计算式如下：
[ E_N(k) = \frac{\sum_{i \in G_k} P_{G,i} E_{G,i} + \sum_{j\in IN(k)} F_{j,k} \cdot E_N(j)}{\sum_{i \in G_k} P_{G,i} + \sum_{j\in IN(k)} F_{j,k}} ]
其中 (E_{G,i}) 是节点 i 上发电机的碳排放强度，(F_{j,k}) 是从节点 j 流入节点 k 的有功功率。这个方程组是线性的，可以整理成矩阵形式：
[ (\mathbf{P}{in,diag} - \mathbf{F}{in}) \cdot \mathbf{E}N = \mathbf{R}{G,node} ]
其中 (\mathbf{P}{in,diag}) 是对角矩阵，对角线元素是各节点总流入有功；(\mathbf{F}{in}) 矩阵的 ((j,k)) 元素是从 j 流入 k 的支路有功（因此每个节点的流入项里减去对应支路碳流）；(\mathbf{R}_{G,node}) 是各节点上发电机的注入碳流率向量。
用Matlab求解这个线性方程组极其简单：
% P_in_total: n×1 节点总流入有功（发电机净出力+所有正流入支路功率） % F_in: n×n 流入矩阵，F_in(j,k)表示从j到k的支路有功（按实际方向，非负） % R_G_node: n×1 发电机碳流率，R_G_node(k) = P_G(k) * E_G(k) A = diag(P_in_total) - F_in; b = R_G_node; E_N = A \ b;这里需要说明一点：如果直接用 (A \backslash b) 求解，前提是矩阵 A 可逆。实际电网拓扑下该矩阵满秩成立，但如果你改了系统数据出现“矩阵奇异”警告，多半是某些节点没有连接到任何电源也没有任何支路注入，属于孤立节点，需要先检查数据。
3.4 支路碳流密度与碳流率的批量计算
得到节点碳势向量 E_N 后，支路碳流密度直接按送端节点碳势赋值，每条支路送端碳势对应其首端节点碳势：
rho_branch = zeros(nb, nb); for k = 1:size(branch, 1) from_bus = branch(k, 1); to_bus = branch(k, 2); % 实际潮流方向 if F(from_bus, to_bus) > 0 rho_branch(from_bus, to_bus) = E_N(from_bus); R_branch(from_bus, to_bus) = F(from_bus, to_bus) * E_N(from_bus); else rho_branch(to_bus, from_bus) = E_N(to_bus); R_branch(to_bus, from_bus) = F(to_bus, from_bus) * E_N(to_bus); end end节点满足基尔霍夫电流定律，碳势则用比例分担原则计算。全网碳流同样满足流守恒：发电机注入的碳流率之和，应当等于所有负荷消耗的碳流率加上网络损耗对应的碳流率。这句话是检验整个计算过程是否正确的黄金标准。
4. 结果可视化与分析：光有数据可不够
4.1 节点碳势分布：一张图看懂系统“谁更绿”
数值计算完成后，第一个动作是画节点碳势柱状图。下面这段代码生成各节点碳势的直方图：
figure; bar(E_N, 'FaceColor', [0.2 0.6 0.3]); xlabel('节点编号'); ylabel('节点碳势/(tCO2/MWh)'); title('IEEE 14节点系统节点碳势分布'); grid on;在样例参数下，节点1自身碳势接近0.9，节点2、3的碳势接近0.4~0.45，而负荷节点（4~14）的碳势会介于三者之间，具体数值取决于它们从哪些支路获得功率。环网支路越多，节点之间的碳势差异越小，因为功率混合得越充分。
工程上，节点碳势高的区域就是“高碳电”集中落地的区域，在这个区域新增负荷会带来更高的碳排放增量，适合以这个指标指导低碳调度和低碳规划。
4.2 支路碳流图：把碳“流”画出来
柱状图只是第一步，我建议把碳流叠加到系统拓扑图上，用线条粗细表示碳流率大小。这样做的好处是能直观看出碳流的主要通道和阻塞点。
实现思路是：先绘制IEEE 14节点的地理接线图（可以用原始坐标或手摆坐标），然后把支路碳流率映射到线宽，用颜色映射碳流密度，这样“哪条线在输碳、碳有多浓”一目了然。
% 简化示意：线宽与碳流率成正比 for k = 1:size(branch, 1) from_bus = branch(k, 1); to_bus = branch(k, 2); if abs(R_branch(from_bus, to_bus)) > 0 line_width = max(0.5, abs(R_branch(from_bus, to_bus)) / 2); plot([X(from_bus), X(to_bus)], [Y(from_bus), Y(to_bus)], ... 'LineWidth', line_width); end end画这样的图能帮助你发现很多矩阵数据看不出来的问题，比如某条支路的碳流率异常偏大，往往是潮流方向反转导致的逻辑bug。
4.3 负荷碳流率对比：算清楚“谁该为碳买单”
碳排放流研究的最终目的大多落在碳责任分摊上，所以我会额外输出一张负荷碳流率表：
% 节点k的负荷碳流率 R_load = P_load .* E_N; fprintf('节点%d 负荷功率 %.2f MW，碳流率 %.4f tCO2/h\n', k, P_load(k), R_load(k));从这张表里能非常直观地看到，某个负荷节点尽管用电量不大，但因为节点碳势高，其碳排放责任反而比另外一个用电量更大的节点更重。这在碳配额分配和绿电消费认证中是非常关键的信息。
4.4 参数敏感性：碳强度变化对碳势分布的影响
我这里还建议做一个简单扩展：把节点2的燃气机组碳强度从0.4改成0.2，再跑一遍，观察全网负荷节点碳势的降幅。这个操作本质上是在模拟“低碳机组替代高碳机组”的调度效果，也能帮你验证程序对输入参数的敏感性，确保算法逻辑没有把碳强度设进去却没有任何反馈。
如果调整碳强度后只有个别节点碳势变化，而其他节点几乎不动，说明流程里可能漏掉了该机组对应支路与负荷节点之间的拓扑关系，要回到流入矩阵的构造上排查。
5. 常见问题与排查技巧实录
5.1 矩阵维度对不上，索引错位
这是复现碳流计算时最常遇到的一类问题。很多人直接把发电机所在的节点编号当作数组下标来访问，但发电机节点是“1、2、3”，计算矩阵可能是按“1~14”全节点编号构建的，两者在该节点没有发电机的位置上就出现了错位。
我的建议是全部用“节点编号映射表”来处理：
node_idx = zeros(1, nb); % node_idx(bus_number) = 数组序号 for k = 1:nb node_idx(k) = k; end % 发电机循环时用 node_idx(gen_bus(k)) 替代 gen_bus(k)同时，所有矩阵的尺寸统一为 nb×nb，再小的功能也不要另造尺寸，减少错位概率。
5.2 节点碳势出现负值或不合理数值
负碳势基本可以断定是流入矩阵构造错了。最常见的错误是把支路负方向功率直接当作正向流入了，导致某些节点的注入功率被抵消，总流入甚至出现负值，方程组求解出来的碳势就会出现离谱数值。
排查方法很简单：把 F_in 矩阵打印出来，逐行核对每个节点的流入支路和潮流方向；尤其注意平衡节点，它的净注入有功很大，但流入矩阵里不能把“净出力”和“支路流入”重复叠加。
5.3 潮流计算不收敛或结果精度不够
MATPOWER对IEEE 14节点通常不会出问题，但如果自己写牛顿-拉夫逊潮流，不收敛的原因多半在初值或参数单位。IEEE标准数据里功率基准是100 MVA，有些原始数据文件已经转成了标幺值，有些还是有名值，混用时支路导纳矩阵会差好几个数量级。
另外，变压器支路的变比默认为1.0时要检查原始数据是否为“非标称变比”，如果漏掉变比数值，环网潮流方向可能完全反过来，碳流结果自然全错。
5.4 碳平衡校验不通过
校验方法如下：
total_generation_carbon = sum(R_G_node); total_load_carbon = sum(R_load); total_branch_loss_carbon = abs(sum(sum(R_branch))) ... % 实际用流入流出差 % 理想情况下: total_generation_carbon ≈ total_load_carbon + network_loss_carbon如果两边差距明显，优先检查是否有负荷节点被漏算，或者发电机净出力与总线注入功率之间是否存在局部消耗。IEEE 14节点的网损在几个MW量级，其对应碳流率一般在总碳流里的占比很小，但如果发现网损碳流异常高，反而说明碳势向量本身可能已有偏差。
5.5 常见问题速查表
| 问题现象 | 可能原因 | 排查思路 |
|---|---|---|
| 碳势矩阵求解时报奇异 | 孤立节点、流入矩阵缺失 | 检查F_in每行是否有非零注入 |
| 部分节点碳势为0 | 该节点无支路流入且无发电机 | 检查是否漏建支路 |
| 某节点碳势超过所有机组值 | 支路方向错误叠加 | 打印F_in，核对潮流方向符号 |
| 全网碳平衡偏差过大 | 负荷或发电机节点漏算 | 逐一比对潮流输入输出 |
| MATPOWER结果与手算不一致 | 数据单位、变比问题 | 校验潮流功率平衡再继续碳流计算 |
写在最后的一点体会
这个算例我从理论推导到代码跑通，前后折腾了两三天，最难的不是公式，反而是把“潮流结果”转成“碳流输入”时那些索引和数据对齐细节。只要把数据结构和电场校验逻辑做好，碳排放流计算本身的代码量并不多，核心求解就是解一个线性方程组，所有复杂度都在数据预处理和结果校验上。
如果你也是刚接触碳流计算，建议不要一上来就追求完整复现论文里的复杂场景，先把IEEE 14节点这套流程跑通，把碳平衡校验做到1%以内，再扩展到更大系统。后续有条件的话，还可以在现有框架上加入网损碳流分摊、储能充放电碳流分析甚至碳-电耦合市场结算，这套基础代码都能直接作为底层支撑扩展。
最后送一个小技巧：计算前先手动估计一下各节点碳势的合理区间，比如煤电节点接近0.9、燃气节点接近0.4、纯负荷节点不高于上游最高碳势，这样程序一跑出来，结果合不合理你心里马上有数。
