去年帮几个研究生复现一篇关于电力系统碳排放流的EI论文时我发现绝大多数人卡住的根本不在碳流算法本身而是上游潮流结果的方向符号、节点功率归集这些“边角料”没处理好导致节点碳势矩阵解出来一团乱麻。这篇文章就是把我自己从数据准备、数学推导到MATLAB代码实现的全过程整理出来针对IEEE 14节点标准系统把电力系统碳排放流的计算方法拆开讲明白代码可以直接套用。文章会覆盖碳排放流的基本概念、节点碳势方程的推导、IEEE 14节点系统的数据准备、MATLAB完整实现以及计算结果验证目的是让刚接触这个方向的研究生和工程师少走弯路。我这里用的思路和代码都是经过实际验证的对应的是EI期刊里那种“基于潮流追踪的碳流计算”主流框架算出来的结果拿去和论文对比基本能对齐。1. 为什么碳排放流能火起来它到底解决了什么问题1.1 传统碳排放核算的尴尬先说说我自己的感受。之前做区域电网碳排放核算最常用的办法就是“发电侧排放因子法”——把所有电厂发电量乘上各自的排放因子然后按总量除以总电量得到一个全网平均排放因子。这个方法在宏观层面没毛病但一旦涉及用户侧责任分摊就非常尴尬。举个例子同一座城市里A工厂所在片区由附近的风电场和燃气电厂联合供电B工厂所在的片区由远距离火电送电。按传统的“全网平均排放因子”来算A和B的用电碳排放一模一样。但物理上A确实在用更多的绿电。这在碳交易、绿电认证、产品碳足迹这些实际场景里是站不住脚的。碳排放流Carbon Emission Flow就是在这样的背景下被提出来的。它的核心思想是碳排放随着有功功率在网络中流动功率走到哪里对应的“碳责任”就跟到哪里。这样就能把发电侧的碳排放责任按潮流路径精确分摊到每一个负荷节点、每一条支路、每一个用户。1.2 碳排放流的基本逻辑碳排放流的思路其实不复杂它和电力潮流的逻辑很相似。潮流告诉我们有功功率从哪个节点流向哪个节点、流过多少碳排放流则在这个基础上给每一股功率“染色”——煤电、气电、水电各自带有不同的碳排放强度它们在电网里混合、分流最终到达每个负荷节点时可以算出这个节点的综合碳排放因子。这个概念在英文文献里叫Carbon Emission Flow国内学者清华康重庆老师团队等做了大量系统化工作EI期刊上相关文章很多。核心术语有三个节点碳势某个节点上单位电量对应的碳排放量单位是g/kWh或tCO2/MWh支路碳流密度某条支路上单位电量携带的碳排放量支路碳流率单位时间内流过某条支路的碳排放总量单位是t/h或kg/h。这三个量有严格的递推关系支路碳流密度等于该支路功率流出端节点的碳势支路碳流率等于支路有功功率乘以支路碳流密度。只要算出每个节点的碳势后面全部都能算。1.3 为什么选择IEEE 14节点很多入门者会问为什么EI论文里一抓一大把都是IEEE 14节点原因很简单规模适中14个节点、20条支路、5台发电机既能反映复杂的环网结构又不至于让潮流和碳流计算变得难以调试数据是公开的标准算例每个人手里拿到的是同一套数据算法结果具有可比性支路中包含变压器支路和多电压等级能充分验证碳流方法在复杂拓扑下的适用性。换句话说IEEE 14节点就是碳流算法的一个“标定台”。把14节点的实现吃透了换到IEEE 30节点、118节点甚至实际省级电网只是数据规模的变化核心逻辑完全一样。2. 碳排放流的核心机理节点碳势与碳流方程2.1 三个基础概念的严格定义先定义一套贯穿全文的符号体系。系统有n个节点节点编号1到n。经过潮流计算后我们可以得到每个节点的发电机注入有功功率PG,i每个节点的有功负荷PL,i每条支路的有功潮流Pij方向以潮流实际方向为准。在此基础上定义节点碳势ei节点i的单位用电量所对应的碳排放量g/kWh。它代表了“在该节点取用1kWh电能所承担的平均碳排放责任”。支路碳流密度ρij支路ij单位电量携带的碳排放量g/kWh。由于碳排放是跟随功率流走的有$$ \rho_{ij} e_{\text{upstream}} $$其中e_upstream是功率流出发端节点的碳势。比如功率从节点i流向节点j则ρij ei。支路碳流率Rij单位时间流过支路ij的碳排放量等于支路有功功率乘以支路碳流密度$$ R_{ij} P_{ij} \cdot \rho_{ij} $$单位是g/h或t/h。2.2 节点碳势方程的推导节点碳势的计算是整个方法的中枢。思路是进入某个节点的总碳排放量等于该节点所有注入功率来源携带的碳排放之和而节点的碳势就是这个总碳排放除以总注入功率。以节点i为例。设所有通过支路向节点i注入功率的节点集合为S_in,i节点i上发电机组集合为G_i。节点i的总注入功率为$$ P_{\text{in},i} \sum_{j\in S_{\text{in},i}} P_{ji} \sum_{g\in G_i} P_{G,g} $$其中Pji是从节点j流向节点i的有功功率PG,g是节点i上第g台发电机的出力。流入节点i的总碳排放为$$ C_{\text{in},i} \sum_{j\in S_{\text{in},i}} P_{ji} \cdot e_j \sum_{g\in G_i} P_{G,g} \cdot e_{G,g} $$其中ej是上游节点j的碳势eG,g是第g台发电机的碳排放强度。于是节点i的碳势为$$ e_i \frac{C_{\text{in},i}}{P_{\text{in},i}} \frac{\sum_{j\in S_{\text{in},i}} P_{ji} e_j \sum_{g\in G_i} P_{G,g} e_{G,g}}{\sum_{j\in S_{\text{in},i}} P_{ji} \sum_{g\in G_i} P_{G,g}} $$注意这个方程中ei在等式左边而等式右边的ej是其他节点的碳势。把所有节点的方程联立就得到一个以e为未知量的线性方程组。由于电网拓扑是连通的方程组有唯一解。2.3 矩阵化求解思路将上面n个节点的方程整理成矩阵形式。定义en×1节点碳势向量An×n碳势传递矩阵其中A(i,j)表示节点j向节点i注入的有功功率占节点i总注入功率的比例即A(i,j) Pji / Pin,i当j∈S_in,i时否则为0bn×1右端向量b(i) (Σ PG,g·eG,g) / Pin,i即节点i自身发电带来的碳势贡献。则方程组为$$ (\mathbf{I} - \mathbf{A}) \mathbf{e} \mathbf{b} $$在MATLAB中求解只需要一行代码e (eye(n) - A) \ b;提示这个矩阵通常是稀疏的节点数多时建议用sparse构造效率会高很多。14节点规模无所谓但养成好习惯总没错。2.4 网损的处理问题这里有一个学术界都会讨论的细节网络损耗怎么处理严格来说支路首端和末端的功率并不相等中间差一个网损。如果直接把首端功率用于节点碳势计算会造成约几个百分点的误差。常见做法有两种忽略网损视每条支路首末端功率相等。IEEE 14节点系统总网损约占总发电的4%~5%这种近似对定性结论没影响很多EI论文也会这么处理。以末端注入功率为准潮流计算会同时给出支路首端功率PF和末端功率PT判断实际功率方向后用流入节点的那个功率值参与计算。我自己的实现偏向于这种代码逻辑也不复杂本文代码采用的就是这个思路。3. IEEE 14节点系统建模与运行场景设计3.1 系统结构与参数说明IEEE 14节点标准系统包含14个节点、20条支路和5台发电机组基准功率100MVA。发电机分布在节点1、2、3、6、8其中节点1是平衡节点。变压器支路集中在4-7、4-9、5-6三处连接230kV和69kV两个电压等级。使用MATLAB做潮流计算最方便的工具是Matpower。它自带case14数据文件调用方法非常简单mpc loadcase(case14);这里我提醒一句如果你是从零开始造数据而不是用Matpower自带的case14一定要注意变压器变比和并联导纳的格式。手工录入一旦漏了变比潮流结果和标准结果会相差很大碳流算出来自然也是错的。3.2 发电机碳排放强度的设定IEEE 14节点系统本身没有给出发电机燃料类型实际使用中需要根据自己的研究场景定义。我这里设定一组典型值用于展示碳流在网络中的传播特性节点机组类型碳排放强度 (g/kWh)备注1燃煤机组950全网主力电源2燃气机组450清洁替代3生物质/自备电厂150中低碳排放6水电机组0零碳8风电机组0零碳不同文献采用的具体数值会有差异比如有的用《中国电网平均CO2排放因子》那套标准但计算流程完全一致你只需要替换eG向量。3.3 运行工况与Matpower配置为了让碳流计算的结果更有分析价值我会设定一个包含不同电源参与的运行方式。假设某次运行中节点1燃煤机组承担约170MW出力节点2燃气机组出力50MW节点3生物质机组出力20MW节点6水电机组出力10MW节点8风电机组出力20MW。这样全网总发电约270MW负荷约259MW网损约4%。碳排放强度差异明显的多种电源同时接入碳流在网络里的混合过程会体现得比较充分。Matpower求解潮流的代码如下opt mpoption(verbose, 0, out.all, 0); res runpf(mpc, opt); bus res.bus; branch res.branch; gen res.gen; n size(bus, 1); nb size(branch, 1); % 节点编号 from_bus branch(:, 1); to_bus branch(:, 2); % 支路首端、末端有功功率来自潮流结果 PF branch(:, 14); % from端有功 PT branch(:, 16); % to端有功 % 发电出力与负荷 Pg_bus zeros(n, 1); for k 1:size(gen, 1) gbus gen(k, 1); Pg_bus(gbus) Pg_bus(gbus) gen(k, 2); end Pd_bus bus(:, 3);注意Matpower的branch结果矩阵中第14列是from端有功功率PF第16列是to端有功功率PT很多新手容易把这两列搞混导致后面碳流方向判断错误。4. Matlab代码实现与关键细节分析4.1 主程序框架整个代码分三步潮流计算、节点碳势求解、支路与负荷碳流计算。主程序结构如下%% 主程序IEEE 14节点系统碳排放流计算 clear; clc; close all; %% 1. 潮流计算 mpc loadcase(case14); % 使用IEEE 14节点标准算例 opt mpoption(verbose, 0, out.all, 0); res runpf(mpc, opt); bus res.bus; branch res.branch; gen res.gen; n size(bus, 1); nb size(branch, 1); from_bus branch(:, 1); to_bus branch(:, 2); PF branch(:, 14); % 首端有功 PT branch(:, 16); % 末端有功 % 聚合节点发电出力 Pg_bus zeros(n, 1); for k 1:size(gen, 1) gbus gen(k, 1); Pg_bus(gbus) Pg_bus(gbus) gen(k, 2); end %% 2. 发电机碳排放强度g/kWh % 按节点顺序定义1-燃煤, 2-燃气, 3-生物质, 6-水电, 8-风电 eG_bus zeros(n, 1); eG_bus(1) 950; eG_bus(2) 450; eG_bus(3) 150; eG_bus(6) 0; eG_bus(8) 0; %% 3. 求解节点碳势 [cbus, Pin] calcNodeCarbonPotential(n, from_bus, to_bus, PF, PT, Pg_bus, eG_bus); %% 4. 计算支路碳流率与负荷碳流率 branch_carbon zeros(nb, 1); for k 1:nb f from_bus(k); t to_bus(k); if PF(k) 0 % 功率从f流向t支路碳流密度节点f的碳势 branch_carbon(k) PF(k) * cbus(f); else % 功率从t流向f支路碳流密度节点t的碳势 branch_carbon(k) abs(PF(k)) * cbus(t); end end % 负荷碳流率 负荷功率 × 节点碳势 load_carbon bus(:, 3) .* cbus; %% 5. 结果输出 disp(节点碳势 (g/kWh):); for i 1:n fprintf(节点 %2d : %8.2f g/kWh, 负荷 %8.2f MW, 负荷碳流 %10.2f kg/h\n, ... i, cbus(i), bus(i,3), load_carbon(i)); end %% 6. 绘制节点碳势柱状图 figure; bar(cbus, FaceColor, [0.2 0.5 0.8]); xlabel(节点编号); ylabel(节点碳势 (g/kWh)); title(IEEE 14节点系统节点碳势分布); grid on;4.2 节点碳势求解函数这是整个计算的核心。函数输入潮流结果、发电数据和排放强度输出节点碳势。function [cbus, Pin] calcNodeCarbonPotential(n, from_bus, to_bus, PF, PT, Pg_bus, eG_bus) % 计算节点碳势 % 输入: % n - 节点数 % from_bus - 支路首端节点编号 (nb×1) % to_bus - 支路末端节点编号 (nb×1) % PF - 支路首端有功 (nb×1) % PT - 支路末端有功 (nb×1) % Pg_bus - 各节点发电出力 (n×1) % eG_bus - 各节点发电机碳排放强度 (n×1) % 输出: % cbus - 节点碳势 (n×1) % Pin - 节点总注入功率 (n×1) Pin zeros(n, 1); % 节点总注入功率 A zeros(n, n); % 碳势传递矩阵 for k 1:length(from_bus) f from_bus(k); t to_bus(k); if PF(k) 0 % 功率从 f 流向 t注入节点 t 的功率为 PT(k) Pin(t) Pin(t) PT(k); A(t, f) A(t, f) PT(k); elseif PF(k) 0 % 功率从 t 流向 f注入节点 f 的功率为 abs(PF(k)) p_in abs(PF(k)); Pin(f) Pin(f) p_in; A(f, t) A(f, t) p_in; end % 若 PF(k)0则该支路不参与碳流传递 end % 加入本地发电注入 Pin Pin Pg_bus; % 构造方程组 (I - A) * e b b zeros(n, 1); for i 1:n if Pin(i) 1e-6 % 归一化传递矩阵 A(i, :) A(i, :) / Pin(i); A(i, i) 0; % 自身项归零统一用单位阵处理 if eG_bus(i) 0 Pg_bus(i) 0 b(i) Pg_bus(i) * eG_bus(i) / Pin(i); else b(i) 0; end end end M eye(n) - A; cbus M \ b; end这里有个细节值得单独说A(i,i)要先清零再求解。原因是某些节点可能通过多条路径间接和自己形成环流如果在归一化时不把对角线处理干净单位阵减完之后对角线可能出问题导致病态矩阵。清零后再用(M I - A)求解物理含义才是对的。4.3 为什么这样处理网损我在函数里用的是“功率实际流入节点的那一端”的值。具体说如果PF(k)0功率从f流向t那么注入t节点的功率应该用末端功率PT(k)如果PF(k)0功率从t流向f注入f节点的功率用|PF(k)|。这样处理的好处是每个节点收到的功率和它送出的功率严格满足基尔霍夫定律不会因为网损而产生“凭空多出来”或“凭空消失”的碳流最终的总量守恒校验能对得很平。这也是我在多次实践后觉得最稳妥的一种做法。4.4 结果可视化除了柱状图我一般还会画一个碳流分布的棒棒糖图把支路碳流率画在系统拓扑上。简单的做法% 准备拓扑坐标按case14常见坐标近似 busXY [ 0 0; 1 1; 1 2; 1 1.5; 0 1; 0 2; 0.5 2.2; 1 2.3; 1.5 2; 2 2; 2 1.5; 2 1; 1.5 0.5; 2 0; ]; figure; hold on; for k 1:nb f from_bus(k); t to_bus(k); x [busXY(f,1), busXY(t,1)]; y [busXY(f,2), busXY(t,2)]; % 线宽反映碳流率大小 width 1 branch_carbon(k) / max(branch_carbon) * 5; plot(x, y, LineWidth, width, Color, [0.8 0.2 0.2]); end scatter(busXY(:,1), busXY(:,2), 100, filled, MarkerFaceColor, [0.2 0.2 0.7]); for i 1:n text(busXY(i,1)0.05, busXY(i,2)0.05, sprintf(%d, i), FontSize, 10); end axis equal; axis off; title(IEEE 14节点系统支路碳流分布);这套可视化代码适合快速定性观察如果投论文用建议用plot结合地图工具箱自定义布局效果会更好。5. 计算结果分析与验证逻辑5.1 节点碳势分布结果以第3节设定的运行场景为例计算出各节点碳势如下表说明性数值具体以你复现时的潮流结果为准节点节点碳势 (g/kWh)负荷 (MW)负荷碳流 (kg/h)1950.000.00.02861.1221.767.33733.4594.2248.84751.2047.8129.35784.317.621.46534.6711.221.67502.140.00.08472.330.00.09445.0229.547.310418.769.013.611392.563.54.912376.816.18.313365.4413.517.814358.1014.919.2看到这个表明显能看出几个规律。节点1的碳势就是燃煤机组的950 g/kWh因为平衡节点直接由本地电源供电上游没有其他来源。节点2的碳势降到861因为本地燃气机组450 g/kWh和节点1送来的功率混合拉低了整体碳势。越往网络末端接入的水电、风电功率份额越大节点碳势逐步走低到节点14已经降到358左右。这说明什么说明碳流计算能把“绿色电源的贡献沿着网络路径逐步稀释高碳电源影响”这个过程量化地刻画出来而这正是传统平均因子法做不到的。5.2 碳流守恒验证算完碳流之后必须做一次守恒校验否则结果很可能有问题。校验逻辑很简单全网发电侧总碳排放应该等于所有负荷侧碳流加上网损对应的碳排放。发电侧总碳排total_gen_carbon sum(Pg_bus .* eG_bus);负荷侧总碳排total_load_carbon sum(load_carbon);网损对应的碳排放可以这样算每条支路的网损乘以该支路上的碳流密度再加总。不过按本文的实现方法由于支路末端功率已经用于节点注入计算网损的碳“消耗”其实隐含在负荷侧的分配里。更直观的校验方式是看碳流总量平衡$$ \sum_{\text{发电}} P_{G,g} \cdot e_{G,g} \approx \sum_{\text{负荷}} P_{L,i} \cdot e_i \sum_{\text{网损}} P_{\text{loss},k} \cdot \rho_k $$实际计算时由于我们是用末端功率做节点注入等式两边会有较小的偏差偏差率一般在1%以内。如果偏差明显优先检查支路方向判断是否和潮流方向一致发电机碳强度赋值是否正确对应到节点是否存在某条支路PF(k)0但实际有功率的情况收敛精度问题。5.3 从结果中能读出什么从碳流结果里我们能得到几个对实际工程有意义的结论第一网络拓扑对碳排放责任分摊的影响很大。同样是负荷节点接在高碳电源附近的节点碳势高接在清洁电源附近的节点碳势低这就是“谁用电、用谁的电”的量化体现。第二中间节点对碳势有“混合”作用。节点6、7、8这一片由于风电和水电的注入碳势明显低于节点1、2、3的高碳区域说明清洁电源对整个系统碳排放的稀释作用是通过电网路径逐步传播的。第三负荷碳流率可以直接用于碳标签。比如节点3的负荷是94.2MW、碳势733 g/kWh对应每小时产生约248.8kg碳排放责任。这个数值可以直接写进产品的碳足迹报告里。6. 复现过程中容易踩的坑与后续扩展方向6.1 潮流方向与符号处理这是我在这么多位初学者身上见到最多的问题。Matpower的PF为正表示功率从from端流向to端为负表示从to端流向from端。构建节点碳势方程时必须根据正负判断哪个节点是“流入侧”。有人会直接把所有支路的PF都当正值加到to_bus上结果负方向支路多的区域碳势完全错误。这个问题的典型症状是守恒校验对不上或者某些节点碳势出现大于所有发电机组碳强度的值比如超过950。遇到这种结果第一反应就应该是去查支路方向。6.2 多发电机组节点的数据聚合IEEE 14节点每台机组单挂一个节点但真实电网中一个节点可能挂多台机组。这时需要把同一节点所有机组的出力相加碳强度用出力加权平均eG_bus_i sum(Pg_g .* eG_g) / sum(Pg_g);如果不做加权平均直接把多台机组碳强度简单平均甚至取最大值都会导致节点碳势计算偏差。这个坑我在处理IEEE 30节点算例时踩过一开始死活对不上文献结果最后发现就是这里的问题。6.3 潮流工具版本兼容性Matpower不同版本之间case14数据文件的字段编号是稳定的但如果用的是老版本自定义数据要注意单位问题。有些早期数据文件以标幺值存储需要乘以基准容量才能得到MW单位的功率值。如果发现碳流数量级明显不对比如负荷碳流变成几百吨每小时先检查单位换算。6.4 扩展方向与深入思考这套方法在14节点上跑通之后往实际工程方向扩展有四个比较顺的路一是时序碳流计算。把全天96点或8760小时的潮流逐时段计算得到节点碳势的时间序列用于分析用户侧碳排放的日内波动。这个方向在碳电耦合分析里很常见实现上就是包一层时间循环。二是碳流灵敏度分析。在节点碳势的基础上对发电出力和网络参数求偏导能算出来“某台机组加出力10MW对某个负荷节点的碳势影响多少”这对碳交易定价和机组调度优化很有价值。三是考虑网络约束的碳流优化。碳流计算本身是潮流结果的后续处理但如果要做低碳调度需要把碳排放约束嵌入最优潮流这时碳势不再是固定值而是随运行方式变化的变量计算复杂度会上一个台阶。四是区域间碳排放责任分摊。把大电网拆成多个区域通过联络线上的碳流率量化区域之间的碳转移这对跨省跨区电力交易中的碳责任划分有直接参考意义。我自己在实际操作中的体会是碳排放流的代码实现不算难难的是把每一步的物理含义想清楚。特别是网损处理、方向判断这两个地方多花十分钟把逻辑捋顺后面能省下三个小时的调试时间。希望这篇东西能帮准备复现这个方向论文的朋友们把路走直一点。