1. 碳排放流到底在算什么一个容易误解的切入点做电力系统碳排放分析最容易踩的坑就是把“碳排放因子乘总发电量”当成完整结论。这种总量核算应付报告可以一旦想回答“某个城市、某个工业用户、某条联络线上到底承担了多少碳”立刻失效。我这次在IEEE 14节点系统上复现电力系统碳排放流Carbon Emission Flow计算核心目的就是把“总量碳排”拆成“逐节点逐支路的碳流分布”再用Matlab代码把整套流程跑通。这个方向最近在EI期刊里出现频率很高因为双碳目标下电网侧需要的不再是一笔总账而是能够物理溯源的电碳耦合关系。碳排放流的思想其实不复杂电力潮流在网络里怎么走碳就跟着怎么走。发电机发一兆瓦时电如果碳排放强度是0.8吨二氧化碳每兆瓦时那么这0.8吨碳就“附着”在这度电上顺着输电路径流向负荷。只要电力系统潮流的分布已知碳流分布就完全由源端的发电碳强度和网络拓扑决定。换句话说碳排放流不是独立于潮流场之外的另一种物理场而是叠加在潮流场上的伴随场。正因如此计算的输入不依赖碳数据去重新仿真电力系统只需要获得稳态潮流结果再按我下面要讲的方法做一次矩阵化推演。我见过不少刚接触这个方向的同学一开始会纠结“碳排放流”是不是需要重新建立热力学模型甚至想去做机组燃烧侧的详细机理仿真。实际完全不需要。在输电网尺度下我们关心的是电量归属而不是烟囱里的化学反应。只要把每个发电节点的碳排放强度当成“源强”再把网络当成“运载管道”剩下的就是线性代数问题。IEEE 14节点系统之所以成为标准验证算例就是因为它规模适中既有环网、多发电机、多负荷又足够小到可以逐条支路核对结果。用Matlab实现时既可以把每一步矩阵打印出来验证也能快速调试这是拿它做EI级复现的重要原因。2. IEEE 14节点算例的数据基石从标准参数到碳强度配置2.1 标准数据怎么读才不出错IEEE 14节点系统的公共数据表大家都能找到一般包含四张表母线bus参数、支路branch参数、发电机gen参数、以及负荷数据。但很多人第一次直接用Matpower时会被数据列顺序搞晕。比如支路表里第1列是起始节点第2列是终止节点第3、4列是电阻和电抗第5列是充电导纳第6列是变压器变比使用时不能只看数字一定要对照文档确认单位。14节点系统中变压器支路不止一条尤其在节点5-6、4-7、4-9等位置变比参数默认不为零直接影响潮流分布进而影响碳流分布。另外标准IEEE 14节点算例的基准功率是100 MVA所有电抗、电阻都是标幺值。Matlab复现时如果忘了基准值换算后续碳流结果会整体偏移。我在做复现时习惯先把所有数据读成结构体再统一乘上基准功率得到有名值。例如支路有功潮流从标幺值转换为MW需要乘baseMVA而发电机出力同样要转换。碳流密度单位是tCO2/MWh这个单位天然是“和功率相乘再相加”的不会因为基准值的改变而变化但潮流功率必须和碳排放强度处在同一时间尺度上。这一点我会在避坑章节专门展开。2.2 机组碳排放强度该怎么设置EI论文里做IEEE 14节点复现时通常不会直接给出发电机碳强度而是自己设置一套场景。最常用的做法是把节点1、2上接的常规机组设为燃煤机组碳排放强度取0.9 tCO2/MWh左右节点3、6、8上的机组设置为天然气联合循环或零碳机组碳强度取0.4~0.6 tCO2/MWh如果做绿色电力场景也可以把水电、光伏设为0。问题在于碳强度设成多少直接决定节点碳势绝对值但不会改变碳流的相对分布。因此复现对标时重点看的应该是节点碳势排序、支路碳流密度分布规律以及系统总碳排是否守恒而不是纠结绝对数值必须和某篇论文完全一致。我在实际配置中还考虑到了一个容易忽略的细节发电机接入节点可能有多个机组但IEEE 14节点标准数据里每个发电机节点只给了一台机的出力上限和下限。碳强度赋值按节点维度处理就足够了不必细分到机组台数。因为碳排放流算法的输入是“节点注入功率乘碳强度”同一节点上若有多个机组其等值碳强度是各机组出力加权平均。如果你要复现的论文里特别标明了每台机组的碳强度那么就需要在构建发电机矩阵时按行对应而不是按节点对应。这一点做错算出来的节点碳势会直接失真。2.3 数据矩阵化从表格到计算把数据读进来之后我建议先把原始表格转成三个核心矩阵节点注入功率向量P_inject发电出力减负荷单位MW、发电机碳强度向量E_gen维度与发电机数量相同、以及支路连接表branch_table包含起始节点、终止节点和潮流结果。这样设计是因为后续碳流计算主要围绕“哪条支路流入哪个节点”来判断而不是直接操作整个导纳矩阵。我写的Matlab代码里数据准备部分大概长这样% 读取IEEE14节点标准数据需要提前准备好bus、branch、gen表 baseMVA 100; bus bus14; % n x 131列节点编号3列有功负荷 branch branch14; % l x 9前两列为起止节点 gen gen14; % m x 211列节点编号2列有功出力 % 构建节点发电向量 P_gen zeros(size(bus,1), 1); E_gen zeros(size(bus,1), 1); % 按节点存储发电机碳强度 for i 1:size(gen,1) idx gen(i,1); P_gen(idx) gen(i,2) * baseMVA; % 标幺值转MW E_gen(idx) carbon_factor(idx); % 自定义碳强度 tCO2/MWh end % 构建节点负荷向量 P_load bus(:,3) * baseMVA; % 支路表后续用于方向判断 branch_info [branch(:,1), branch(:,2)];这里有两个小提醒。第一是节点编号必须是连续整数因为后续要用节点编号作为数组索引如果数据里有0或者跳号所有矩阵索引都会错位。第二是如果直接用Matpower的loadcase函数一定要先mpc loadcase(case14);再取mpc.bus、mpc.branch、mpc.gen因为这些结构的列顺序和裸表格稍微不同不熟悉的人容易少取一列。3. 碳排放流计算核心算法矩阵化推导与步骤拆解3.1 节点碳势方程从“流入权重”看懂碳在节点上如何混合碳排放流的核心变量是节点碳势物理含义是“该节点流出的每单位电量里携带了多少碳”。我习惯把它理解成“电的含碳浓度”。发电节点含碳浓度等于发电机运行时的碳排放强度而中间节点没有发电含碳浓度等于所有上游支路流入电量的加权平均碳浓度。如果有多个电源从不同路径汇聚到同一个节点电能会完全混合所以节点碳势只有一个值。用公式表达对任意节点i如果其功率注入方向是正的即该节点作为汇那么节点碳势满足节点i碳势 (该节点所有上行支路碳流率之和 该节点发电碳排率) / (该节点所有上行支路流入功率之和 该节点发电出力)这里的“上行支路”需要根据潮流实际方向判断。比如支路从j流向i那么这条支路就属于节点i的上行支路其碳流密度等于节点j的碳势支路碳流率等于密度乘支路有功潮流。如果节点i同时有负荷负荷只是将电从节点抽出不改变节点碳势因为负荷不向节点注入能量。因为网络存在环网和多个电源这个方程是隐式的节点碳势依赖于其他节点碳势其他节点碳势又可能反过来依赖它。所以在实现时不能用简单的顺序遍历必须要么构建线性方程组求解要么用迭代法逼近。IEEE 14节点规模小用迭代法足够稳定我用的是Gauss-Seidel式的更新方式实测迭代不到10次就能让碳势变化小于1e-6。3.2 支路碳流密度与负荷碳流率的计算得到所有节点碳势后支路层面的计算就非常简单。在稳态潮流下如果支路有功潮流方向是从节点i流向节点j那么这条支路的碳流密度就是节点i的碳势。碳流率等于碳流密度乘以支路有功功率。这里不需要再对每条支路单独建模因为电能是均匀混合的“从高碳节点流出的电天然带高碳属性”。负荷碳流率则是节点碳势乘以该节点的有功负荷功率。系统总碳排放量自然地分解到每个负荷节点上这就是“碳排放流”最直接的价值给每个用户算‘碳排放责任’。当然输电线路上的网损也会带走一部分碳这部分碳既没有到达负荷也没有到达对端母线而是耗散在线路上。如果追求完整守恒需要单独计算网损碳流率。对于IEEE 14节点这样规模较小的系统网损约占总发电功率的3%左右对应的碳排占比更小。但在顶刊复现中网损碳流通常会被单独列出一项以保证“总发电碳排 所有负荷碳流之和 所有线路网损碳流之和”。碳流计算的最核心矩阵方程我整理成了下面的伪代码% 输入P_flow - 支路有功潮流方向已按从起点到终点标准化 % branch_from, branch_to - 支路首末端节点编号 % P_inject - 节点净注入功率发电-负荷 % P_gen, E_gen - 发电功率和碳强度向量 % 输出e_node - 节点碳势rho_branch - 支路碳流密度C_load - 负荷碳流率 n length(P_inject); e_node zeros(n, 1); e_node_new zeros(n, 1); % 判断每条支路的实际流向得到每条支路上的“源节点”索引 source_idx zeros(length(P_flow), 1); for k 1:length(P_flow) if P_flow(k) 0 source_idx(k) branch_from(k); else source_idx(k) branch_to(k); P_flow(k) abs(P_flow(k)); end end % 迭代计算节点碳势 max_iter 100; for iter 1:max_iter for i 1:n % 找到流入节点i的所有支路 inflow_power 0; inflow_carbon 0; for k 1:length(source_idx) if branch_to(k) i || branch_from(k) i % 这里需要判断该支路是否流入i根据source_idx和另一侧判断 % 若source_idx(k) ~ i 且支路另一端是i则说明功率流入i end end % 加上本地发电机 inflow_power inflow_power P_gen(i); inflow_carbon inflow_carbon P_gen(i) * E_gen(i); if inflow_power 1e-9 e_node_new(i) inflow_carbon / inflow_power; else e_node_new(i) 0; end end if norm(e_node_new - e_node, inf) 1e-9 break; end e_node e_node_new; end这段伪代码只是把框架列出来真正可运行的代码还要处理支路流入判别的索引逻辑不能简单比较branch_to i || branch_from i因为这个条件同时包含流入和流出。正确写法是对于支路k取该支路两端的节点id判断其中一端是否是节点i且另一端不是源节点。实现时可以先建立一个数组每条支路有一个“源节点”那么当支路的另一端非源端点等于节点i时这条支路才流入节点i。边界条件处理好整个迭代就没问题。3.3 为什么用迭代法而不是直接解矩阵理论上节点碳势可以写成形如e B \ c的线性方程组但我实践下来觉得迭代法更可控原因有两个。第一自动适应潮流方向变化。支路有功方向可能在迭代过程中发生小波动如果直接构造矩阵方向变化后矩阵需要重建而迭代法每轮更新时就按当前方向计算天然鲁棒。第二容易加入“注入功率接近0”的特殊节点。这些节点可能会出现分母为零的情况迭代法里单独设置阈值即可而矩阵法需要额外处理奇异。对于IEEE 14节点迭代法完全够用而且代码逻辑清晰方便和论文里的手算过程对照。4. Matlab代码实现从潮流结果到碳流分布4.1 潮流计算自写牛拉法还是用Matpower碳排放流计算的前置条件是稳态潮流结果。如果你是教学或自用建议直接用Matpower的runpf函数它内部实现了牛拉法对IEEE标准算例开箱即用。我这次的代码就是先调用runpf(mpc)得到result.branch(:,14)有功潮流、result.bus(:,3)负荷功率和result.gen(:,2)发电机出力。用Matpower的好处是能自动处理变压器变比、无功、机端电压控制等问题避免自己写潮流时因为无功不收敛而卡住。如果你希望完全脱离Matpower自己实现牛拉法也不难但需要注意IEEE 14节点系统中PQ节点、PV节点的编号约定以及平衡节点的处理。牛拉法的雅可比矩阵在14节点系统里是36阶左右不计平衡节点代码量不大但调试成本比较高。我建议第一次做碳排放流复现时直接把重心放在碳流计算本身潮流就用成熟工具。等整套流程跑通再回头替换/实现自己的潮流算法也不迟。4.2 碳流计算主函数一份可直接复用的代码骨架下面这段是我整理后可以复用的碳流计算主函数骨架输入输出都做了注释。为了兼容不同Matlab版本我建议所有变量命名用英文代码里不要写中文注释避免编码问题后面避坑章节会细说function [e_node, rho_branch, C_load, C_total] calc_carbon_flow(P_flow, branch_from, branch_to, P_load, P_gen, E_gen) % CALC_CARBON_FLOW Calculate carbon emission flow for given power flow. % Inputs: % P_flow - branch active power flow (MW), lx1 % branch_from - from-bus vector, lx1 % branch_to - to-bus vector, lx1 % P_load - nodal load vector, nx1 (MW) % P_gen - nodal generation vector, nx1 (MW) % E_gen - nodal generation carbon intensity, nx1 (tCO2/MWh) % Outputs: % e_node - nodal carbon potential (tCO2/MWh) % rho_branch - branch carbon flow density (tCO2/MWh) % C_load - load carbon flow rate (tCO2/h) % C_total - total generation carbon emission rate (tCO2/h) n length(P_load); l length(P_flow); e_node zeros(n,1); % Standardize branch direction: P_flow positive means from - to source zeros(l,1); P abs(P_flow); for k 1:l if P_flow(k) 0 source(k) branch_from(k); else source(k) branch_to(k); end end % Build inflow branch list for each node: branch index, source node, power for iter 1:100 e_new zeros(n,1); for i 1:n inflow_power 0; inflow_carbon 0; for k 1:l % branch k is inflow to i if its non-source end equals i other branch_from(k); if other source(k) other branch_to(k); end if other i inflow_power inflow_power P(k); inflow_carbon inflow_carbon P(k) * e_node(source(k)); end end inflow_power inflow_power P_gen(i); inflow_carbon inflow_carbon P_gen(i) * E_gen(i); if inflow_power 1e-6 e_new(i) inflow_carbon / inflow_power; else e_new(i) 0; end end if max(abs(e_new - e_node)) 1e-9 e_node e_new; break; end e_node e_new; end % Compute branch carbon flow density and rate rho_branch zeros(l,1); for k 1:l rho_branch(k) e_node(source(k)); end C_branch rho_branch .* P; % tCO2/h % Compute load carbon flow rate C_load e_node .* P_load; % tCO2/h % Total generation carbon emission C_total sum(P_gen .* E_gen); end这里有几个值得注意的实现细节。第一个是支路流入判定的逻辑other变量的写法确保了对每条支路都能找到“不是源节点”的那一端。如果该端是节点i则功率从源节点流向节点i确实形成流入。第二个是迭代收敛判据我用无穷范数阈值取1e-9实际运算一两秒就完事。第三个是P_flow方向标准化这一步必须在支路潮流矩阵进入函数前完成或在函数内部做abs处理否则碳流密度符号会错。4.3 结果可视化让碳流分布看得见算完碳流之后纯数字表格很难发现问题。我习惯用两种图来验证节点碳势柱状图和支路碳流分布图。节点碳势柱状图可以直接用bar(e_node)然后把节点号标上去。因为IEEE 14节点系统发电机集中在1、2、3、6、8如果碳强度配置合理这五个节点的碳势会呈现明显差异而下游负荷节点的碳势应该介于相邻电源碳势之间。支路碳流分布图可以用quiver或者plot配合支路坐标绘制但更简单的方法是画成带颜色的线线宽正比于碳流率。在这个图上你能直观看到碳从高碳机组所在节点向周围扩散高碳支路线特别粗低碳区域则很细。如果某些支路碳流方向和预期相反多半是潮流方向没有正确标准化。先看图再查表能省大量调试时间。5. 与EI论文复现结果对标节点碳势排序与误差讨论5.1 到底要对标哪些指标很多同学做“EI完美复现”时以为把论文里的节点碳势图跑得一模一样就算成功。我实际对比后认为需要重点对标的指标是三个节点碳势排序、支路碳流密度相对大小、系统总碳排守恒关系。前两个反映方法实现是否正确第三个反映是否有漏算或重复计算。拿我设定的场景举例节点1和2的燃煤机组碳强度0.9节点3的燃气机组0.5节点6和8的零碳机组0。潮流计算后节点1、2附近的碳势明显高于节点6、8附近。如果论文里也是类似的配置那么节点碳势的排序应该基本一致高碳区域在节点1、2所在的上游低碳区域在节点6、8及其下游。不过要注意如果论文同时接了火电、水电并明确设定了不同出力排序关系可能会因为潮流方向不同而改变不能机械照搬。5.2 常见误差来源分析误差主要来自四个方面。第一是碳排放强度单位不统一。有的论文用gCO2/kWh你的Matlab代码里用的是tCO2/MWh相差1000倍如果忘了换算节点碳势会整体错一个量级。第二是支路潮流方向定义。Matpower的result.branch(:,14)是有正负号的正号表示从from端流向to端负号相反。碳流计算函数里必须基于这个符号构造source向量。第三是网损处理。如果不把网损碳流单独拿出来用“总碳排除以总负荷”来反推平均碳势结果会和节点碳势加权平均值差一个网损系数。第四是负荷功率单位。从mpc.bus里取的负荷列是标幺值需要乘baseMVA否则碳流率会偏小100倍。5.3 如何确认自己的实现是“完美复现”我个人的复现流程是这样先不看论文结果自己把代码跑通然后打印系统总碳排、节点碳势总和、负荷碳流率总和先做守恒校验。如果sum(C_load)加上网损碳流与C_total的相对误差在1%以内说明计算逻辑基本正确。然后再对照论文给的曲线重点比较趋势。如果趋势一致但绝对值有偏移大概率是碳强度设置不同调整后就能完全对齐。如果趋势都不一致优先检查潮流方向标准化的代码因为这是碳流计算最容易埋雷的地方。守恒校验其实是一条非常有用的经验IEEE 14节点系统里无论网络怎么连接发电碳排总量必须等于所有负荷碳流率之和加上所有网损碳流率之和。如果这个等式不成立说明不是漏算了一条支路就是在节点碳势迭代里把流入/流出方向搞错了。我在第一次跑通代码时就发现迭代逻辑里把流出支路当成了流入结果节点碳势全部偏低总碳排比负荷碳流率大不少。当时怎么查都没发现问题最后就是靠守恒校验定位到方向判定的bug。所以新手在实现前先把这个守恒关系写在注释里调试时会少走很多弯路。6. 实操避坑与经验总结这些细节决定计算结果是否可信6.1 数据格式和单位换算的常见坑Matlab里读取IEEE14节点标准数据最常见的坑有三个。第一个是节点编号从1开始如果原始数据里编号从0开始必须在读入后加1。第二个是branch表中的变压器支路有“初始相角差”反映在变比列中这和Pu标准有关。第三个是发电机表里的Qmax、Qmin列可能为0不影响有功潮流但如果你是直接用gen(:,2)取有功出力千万不要把列索引搞错。另外我在处理碳强度向量时一开始用的是carbon_factor zeros(n,1)然后在发电机循环里赋值。这样没问题但要注意E_gen向量的维度是按节点数还是按发电机台数。如果后续代码写成E_gen(gen(:,1))其实Matlab允许向量按索引赋值但这个写法容易掩盖发电机重复接入同一节点的情况。最稳妥的做法是显式循环每个发电机节点单独赋值并检查是否有重复索引。6.2 潮流不收敛与支路功率过小的问题使用Matpower时runpf通常直接收敛但如果你修改了负荷或出力可能触发不收敛。这时候不要急着改算法建议先检查负荷是否超过该节点所在区域最大传输能力。IEEE14节点系统虽然小但如果把某个节点负荷调到300 MW以上潮流计算会失败因为标准算例的初始状态没有那么大传输裕度。这个问题在做碳流灵敏度分析时尤其容易出现。还有一个容易被忽略的情况某些支路在轻载时有功潮流绝对值可能非常小比如小于0.001 MW。此时abs标准化后的方向可能受潮流计算数值误差影响一会儿是正一会儿是负。我处理的办法是设置一个功率阈值P_min 1e-4低于阈值的支路直接忽略碳流贡献或者保持上一轮迭代方向。对IEEE14节点这样的小系统影响很小但对实际大系统很有必要。6.3 Matpower版本差异和Matlab编码问题Matpower从6.0到7.xrunpf的返回结构基本稳定但个别版本中branch列含义有扩展比如增加了故障前、故障后数据。如果直接按固定的列索引取值建议先从result.branch表头确认或者用mpc.branch的字段名取数。我见过有人因为用了Matpower 7.0后把原来的14列索引当成15列导致潮流结果完全不对。在碳流计算主函数里最安全的做法是从result结构中手动指定字段而不是依赖列数。Matlab的编码问题也值得单独说。近两年用Matlab 2023b、2025等版本时很多人写代码喜欢加中文注释结果换到另一台机器就乱码。原因多半是中文编码格式不统一尤其从网页复制代码时常常自带乱码。如果你像我一样需要分享或复现建议所有代码文件统一使用英文注释或者把中文注释放到Markdown说明文档中。这也是我在这套代码里刻意不用中文注释的原因——代码能直接被别人无痛运行比多几行中文注释重要得多。6.4 从IEEE 14节点扩展到更大系统的方法如果你后续要跑IEEE 39节点、118节点甚至实际省级电网碳流计算主函数不用改只需要保证输入数据格式统一。需要注意的是迭代方法的收敛速度会随着系统规模增大而变慢这时可以改成线性方程组直接求解。另外大系统中支路方向可能因环网而复杂建议在碳流计算前用拓扑搜索建立所有节点的上行支路集不要每次都遍历全部支路否则几百条支路计算量虽不大但调试麻烦。对于实际电网系统碳排放强度不再是一个固定值而是随着机组出力变化而变化的动态参数这时需要把机组碳排放强度建模成出力或煤耗的二次函数在每次潮流计算后更新一次碳强度向量再进入碳流计算。这个方向其实比“完美复现IEEE14节点”更有工程价值但底层逻辑和我在本文梳理的完全一致。把标准算例吃透后再做扩展效率会高很多。我个人实际操作下来最大的感受是碳排放流计算方法本身并没有门槛真正的门槛是把所有细节规范化尤其是潮流方向、单位换算、数据列索引、守恒校验这四个环节。只要这四个环节不出错IEEE 14节点算例基本一遍跑通。整套Matlab代码运行时长不到半秒剩下的时间都花在校验和理解结果上。如果你正在被这个方向困扰建议先把我前面讲的方向标准化逻辑写进代码再跑一遍系统总碳排守恒校验很多问题会立刻清晰起来。