说实话，刚拿到这个任务的时候，我心里想的是：一篇EI论文的复现项目，电力系统碳排放流计算方法，IEEE 14节点，Matlab实现——这几个关键词凑在一起，难度应该不大。毕竟碳排放流本质还是基于潮流结果做碳分摊，潮流计算有现成工具，公式无非是矩阵运算。真正动手以后我才发现自己错得离谱。
论文里关于碳流密度的定义、网损的归算方式、比例分摊假设，每一条都暗藏选择。选错了，结果对不上论文，你甚至会怀疑是不是自己的代码写错了。这篇文章把整个复现过程完整拆开：从碳排放流原理、IEEE 14节点数据准备、Matlab实现细节，到最终结果验证和踩坑记录，给想做类似复现工作的同学一条可以直接参考的路径。
1. 复现前先理清：这篇论文到底要解决什么问题
1.1 为什么传统碳排放统计不够用
传统电力系统碳排放核算的方式，是在发电厂关口装表计量，用发电量乘以机组排放因子，得到全厂排放总量。这种方法只能给出一个“总盘子”，回答不了几个非常现实的问题：某条10千伏馈线送出去的电里到底包含了多少碳？工业园区的负荷用电，对应的碳责任应该算在哪个电厂头上？电网的运行方式改变之后，碳排放会怎样重新分布？
碳排放流方法的核心，就是把“碳”当成一种看不见的物质，让它附着在有功功率上，跟着潮流的路径从发电机流向负荷。这个思路和电力潮流的逻辑完全一致。电流从多个电厂汇入电网，在节点上混合，又从不同支路分别流向负荷，物理上我们区分不出哪些电子来自哪个电厂，但碳排放流通过一套分摊规则，可以把每一度电的碳足迹从源头一直追到末端负荷。
这套方法一出来，电网层面的碳排放核算就从“厂界总量”变成了“全链路可观测”，每个节点的碳势（单位电量所含碳量）都可以算出来，每条线路的碳流率也可以算出来。这是传统方法做不到的。
1.2 复现论文时最大的障碍不是公式，是选择
我拿到的这篇EI论文，结构非常典型：先定义碳排放流的基本概念，再给出无损网络下的计算方法，然后推广到有损网络，最后在一个标准算例上演示结果。看起来逻辑一环扣一环，实际动手写代码时，问题几乎全藏在“推广”两个字里。
无损网络下的公式很漂亮：每个节点的碳势，等于流入该节点的总碳流率除以流入该节点的总有功功率。可一旦网络有网损，公式里就要额外处理“损耗这部分能量对应的碳到底算给谁”。论文在这个地方往往只有两三句话，不会把细节展开，甚至会省略推导过程。你在实现时就必须自己做一个决定：用平均功率做近似？还是把网损等效成虚拟负荷？还是上迭代修正？
论文给出的其他参数也经常不完整。比如发电机碳排放强度，论文可能只给一句“参数见附录”，而附录里只有表格没有单位；机组的类型倒是写了“燃煤”“燃气”“水电”，但具体的排放因子需要你自己找参考值。这些细节每个单独看不难，组合在一起就容易让人反复推翻重来。
团队里经常有同学把复现理解成“翻译代码”，拿到论文就把公式一个个改成Matlab语句。这是最危险的思路。公式后面站着假设，假设背后还有物理意义。如果你不理解作者为什么这样定义碳势，你连错误都不知道出在哪一步。
2. 碳排放流的核心原理：节点碳势、分支碳流与比例分摊
2.1 三个关键量的定义
碳排放流里最核心的三个量，我在项目里天天跟它们打交道。
第一个是节点碳势，记作 rho，含义是单位电量承载的碳排放量，常用单位是 kgCO2/MWh。它反映的是“这个节点上所有流入电能混合后的平均碳浓度”。好比很多条水流汇进一个水池，每股水里的墨水浓度不一样，充分混合之后池水的浓度就是节点碳势。
第二个是支路碳流率，等于支路首端有功功率乘以首端节点碳势，单位是 kgCO2/h。它表示一小时内通过这条支路“输送”了多少碳。注意这里用的是首端功率，因为碳跟着有功走，首端碳流在传输过程中会因为网损有所衰减，但那是随后要处理的损耗分摊问题。
第三个是负荷碳流率，等于负荷有功功率乘以负荷所在节点的碳势。这个值特别有意义，它直接把用户的用电量换算成了对应的碳排放责任。工厂、数据中心、居民小区的月度电费单如果配上碳流单，看到的就是这个数字。
我建议做复现时先把单位统一好。论文里有的用 tCO2/MWh，有的用 gCO2/kWh，相差三个数量级。我自己全程统一用 kgCO2/MWh，读数据时顺手转换，避免后面核验时被数量级迷惑。
2.2 比例分摊原则是碳流追踪的基本规则
碳排放流不是物理上可测量的量，它是我们定义出来的，因此“怎么分摊”必须有合理规则。目前主流方法用的是比例分摊原则。
假设有一个节点，上游来了三股功率流，下游分两条支路送电，同时还有一个本地负荷。碳流法这样处理：节点碳势是流动混合的结果，所有下游支路和负荷，都按各自功率在总流出功率中的占比，分享这个节点的总碳流率。功率大，分到的碳就多；功率小，分到的碳就少。这就是比例分摊，简单说就是“按比例分碳”。
这个原则翻译成公式，节点碳势等于流入节点的总碳流率除以流入节点的总有功功率。注意分子是各上游支路碳流率和发电机碳注入的总和，分母是各上游支路有功和发电机有功的总和。它天然满足了比例分摊的假设，同时可以写成线性方程组，用矩阵一次性求解所有节点的碳势。
我用一个生活化的场景帮大家记忆：一条主管道里的水带着碳墨水，流到一个分叉口，分成了粗细两根管子。粗管流量大，带走的墨水自然多，细管流量小，带走的墨水就少。每根管子里的墨水浓度和主管道完全一样。这个“浓度一样，流量不同”就是比例分摊的直观图像。
2.3 有损网络的网损归碳方式决定结果精度
如果系统完全无损耗，上面的公式可以直接用，节点碳势定得死死的。但真实电网不是这样，IEEE 14节点系统里的线路有电阻，重载时线损可能占负荷的几个百分点。网损对应的那部分能量也有碳排放责任，总得算到某个对象头上。
实际实现里常见两种做法。第一种是等效负荷法：把每条支路的功率损耗看作支路末端的一个额外负荷，参与碳势分摊。这样做物理意义明确，实现起来也不难，只要在每一条有损支路末端增加一个负的注入功率项，让碳势方程自动把损耗对应的碳分摊给所有下游负荷。第二种是迭代修正法：先用无损模型算一次节点碳势，然后用支路首末两端功率的平均值代替首端功率，重新计算碳势，反复迭代几次。实测下来，迭代三次左右碳势变化就小于万分之一。
我复现时优先采用了迭代修正法。原因是它不需要改变网络拓扑结构，代码改动量小，而且和某篇论文描述的结果对齐时误差可控。如果论文里没说清楚网损处理方式，强烈建议你检查它结果表里的总数，否则你的结果和论文结果会有系统性的偏差，而且偏差不是在某个节点，是所有节点一起偏。
3. IEEE 14节点系统的数据准备与潮流计算
3.1 标准测试系统的数据构成
IEEE 14节点系统是电力系统领域使用频率最高的标准算例之一，规模不大但结构完整。全系统一共14条母线，20条交流支路，其中包含输电线路和变压器支路，5台发电机分别挂在1号、2号、3号、6号和8号母线，负荷分布在11个不同的母线上。就验证碳排放流算法来说，这个系统非常合适：规模不算大，手动核对节点方程不现实；但也足够复杂，能逼你把矩阵构建、方向处理、网损归算这些细节全部考虑到位。
需要特别提醒的是，网上流传的IEEE 14节点数据有多个版本。不同版本的变压器变比、支路电阻电抗、负荷有功无功，甚至基准电压都可能不一样。这直接导致潮流结果不同，碳流结果自然跟着变。如果你直接拿一个网上下载的case14数据去复现论文结果，大概率对不上。
我的处理办法是：先把论文提供的系统参数表完整梳理一遍，逐条录入Matlab数据结构。如果论文没给全，就选用公认的标准版本数据，并在最终文档里明确标注数据来源。这一步看起来繁琐，但能帮你避免后面所有“结果对不上”的后续排查。
3.2 用Matlab求潮流并输出碳流计算所需边界条件
碳流计算的前提，是有一个收敛的稳态潮流结果。这里我直接用了Matpower工具包里的标准14节点算例做底子，然后根据论文需要调整负荷和发电机出力。如果是自编潮流，也可以用牛顿-拉夫逊法，但完全没有必要重新造轮子，把时间省下来放在碳流模块上更划算。
跑潮流之前先确定几个边界条件：平衡节点是1号母线，它负责吸收全系统的功率不平衡量；发电机节点给定有功出力和机端电压幅值；负荷节点给定有功和无功。这些数据和标准算例里是一致的，调整参数时不要让发电机出力之和与负荷之和差得太大，否则潮流很难收敛。
潮流跑完以后，要提取的边界条件包括：每台发电机的有功出力、每个负荷节点的有功功率、每条支路首末两端的有功功率和线路功率损耗。这里有一个初学者特别容易犯的错：把发电机的额定容量当作注入功率送去算碳流。额定容量只是设备上限，实际出力才是碳流方程的输入。不同的运行方式下，同样一台发电机出力不同，它对全网碳势的贡献也完全不同。
4. Matlab代码实现：从碳流矩阵到结果输出
4.1 程序框架与模块划分
拿到算例数据和潮流结果后，我按照功能把整个工程拆成了六个模块，每个模块单独调试，最后串起来跑主程序。
main_carbon_flow_14.m % 主程序入口 ├── load_system_data.m % 读取母线、支路、发电机、负荷数据 ├── run_power_flow.m % 调用潮流计算并返回支路功率 ├── build_carbon_matrix.m % 构建碳势线性方程组 ├── solve_carbon_potential.m % 求解各节点碳势 ├── compute_carbon_flow.m % 计算支路碳流与负荷碳流 └── plot_results.m % 输出可视化图表模块化带来的最大好处是定位问题快。当结果不对时，我可以单独打印中间变量，检查矩阵A的结构、右端项b的数值，很快就能判断是潮流环节出错还是碳流方程构建出错。如果全部堆在一个脚本里，排错成本会高很多。
在build_carbon_matrix模块里，最重要的事情是建立支路功率的方向约定。我用的是“首端送末端”的正向假设，支路数据处理时先把潮流结果排序，保证branch_pf是首端送到末端的正向有功。方向一旦混了，碳势方程直接乱套，而且不容易察觉，因为解出来的数看起来“挺合理的”。
4.2 核心代码：构建碳势方程并求解
下面这段代码是碳流计算最核心的部分，我把关键逻辑简化后贴在下面。
% 构建碳势方程: A * rho = b % 物理含义: 流入节点的碳流量 + 发电机碳注入 = 流出节点的碳流量 nb = size(bus, 1); % 节点数 nl = size(branch, 1); % 支路数 A = sparse(nb, nb); b = zeros(nb, 1); % 1) 遍历每条支路，把"首端节点碳势 * 支路有功"作为流入项加到末端节点行 for k = 1:nl i = branch(k, 1); % 首端节点编号 j = branch(k, 2); % 末端节点编号 Pij = branch_pf(k); % 首端送末端的有功功率 if Pij > 0 A(j, i) = A(j, i) + Pij; else % 如果潮流反向，说明实际方向是 j -> i A(i, j) = A(i, j) - Pij; end end % 2) 每个节点行加上流出项: 负荷功率 + 该节点所有出线功率之和 for i = 1:nb P_load_i = load_p(i); % 节点负荷有功 P_out_i = sum(branch_pf(branch_from == i)); % 流出支路功率总和 A(i, i) = A(i, i) + P_load_i + P_out_i; % 发电机注入碳流量作为右端项 b(i) = P_gen(i) * rho_G(i); end % 3) 求解线性方程组 rho = A \ b;这段代码里，方程组的第i行逻辑是：所有上游支路送到节点i的碳流率之和，加上节点i本地发电机的碳注入，等于节点i流向下游各支路的碳流率之和，加上节点i本地负荷带走的碳流率。把所有已知的发电机碳注排放到右端项b，把包含未知节点碳势的项放在左边系数矩阵A中，正好就是标准形式A*rho=b。
求解用左除运算符即可。对于IEEE 14节点这种规模的系统，稀疏矩阵求解瞬间完成。如果要扩展到几百节点的系统，这个表达式依然很快。
对于进入迭代修正的情况，我在主循环里把支路功率更新为平均功率，重新组装A矩阵，循环三次后rho变化已经非常小。这样做的好处是不需要改变矩阵结构，只改动Pij为0.5*(Pij首位+Pij末端)。
4.3 结果输出与可视化
碳势求解完成之后，计算结果包括三个部分：节点碳势向量rho、支路碳流率、负荷碳流率。为了快速检查，我用柱状图把节点碳势画出来，看高低排序是否符合直觉——靠近大燃煤机组的节点碳势应该偏高，水电和负荷中心附近的节点碳势应该偏低。
figure; subplot(2,1,1); bar(rho); xlabel('节点编号'); ylabel('节点碳势 (kgCO2/MWh)'); title('IEEE 14节点系统节点碳势分布'); subplot(2,1,2); bar(load_carbon); xlabel('负荷节点编号'); ylabel('负荷碳流率 (kgCO2/h)'); title('各负荷节点碳流率');如果要复现论文里的网络拓扑彩图，可以把节点坐标整理好，用plot函数画网络图，节点颜色映射碳势高低、支路粗细映射碳流率大小。这一步是结果展示的点睛之笔，也是论文图表最容易对照复核的部分。我建议先出数值表，再画图，避免图表好看但数字经不起推敲。
5. 复现结果怎么验证：守恒校验与算例对比
5.1 第一道校验：全系统碳守恒
碳流计算完成后，第一件事不是画图，而是做全系统的碳平衡校验。
整个系统里，碳排放的来源只有发电机，所以发电侧总碳排等于所有发电机有功出力乘以各自碳排放强度再求和。碳的流向只有三个出口：负荷碳流、网损碳流、线路储能变化（稳态下为零）。因此必须满足：
发电机总碳排放 = 所有负荷碳流率之和 + 所有支路损耗对应的碳流量之和
这条等式如果误差超过0.1%，说明碳流方程构建有误。最常见的错误是支路方向反了、网损没有归到任何对象、负荷节点漏算。我把这个校验做成主程序里的一道断言，不通过就直接报错，避免拿错误数据去做后面的分析。
在IEEE 14节点系统里，线路损耗对应的碳流虽然占比不大，但完全不计会让生成侧总排和负荷侧总排对不上。如果只算到负荷碳流为止，发电机总排会比负荷碳流大，多出来的量正好是网损部分。这个差值本身就是一个快速检验指标。
5.2 第二道校验：最小系统测试
在跑完整系统之前，我建议先构造一个极简单的两节点系统做单元测试：一台燃煤发电机直接带一个纯有功负荷，线路电阻为零。
这个系统里没有任何拓扑复杂性，理想情况就是负荷节点碳势等于发电机碳势，负荷碳流率等于负荷功率乘以发电机碳排放强度。如果这个用例在代码里算出来不对，那问题基本出在矩阵构建逻辑本身，而不是IEEE 14节点系统的细节问题。
我复现时先用的是三节点辐射状系统，一台发电机送两路负荷，这样可以验证比例分摊是否正确。比如某台发电机碳势400kg/MWh，送出总功率100MW，两个负荷各50MW，那么两个负荷的碳流率应该各是20000kg/h，且两者之和等于发电机总碳排40000kg/h。这种手工可算的用例，比任何调试器都管用。
5.3 第三道校验：与论文结果图表对比
基础校验通过之后，才轮到和论文结果比对。比对分几层：数值表、柱状图、网络拓扑彩图。
论文里如果给了各节点碳势表格，你把两个向量放在一个表格里对比，记录相对偏差。如果论文给的是负荷碳流率饼图，需要根据占比推算数值。还有一种情况是论文只给了一张节点碳势分布图，没有给数字，那就要靠目视对比量级和大小排序。
需要特别说明：如果你的网损处理方式和论文不一致，数值上会有几个百分点的偏差，这属于正常现象。重点要看的是节点的排序趋势和整体量级是否一致。如果论文里某个节点碳势很高，你的结果里却偏低，那大概率是发电机碳势参数设置不同，或者该节点的上游支路方向判断错了。我自己遇到过一次，后来发现是把一条变压器支路的首末端写反了，导致下游整个区域的碳势全错了。
6. 踩坑记录与扩展建议
6.1 数据版本不一致导致的坑
IEEE 14节点系统的公开数据有多个版本，不同版本之间最明显的差异是变压器支路的参数和负荷功率水平。某些版本的1号母线和2号母线之间的变压器变比略有不同，导致潮流分配变化，碳势结果整体偏移。
我的做法是优先使用论文作者公开的数据文件，如果没有，就固定使用一个权威数据源，并做版本备注。跑复现项目最忌讳的是左右换数据，换一次数据结果全变，最后连你自己都不知道哪个版本对应论文的哪张图表。
6.2 平衡节点碳势的初始值问题
在碳势方程里，发电机的碳注入是已知右端项，平衡节点的处理相对特殊。平衡节点的有功出力不是预先固定的，它由潮流计算结果决定，所以在碳流计算中如果把它当作普通节点一样代入，可能出现方程数不够或者冗余的情况。
我的处理方式是把平衡节点的碳势方程单独处理成给定值，或者干脆从方程组里移除平衡节点的行，最后再用已经解出的碳势去算平衡节点的碳流率。这种处理不会影响其余节点的碳势结果，因为平衡节点的碳势是由其连接关系决定的，它的碳注入量在发电侧总排里照常计入。
6.3 有损网络的迭代收敛问题
用迭代修正法处理网损时，迭代初值用无损模型的节点碳势，然后用平均功率重新组装矩阵。这里有个小技巧：平均功率不是简单的算术平均，而是应该对每条支路分别取首端和末端功率的平均。重载支路的损失占比大，迭代修正的效果更加明显。
我在IEEE 14节点系统上实测，迭代三次之后节点碳势的相对变化就低于1e-4。如果你发现迭代一直不收敛或者振荡，检查一下是不是把支路损耗当成负荷注入到了错误的位置，这种错误会让方程失去对角占优性，数值上表现为震荡或负碳势。
6.4 发电机碳势参数怎么定才合理
论文如果只给机组类型，不给具体排放因子，你可以参考行业常用缺省值：燃煤机组大约800到1000 kgCO2/MWh，燃气机组大约400到500 kgCO2/MWh，水电机组和风光机组按零碳处理。这些数值不是绝对的，论文算例的重点是验证方法，不是给出精确的碳排放核算报告，所以只要你把参数设定写清楚，结果就是可复现的。
我还做过一个灵敏度实验：把燃煤机组排放因子从850调到950，系统负荷节点碳势几乎线性跟着变，但节点之间的相对排序完全不变。这说明方法本身的稳健性很好，参数的绝对取值影响的是量级，不影响趋势。复现时如果和目标论文差值大，先看排放因子，再看网损处理方式。
6.5 后续可以怎么扩展
这套代码做完之后，扩展空间很大。把14节点改成IEEE 30节点或者118节点系统工程，只需要替换数据文件和调整节点编号映射，碳流计算模块基本不用动，矩阵求解依然秒级完成。我后来就把这套逻辑移植到了一个省级电网模拟数据上，效果很稳定。
更进一步可以做时序碳排放流分析，把一天24小时的负荷曲线放进去，得到节点碳势的日变化曲线，能看出峰荷时段哪些区域碳势升高，这对碳电价的制定有直接的参考价值。也可以做碳流灵敏度分析，研究某条线路扩容后，不同区域负荷碳势的变化幅度，用来指导电网低碳改造的投资优先级。
最后分享一个做这类复现项目最深的体会：复现EI论文的重点从来不是把代码跑通，而是把论文作者在公式背后的每一个选择都还原出来。碳排放流这个方向，公式看起来大同小异，但实现细节里的一点点偏差，就会造成结果的系统性偏移。建议你拿到论文后先准备一个3节点手算算例，把每一步都用手算核对一遍，再上IEEE 14节点系统，你会感谢这个习惯的。
