最近在做一个电力系统碳排放流的Matlab复现项目目标算例是标准的IEEE 14节点系统。论文里的碳排放流公式看起来并不复杂——无非是把潮流结果拿过来按比例共享原则追一遍功率来源——但真正落到代码上节点编号错一位、支路首末端取反、网损不处理结果就完全没法看。这篇博文把整个复现过程写透从算例数据准备、核心矩阵推导到代码模块划分和避坑点适合正在跟碳排放流论文、想做电碳耦合核算或者刚接触IEEE节点算例的同行参考。1. 碳排放流到底在算什么把发电侧的“碳责任”摊到每一个用电节点1.1 为什么需要碳排放流从“发电机核算”到“用电侧追溯”传统的电力碳排放核算基本是“谁发电谁负责”燃煤机组烧了多少煤按燃料含碳量折算成碳排放记到发电机头上。这样算简单直接但有个问题——电在网上是混着走的用户从插座取电时根本不知道取到的电来自哪台机组。一台燃煤机组和一个风电场同时并网负荷侧实际消纳的是两者混合后的电不能简单地把所有排放都算给煤电。碳排放流解决的就是这个问题。它基于物理潮流的分布结果把发电机侧的碳排放沿着支路“分摊”到每个节点、每个负荷身上。算完之后可以回答节点5的负荷每用一度电对应的碳排放是多少克某条支路上传输的电量中有多少碳排放责任是从燃煤机组来的。这就是“碳势”的概念——相当于给每一度电贴了一个“碳排放浓度标签”。这个需求在绿电交易、产品碳足迹核查、园区级碳排放监测场景下非常实用。你不可能在每个用户门口装一个二氧化碳检测仪但只要你掌握了系统潮流和发电侧排放强度就能通过碳流计算把排放责任拆解到每个负荷节点精度还很高。1.2 比例共享原则流量混合后按比例追溯是最合理的假设碳排放流的核心理论基础是比例共享原则。想象两根水管汇成一根一根是热水一根是冷水汇合之后你取一瓢水你没法精确说出这瓢里到底有多少热水、多少冷水——你唯一的合理假设就是按两股水的流量比例分配。电力潮流也一样多个支路功率流入一个节点后在节点母线处完成电气混合之后从该节点流出的任何功率只能认为是从上游各支路按流入比例带过来的。比例共享原则具体表述是从某节点流出的功率中来自某条注入支路的比例等于该注入支路功率占该节点总注入功率的比例。这不是精确的物理定律而是电力系统碳流分析公认的近似假设。它看上去简单却是整个碳排放流矩阵方程成立的地基。这个原则也意味着碳排放流计算不能改变原有潮流结果——物理潮流是多少就是多少碳流只是在潮流之上做标定和追溯。所以整个计算链路里第一步永远是先跑一个准确的潮流解后面的碳流计算全部依赖这个结果。1.3 碳排放流与传统潮流计算的关系不改变物理潮流只改变责任追溯很多人第一次接触碳排放流会误以为要重新解一次方程组。实际上完全不是这样。传统的牛顿-拉夫逊潮流求解的是电压幅值和相角得到节点注入功率和支路潮流碳排放流是利用潮流结果做线性追溯计算。只要系统拓扑和运行点确定潮流结果定了碳流结果也就唯一确定了。这也带来了一个好处碳流计算可以和现有电力系统分析流程解耦。你想研究某个风电场接入后对负荷侧碳排放的影响只需要重新跑一次潮流再把发电机组碳排放强度更新为新的机组数据碳流计算模块不需要改动。我在代码实现时也特意把碳流求解封装成独立函数输入是潮流结果和发电碳排放强度向量输出是节点碳势、支路碳流率等结果这样后续扩展场景时非常方便。2. IEEE 14节点算例准备Matpower潮流数据与碳流计算的前处理细节2.1 IEEE 14节点系统的结构与关键节点IEEE 14节点系统是最常用的电力系统测试算例之一结构不大但足够验证算法正确性14个节点、20条支路含3台变压器、11个负荷系统里有5台机组分布在节点1、2、3、6、8。其中节点3、6、8上的机组在标准算例里往往作为同步调相机运行也就是有功出力为0只提供无功支撑。这块要特别提醒同步调相机节点虽然不发出有功但它是电压支撑节点在潮流计算中不能省略。在碳排放流计算中调相机有功出力为0对应的发电碳排放强度乘以0后还是0所以它们不参与碳流注入但节点本身作为拓扑节点仍然要保留在矩阵里。标准IEEE 14节点数据可以从Matpower中直接获取Matpower用case14函数就能生成该算例的完整数据。需要注意区分Matpower内置的14节点数据5台机组版本和某些论文里修改过的14节点数据——我看到过一些论文把3、6、8节点换成实际发电机组碳流结果自然会不同。2.2 用Matpower加载与求解标准算例使用Matpower时先加载算例数据然后调用潮流函数求解基态潮流mpc case14; % 加载IEEE 14节点标准算例 result runpf(mpc); % 牛顿-拉夫逊法求解潮流runpf返回的结构体result里包含我们需要的所有数据bus矩阵节点数据、branch矩阵支路数据、gen矩阵机组数据。result.bus里每一行对应一个节点关键列包括节点编号第1列、有功负荷Pd第3列result.branch里每一行对应一条支路包括首端节点第1列、末端节点第2列、支路有功潮流PF第14列Matpower中单位是MWresult.gen里是有功出力Pg第2列。这一步务必做一次数据合理性检查所有节点有功负荷加起来是否等于总负荷所有发电机有功出力减去总负荷再减去总网损是否约等于0。如果潮流不收敛后面全是白算。2.3 从潮流结果中提取碳流输入向量碳流算法需要的输入可以归纳为四个向量/矩阵节点编号集合bus(:,1)一般是从1到14的连续编号但不要假设连续用unique提取更稳妥。节点有功注入向量发电机有功出力含调相机出力0减去该节点负荷。这里要注意如果有多个机组挂在同一节点各自的Pg要合计。支路首末端节点编号数组branch(:,1)和branch(:,2)用于构建拓扑连接关系。支路有功潮流向量branch(:,14)注意Matpower在支路数据里给出的PF是从首端流向末端的功率如果潮流反向PF会是负值。把这些数据从Matpower结构体里抽出来整理成碳流计算专用格式我建议转成列向量并保存成*.mat文件这样后面反复调试时不用每轮都重跑潮流。2.4 节点编号与数据类型的一致性检查这部分是整个前处理最容易出错的地方值得单独说。第一Matpower的bus矩阵本身是按节点编号排序的但你的后续矩阵构造如果直接使用循环基于节点编号定位要确保数组索引和节点编号一致。我习惯的做法是建立一个节点编号到矩阵行号的映射表例如node_index containers.Map(bus(:,1), 1:14)后续通过node_index(node)取行号。第二支路的首末端方向要和潮流方向区分开。首末端是物理拓扑方向PF的正负代表实际潮流方向。在构建“谁向谁送电”的流入矩阵时一定要用PF的实际符号判断而不是拓扑首末端。第三数据类型的坑也不小Matpower里所有功率单位都是MW但很多论文里的碳势和碳流率用的是MWh口径下的tCO2/MWh。单位不统一时最后算出来的支路碳流率可能差出三个数量级。建议从一开始就固定使用MW和tCO2/h这套单位组合发电碳排放强度的单位就用tCO2/MWh。3. 核心公式推导比例共享原则下节点碳势与支路碳流率的矩阵方程3.1 三个核心变量的定义节点碳势、支路碳流密度、支路碳流率碳排放流计算最终输出的核心变量有三个节点碳势e_n节点n上单位有功功率对应的碳排放速率单位tCO2/MWh。直观理解就是节点n的“电的碳浓度”。支路碳流密度ρ_ij支路ij传输的功率中单位电量带有的碳排放单位也是tCO2/MWh。按照比例共享原则支路碳流密度等于该支路起点节点的碳势即ρ_ij e_i。支路碳流率R_ij支路ij单位时间传输的碳排放总量单位tCO2/h计算公式为R_ij ρ_ij × P_ij其中P_ij是支路有功潮流。就这三个变量。整个碳排放流计算本质上就是把发电机的碳排放注入量按潮流分布“稀释”或“混合”到全网的节点和支路上最后算出一个能够表征“每度电碳含量”的空间分布场。3.2 构造节点有功通量谁在向这个节点送电要算节点碳势第一件事是统计每个节点的总流入有功功率。这里要用“流入”口径包括从与该节点相连的支路流入的功率以及该节点本地发电机的注入功率。用数学式表示对于节点nP_in(n) sum(P_ij for all j - n, where P_ij 0) P_G(n)其中P_G(n)是节点n上发电机的总有功出力。如果某条支路的潮流方向是从节点n流出即P_ij 0说明它属于流出侧不应该计入P_in(n)。在我的代码实现里遍历每条支路判断潮流正负号把正向潮流加到首端节点还是末端节点取决于PF正负号对应的方向。这一点和拓扑首末端解耦独立判断可以有效避免方向错误。3.3 矩阵方程与迭代求解两种落地形式在比例共享原则下节点碳势满足如下平衡方程e_n ( sum(R_ij for all j - n) E_G(n) ) / P_in(n)其中E_G(n) e_G(n) × P_G(n)是节点n本地发电机碳排放注入速率e_G(n)是该发电机的碳排放强度。把所有节点的方程联立起来可以写成线性方程组。设A为节点间“碳流转移矩阵”A(n, m)表示从节点m流入节点n的功率占节点n总流入功率的比例再乘上节点m的碳势再加上本地注入项最终得到(I - A) * e b其中e是节点碳势列向量b是归一化后的本地碳排放注入速率向量b(n) E_G(n) / P_in(n)。求解这个线性方程组即可得到所有节点的碳势。矩阵A的维度是14×14规模不大直接用MATLAB的\运算求解数值稳定性很好。另一种实现方式是迭代求解。把上面的式子改写成迭代格式e_new(n) ( sum(P_j_n * e_old(j) for all j - n) E_G(n) ) / P_in(n)初始值e_old全设为0然后反复迭代更新直到两次迭代所有节点碳势的最大绝对差小于1e-10。这个迭代法本质上就是高斯-赛德尔迭代的变体对于无环或少环网络收敛很快通常几十步内就能达到很高的精度。我个人推荐用矩阵法做主计算迭代法做交叉验证。两者结果一致基本可以确保矩阵构造没有方向性错误。3.4 负荷碳流率与系统总碳排放的计算口径有了节点碳势e_n节点n上的负荷碳流率可以直接算出来R_load(n) e_n × P_load(n)其中P_load(n)是节点n的有功负荷MW。把全系统所有负荷节点的碳流率加起来理论上应该等于全部发电机的碳排放注入速率之和。这就是系统级碳平衡校验公式sum(R_load(n)) ≈ sum(E_G(n))注意这个等式存在一个微妙的偏差来源如果全网网损没有做处理支路损耗会导致末端负荷承担的碳流率总和略小于发电侧注入总量因为损耗相当于“凭空消失了”一部分功率和对应的碳。这就是网损处理问题我在后面专门展开。4. Matlab实现链路如何把矩阵公式落成可调试的代码模块4.1 主程序框架数据层、计算层、结果层分离我复现时的代码结构是三个模块数据准备、碳流核心计算、结果输出与校验。这个结构的好处是把容易出错的数据整理逻辑和核心算法隔离开调试时能快速定位问题出在哪一层。% 数据层 mpc case14; result runpf(mpc); [node_data, branch_data] extract_data_from_mpc(result); % 碳流计算层 carbon_result carbon_flow_solver(node_data, branch_data); % 结果层 % 输出表格、校验、绘图数据层返回的node_data至少应包含节点编号、节点有功负荷、节点发电出力、节点编号映射branch_data包含首端节点编号、末端节点编号、支路有功潮流带符号。不建议直接在核心计算函数里读result.bus的原始列因为列索引太多极易写错列号。4.2 构建流入功率矩阵与线性方程组的Matlab实现下面这段是碳流核心计算里最关键的矩阵构造代码我简化了注释但仍可直接运行function carbon carbon_flow_solver(node_data, branch_data) n_bus length(node_data.bus_id); node_index flat_map(node_data.bus_id, node_data.bus_id); % 替代containers.Map % 节点总流入功率与发电碳注速率初始化 P_in zeros(n_bus,1); E_G zeros(n_bus,1); % 发电注入 for k 1:length(node_data.gen_bus) n node_index(node_data.gen_bus(k)); P_in(n) P_in(n) node_data.gen_pg(k); E_G(n) E_G(n) node_data.gen_emission(k) * node_data.gen_pg(k); end % 构建转移矩阵 A 和右端向量 b A zeros(n_bus, n_bus); b zeros(n_bus, 1); for k 1:size(branch_data,1) from_node branch_data.from_bus(k); to_node branch_data.to_bus(k); pf branch_data.pf(k); if pf 0 src node_index(from_node); dst node_index(to_node); else src node_index(to_node); dst node_index(from_node); end % src流向dst功率pf进入dst节点的总流入 P_in(dst) P_in(dst) abs(pf); A(dst, src) A(dst, src) abs(pf); end % 归一化转换先A按P_in归一化需避免除零 for n 1:n_bus if P_in(n) 1e-9 A(n,:) A(n,:) / P_in(n); b(n) E_G(n) / P_in(n); end end % 求解线性方程组 (I - A) * e b e (eye(n_bus) - A) \ b; % 支路碳流密度与碳流率 R_branch zeros(size(branch_data,1),1); for k 1:size(branch_data,1) from_node branch_data.from_bus(k); to_node branch_data.to_bus(k); pf branch_data.pf(k); if pf 0 R_branch(k) e(node_index(from_node)) * pf; else R_branch(k) e(node_index(to_node)) * abs(pf); end end carbon.e e; carbon.R_branch R_branch; carbon.R_load e .* node_data.pd; end这段代码里有一个容易忽略的细节P_in在计入支路潮流之后同时也参与了b的归一化。如果某个节点既没有发电注入也没有支路流入P_in为0那么碳势方程无定义实际系统里这种节点通常是负荷很小且上游断开的孤立节点。处理方式是在归一化时把该节点的碳势直接置0或置为0附近并在结果里标记出来。4.3 迭代求解的初值与收敛设置矩阵法求解基本上一步到位但我仍然建议实现一个迭代法做交叉验证。迭代代码的核心循环如下e_old zeros(n_bus, 1); max_iter 1000; tol 1e-10; for iter 1:max_iter e_new zeros(n_bus, 1); for n 1:n_bus inflow_carbon A(n,:) * e_old; e_new(n) (inflow_carbon * P_in_save(n) E_G(n)) / P_in_save(n); end if max(abs(e_new - e_old)) tol break; end e_old e_new; end这里有个数值经验迭代初值不能全给0吗可以但收敛会稍慢。更好的做法是把初值设置为所有发电机碳排放强度的加权平均值基本十几步就能收敛。另外迭代法和直接求解的差异如果超过1e-6几乎可以肯定是转移矩阵构造有方向性错误这个可以当作一个自检指标。4.4 结果整理表格、误差校验与可视化算完之后把结果整理成可读的表格。我通常输出三张表节点碳势表、支路碳流率表、负荷碳流责任表。表格用MATLAB的table类型列名用中英文都行关键是保留原始节点编号方便对比。可视化方面一个简单但有效的画法是把14节点画成拓扑图节点颜色深浅代表碳势大小支路粗细代表碳流率大小。这样能直观看出碳流从高碳机组所在的节点向低碳区域扩散的趋势。可视化代码不复杂核心是plot加text标注节点号这部分就不展开完整代码了。5. 结果自检与复现中的三个高频错误碳势越界、碳流失衡、网损处理5.1 基准检验统一碳强度下节点碳势必须完全相等做任何算例之前先做一个最简单的基准测试把所有发电机的碳排放强度设为同一个值比如1 tCO2/MWh。此时无论网络拓扑多复杂、负荷分布多不均匀系统内任意节点的碳势都必须等于1。如果算法正确节点碳势应该是完全相同的常量如果某个节点碳势偏高或偏低说明转移矩阵构造或方向判断有误。这个基准测试的优点是不需要手算任何数值就能判断算法对错。我建议把这段测试直接写成一个test_carbon_solver.m脚本每次修改代码后优先跑一遍。第一次实现时我就靠这个测试抓到了一个方向性错误当时把PF正负方向处理反了导致某些节点碳势超过1某些低于1基准测试立刻报警。5.2 碳守恒校验负荷碳流总和应对应发电碳排放总和第二个自检是系统级碳守恒。将全部负荷碳流率相加与全部发电碳排放注入率相加做对比。忽略网损的理想情况下两者应该严格相等。我的实践是计算两者的相对误差error abs(sum(R_load) - sum(E_G)) / sum(E_G)在实际IEEE 14节点算例中如果不处理网损这个误差通常在1%到3%之间具体取决于系统总网损率。如果误差超过5%基本可以判断不是网损问题而是某个负荷节点的碳流漏算或者矩阵行错误。5.3 网损处理的两种约定与它们对结果的影响网损是碳排放流复现里最大的“暗礁”。网络支路存在电阻损耗潮流从首端流入的功率大于末端流出的功率这中间“消失”的一部分功率对应着一部分碳排放。如果不处理负荷侧碳流总和会小于发电侧碳排放总量。我见过论文里常见两种处理方式。第一种忽略损耗即直接用首端潮流功率作为整条支路的传输功率。这样处理简单但系统碳平衡会产生偏差。适合网损率很低、对精度要求不高的分析场景。第二种等效负荷法。把每条支路的损耗功率首端PF减去末端PF注意都是按实际方向计算看作一个挂在支路末端的额外有功负荷参与碳流分摊。这种处理方式能严格保证碳流守恒也是工程复现中推荐的做法。实现时在提取支路数据阶段就额外计算每一条支路的损耗loss_k abs(pf_first) - abs(pf_last); % 按实际方向pf_first0时为正然后把所有损耗功率归入对应末端节点的负荷里。注意如果某条支路潮流方向为负损耗仍然是正值直接把abs(pf_first) - abs(pf_last)取正数就行。5.4 调相机与零有功节点对矩阵奇异性的影响第二个高频错误出在调相机节点和零注入节点。节点3、6、8上的调相机有功出力为0发电碳排放注入为0这些节点的P_in完全来自其他支路流入。如果某条连接它们的支路因为潮流方向判断错误导致P_in计算为0转移矩阵A的某一行就全是0线性方程组可能出现奇异或得到错误解。调试时如果\运算提示矩阵接近奇异或者某个节点碳势出现负值优先检查这些零有功注入节点附近的支路方向。还有一个经验不要把调相机节点从算例里删掉虽然它们不注入有功但作为拓扑连接点删掉会改变潮流分布碳流结果和原始论文就无法对应了。6. 从单算例到实际扩展复现论文级代码时我的实操体会6.1 复现前先统一符号约定不同论文里的符号差异很大有的用e表示支路碳流率有的用e表示节点碳势。我复现时第一步不是写代码而是先把目标论文里的所有变量写在一张纸上统一成下文这套记号节点碳势e_n、支路碳流密度ρ、支路碳流率R、发电碳排放强度e_G。然后把论文中的公式逐个翻译成这套记号再做代码实现。这一步看着麻烦其实能省大量返工时间。我曾见过有人按论文原符号写了一半代码发现某处公式里的矩阵维度对不上回头一查是原论文把碳流密度和碳流率混用了这种情况如果你没有自己的统一记号排查难度极大。6.2 别轻信“完美复现”每篇论文对损耗的处理口径不同标题里的“完美复现”只是宣传话术。实际复现论文代码时你会发现不同论文对网损、调相机、多机组节点并网的处理口径完全不同。某篇EI论文可能在算例里有意简化了损耗处理你按另一篇论文的方法复现结果自然对不上。我的建议是先做基准测试确定自己的算法内核正确再针对目标论文调整外围口径。如果你的结果和目标论文数值不完全一致大概率不是算法错误而是损耗处理、发电机碳排放强度取值或基准容量不一致。复现时要把这些口径差异记录在注释里方便后续查证。6.3 用3节点手推案例验证矩阵构造在直接上手IEEE 14节点之前我强烈建议先构造一个3节点小手算案例两个发电机节点一个负荷节点两条支路负荷处功率刚好是两个发电机各供一半。手算这3个节点的碳势和支路碳流率然后把代码里的算例数据替换成这个小案例跑一遍看是否和手算一致。这个验证步骤能帮你隔离两类问题一是你矩阵公式本身的错误二是从Matpower数据提取过程的错误。如果3节点案例都跑不对就不要浪费时间调试14节点了。这个小技巧是我过去复现各类碳流论文时最有效的提速手段。6.4 代码组织建议把碳流求解放进独立函数把碳流求解封装成carbon_flow_solver.m这样的独立函数输入输出都定义清楚好处是后续换算例、换场景时不用改核心代码。我在原项目里同时做了IEEE 14节点和IEEE 30节点两个算例主程序只改了一行mpc case14改mpc case30其他都不用动靠的就是函数封装。另外输入参数的顺序建议固定为(node_data, branch_data)别搞成(branch_data, node_data)这种小规范能避免同事或者未来的你自己在调用时犯低级错误。6.5 扩展方向连续时段碳势追踪与多区域分摊最后说下这套代码的扩展空间。碳排放流不只是算单一时段实际应用往往需要连续追踪一天24小时或全年8760小时的碳势变化。这种场景下把潮流计算封装进循环每个时段跑一次潮流再算一次碳流配合发电计划数据就能输出节点碳势时序曲线。如果再加上机组启停和发电碳排放强度的时序变化可以进一步分析电网碳排放的时空分布特征这也是当前电碳耦合核算中比较受欢迎的方向之一。我自己在复现过程中最大的体会是碳排放流的数学工具并不深奥真正决定代码能不能用的是对潮流数据细节的把控和对各种“口径陷阱”的敏感度。如果你也在做类似的复现工作建议从3节点手算案例开始把基准测试脚本写好再碰IEEE 14节点就会顺畅很多。
