基于潮流结果的电力系统碳排放流计算：IEEE 14节点Matlab复现全解析

最近在做一个电网碳排放核算相关的项目翻了不少EI论文发现“电力系统碳排放流”这个词出现频率特别高。无论是做源网荷协调优化、碳追踪还是算负荷侧碳责任大家都在用这个方法。我前阵子把一个最经典的版本——基于潮流结果的碳排放流计算方法在IEEE 14节点系统上用Matlab完整复现了一遍。这篇文章就把整个思路、矩阵推导、代码实现和踩坑记录全部写出来正好给正在看论文但被公式卡住、或者想在自己算例上做碳流分析的读者一份可以直接照做的参考。1. 碳排放流是个什么东西为什么大家都在算1.1 从“发电排放”到“用电责任”传统意义上我们核算电力系统碳排放基本都是按发电厂的口径统计某台火电机组装机多少、发了多少电、烧了多少煤乘一个排放因子就得出这家电厂排了多少碳。这种“源头核算”方法本身没问题但它回答不了几个很现实的问题电网里流动的电能到底是谁发的用户用的每一度电到底对应多少碳排放西部风电场发的清洁电送到东部东部负荷的碳排放到底算谁的碳排放流就是用来回答这些问题的。它的核心思想很朴素把发电侧产生的“碳排放”当成一种伴随有功功率流动的物质顺着输电线路、变压器一级一级往下游传播最后分摊到每个负荷节点和每条支路上。这样就能算出每个节点的碳势相当于每个节点的“碳浓度”、每条支路的碳流率相当于“碳流量”以及每个负荷承担了多少碳责任。这套方法和电力系统里最成熟的潮流计算是天然配套的潮流算出功率怎么流碳流就跟着功率怎么走。所以碳排放流本质上是一种基于潮流的后处理分析工具不会反过头去影响系统运行方式计算量也小非常适合做碳排放核算、碳追踪、碳责任分摊这些场景。1.2 为什么拿IEEE 14节点当试验田很多刚接触碳排放流的人会问为什么大家都拿IEEE 14节点系统做复现而不是直接上IEEE 118或者某省实际电网原因有几个。第一IEEE 14节点规模小但有代表性14条母线、5台发电机、20条支路线路加变压器拓扑不简单到一眼看穿也没有复杂到手工算不动刚好能把碳流计算的全流程走通。第二这个算例是公开标准算例Matpower里自带case14数据任何人下载安装Matpower就能直接跑复现门槛极低。第三EI论文里做碳排放流验证的一半以上都用14节点系统作为基础算例复现它等于和文献里的结果有了一个“公共坐标系”后续换更大系统或者改算法对比起来方便。我在实际做的时候还特意把case14里默认的零出力机组节点3、6、8改成有出力用来模拟多电源混合输送的场景。这样算出来的碳流分布更有意思能明显看到清洁电源对下游节点碳势的“稀释”作用。2. 数学模型怎么把碳流算出来2.1 三个基础假设先立规矩碳排放流计算不是凭空来的它建立在几个假设之上。理解这几个假设后面看公式才不迷糊。第一个假设碳排放流只伴随有功功率流动无功功率不承担碳流传递。这个其实很好理解碳排放本质上是和电能生产量挂钩的无功功率不产生实质电能量所以碳流只沿有功潮流路径传播。第二个假设在同一个节点上来自不同电源的碳排放完全混合节点上所有流出功率包括负荷和转出支路的碳势相同。就好比两杯不同浓度的糖水倒进一个杯子里搅匀从这杯水里倒出去的任何一杯糖水浓度都一样。这是碳排放流的“共享分配”原则也是它和“逐源追踪”类方法的最大区别。第三个假设碳排放流在网络中是守恒的。发电机注入多少碳流最终必须全部分摊到负荷和网络损耗上去不多不少。基于这个守恒关系才能在每个节点建立碳流流入等于流出的平衡方程。这三个假设决定了方法的基本框架。实际工程中还有一些扩展做法比如把网损按比例分摊到各负荷或者引入分段碳排放流来处理时序问题但核心思路都离不开这三点。2.2 两个核心变量节点碳势和碳流率碳排放流体系里的核心变量说到底是两个。第一个是节点碳势记作 e_i单位一般是 kgCO2/kWh。直观理解就是节点 i 上“每度电对应的碳排放量”。它有点像电力系统里的电压幅值是整个碳流计算里最关键的中间量。只要所有节点的碳势求出来后面所有碳流率都唾手可得。第二个是碳流率记作 R单位一般是 kg/h 或 t/h。它表示单位时间内伴随功率流动通过某一断面或注入某节点的碳排放量。具体又分几类支路碳流率支路上传播的碳流量、负荷碳流率负荷从电网“取走”的碳流量、发电碳流率发电机注入电网的碳流量。碳流率的大小等于对应有功功率乘以该功率携带的碳势所以它是一个直接反映“碳怎么流”的量。把两个变量放在一起看发电机往节点注入功率 PG 的同时也注入了一个碳流率 PG × e_Ge_G 是发电机的碳排放强度也就是单位发电量对应的排放节点上所有流出功率无论去负荷还是去支路都按节点碳势 e_i 携带碳流。2.3 从两节点推起看懂矩阵方程来源矩阵公式 F × e C 看起来高大上但它的来源用两节点系统就能推得明明白白。假设发电机接在节点1出力 PG1碳排放强度 e_G1节点1通过一条线路向节点2送有功 P12节点2接了一个负荷 PL2。先看节点1流入节点1的功率只有发电机出力 PG1没有其他支路注入所以节点1的碳势就是发电机的碳排放强度即 e_1 e_G1。再看节点2流入节点2的功率是支路有功 P12流动过程中携带的碳流率是 P12 × e_1等于 P12 × e_G1流出节点2的功率是负荷 PL2碳流率是 PL2 × e_2。根据流入等于流出的守恒关系P12 × e_G1 PL2 × e_2所以e_2 (P12 / PL2) × e_G1注意如果 P12 大于 PL2说明线路上存在网损那么 e_2 会比 e_G1 大。这个结果其实非常有物理意义网损相当于白白消耗了一部分碳流剩下的碳流由更少的负荷功率承担所以负荷侧的“碳浓度”被抬高了。这也是碳排放流里一个很重要的结论——线损客观上会增加下游用户的碳责任。把这种“节点流入碳流等于流出碳流”的平衡关系对每一个节点都写出来就组成一个线性方程组。把所有节点功率关系整理成矩阵形式就是后面代码里要解的F × e C其中 F 是节点有功通量矩阵对角线元素是节点总注入功率发电机出力加支路流入非对角线元素是负的支路流入功率C 是各节点的发电碳流率向量等于节点上发电机出力乘对应碳排放强度。解这个方程组一次性得到所有节点的碳势 e。3. Matlab复现全过程3.1 环境准备MATPOWER装好就成功了一半碳排放流计算的第一步是拿到准确的潮流结果。自己写牛顿-拉夫逊潮流代码不是不行但完全没有必要——学术界做电力系统分析事实标准的工具是MATPOWER开源、免费、精度高而且自带IEEE 14节点标准数据。MATPOWER的安装非常简单去官网下载压缩包解压后把文件夹路径添加到Matlab路径或者在Matlab里直接运行文件夹里的install_matpower脚本。装完之后在命令窗口敲 case14能输出一个14节点的数据体就说明安装成功了。我建议再顺带跑一下自带测试mpc case14; res runpf(mpc);如果 runpf 返回的 res 结果里 c5conv 字段是1表示收敛环境就完全OK了。后面所有碳流计算都以 res 里保存的潮流结果作为输入。3.2 跑潮流和提取数据MATPOWER的 runpf 返回结果是一个结构体里面包含 bus、gen、branch 三个最关键的子表。bus 表的第3列是节点有功负荷MWgen 表的前两列是发电机所在节点和出力MWbranch 表的第1、2列是支路首端和末端节点编号第14列是首端有功潮流 Pf第16列是末端有功潮流 Pt。这里的正负号约定特别重要一定要先弄清楚Pf 表示从节点 f 注入线路的有功功率Pt 表示从节点 t 注入线路的有功功率。如果 Pf 大于0说明实际功率从 f 流向 t如果 Pf 小于0说明实际功率从 t 流向 f。Pt 的正负同理只是相对 t 节点而言。很多复现翻车都翻在这一步下面代码里会刻意按这个约定处理。提取基础数据的代码如下mpc case14; res runpf(mpc); bus res.bus; gen res.gen; branch res.branch; nb size(bus, 1); % 节点数 nl size(branch, 1); % 支路数 % 节点发电出力聚合可能有多个发电机在同一节点 PG zeros(nb, 1); for k 1:size(gen, 1) gbus gen(k, 1); PG(gbus) PG(gbus) gen(k, 2); end % 节点有功负荷 PD bus(:, 3);这段代码里给每个节点的所有发电机出力做了聚合因为case14里一台发电机对应一个节点但实际系统里一个节点挂多台机很常见聚合逻辑可以直接复用。3.3 构造碳流矩阵的核心代码构造节点有功通量矩阵是碳流计算的核心环节。原理很简单遍历每条支路判断功率实际流向把流入某节点的功率加到该节点的总注入里同时在矩阵对应位置写上负的流入功率。关键代码如下% 发电机碳排放强度单位 kgCO2/kWh按需自行调整 EG_bus zeros(nb, 1); EG_bus(1) 0.95; % 节点1燃煤机组 EG_bus(2) 0.55; % 节点2燃气机组 EG_bus(3) 0; % 节点3清洁电源 EG_bus(6) 0; % 节点6清洁电源 EG_bus(8) 0; % 节点8清洁电源 % 构造节点有功通量矩阵 F以及发电碳流率向量 C F zeros(nb, nb); for k 1:nl f branch(k, 1); t branch(k, 2); pf branch(k, 14); pt branch(k, 16); if pf 0 % 功率方向 f - t流入 t 的实际功率为 -pt inflow_t -pt; F(t, t) F(t, t) inflow_t; F(t, f) F(t, f) - inflow_t; else % 功率方向 t - f流入 f 的实际功率为 -pf inflow_f -pf; F(f, f) F(f, f) inflow_f; F(f, t) F(f, t) - inflow_f; end end % 加上发电机出力到对应节点对角元 for i 1:nb F(i, i) F(i, i) PG(i); end % 发电碳流率向量 C PG .* EG_bus; % 求解节点碳势单位 kgCO2/kWh e_node F \ C;这里面有个值得注意的细节我判断潮流方向只用 pf 的符号但累加节点流入功率时用的是对端功率 pt取绝对值。原因是有功损耗首端送出的功率经过线路阻抗后会损耗一部分实际到达末端节点的功率比首端小。如果用 pf 累加流入节点的功率会导致节点功率不平衡求出来的节点碳势也会失真。用末端实际流入功率 -pt 来构建矩阵等价于把网损从下游节点的流入功率里扣除保证每个节点“流入流出损耗”的物理关系成立数值上更可靠。解线性方程用的还是 Matlab 反斜杠运算符 F \ C14阶矩阵求解是毫秒级的事情完全不用担心性能。计算完成后还可以顺手做支路碳流率和负荷碳流率的后处理% 支路碳流率单位 kg/h branch_carbon zeros(nl, 1); for k 1:nl f branch(k, 1); t branch(k, 2); pf branch(k, 14); if pf 0 branch_carbon(k) pf * 1000 * e_node(f); else branch_carbon(k) (-pf) * 1000 * e_node(t); end end % 负荷碳流率单位 kg/h load_carbon PD .* 1000 .* e_node;这里功率用MW碳势用kgCO2/kWh乘上1000是因为MW等于1000kW最终算出的碳流率单位是kg/h。如果习惯用t/h把结果再除以1000即可。3.4 结果解析与出图用以上代码跑完case14会得到每个节点的碳势值。复现时观察到的典型结果大致是这样电源节点1燃煤机组碳势约0.95 kgCO2/kWh节点2燃气约0.55上下而接有清洁电源的节点3、6、8碳势会明显偏低清洁电源下游的节点碳势也会被拉低。负荷较重的节点如果上游通道较长、损耗较大碳势通常会被抬高这就是前面两节点推导里“网损抬高下游碳浓度”的体现。论文配图一般画几种类型的图。节点碳势用柱状图最直观figure; bar(1:nb, e_node, FaceColor, [0.2 0.5 0.8]); xlabel(节点编号); ylabel(节点碳势 (kgCO2/kWh)); grid on;支路碳流率可以用有向图表现Matlab的graph对象天然支持G digraph(branch(:,1), branch(:,2), branch_carbon); figure; p plot(G, Layout, force, LineWidth, 2); p.EdgeCData branch_carbon; colorbar;线越粗、颜色越深说明这条支路承载的碳流量越大。这种图放在报告里非常直观评审和同事看了都能一眼抓住重点。4. 复现中踩过的坑4.1 矩阵报奇异的三种可能我第一次跑通之前F \ C 直接报矩阵接近奇异warning刷了一大片。排查下来基本是三个原因。第一支路方向判反导致矩阵结构错误。这在新建节点通量矩阵时最容易出现。比如 pf 为负时如果不判断方向仍然往 F(t,t) 里累加功率就会把“流出”当成“流入”矩阵某一行可能加起来小于等于零。判断方向必须以 pf 符号为准不能想当然按支路编号从f到t。第二节点功率不平衡。如果直接拿 case14 原始数据构建矩阵而不先跑 runpf 取潮流结果那么节点注入和流出功率并不严格满足平衡条件矩阵可能是病态的。记住构建碳流矩阵的功率数据必须来自潮流计算结果不是来自算例原始数据。第三孤立节点或者零注入节点处理不当。IEEE 14节点系统本身拓扑是连通的但如果自己改成其他算例可能存在某些节点没有任何功率流入也没有发电机那 F 对应行全为零矩阵必然奇异。标准做法是先检查每个节点的总注入功率对孤岛节点单独处理或直接删掉。4.2 支路潮流方向处理的经典误区这一条我觉得值得单独拿出来讲因为几乎所有初学碳排放流的人都会在这里出错。Matpower里 pf 为正表示从 f 流向 tpf 为负表示从 t 流向 f。但很多人图省事直接取 abs(pf) 然后默认功率是 f 流向 t结果在环网或者功率倒送的支路上碳流方向和实际完全相反节点碳势全是负的。我建议把方向判断和功率累加写成独立函数并且用输入pf符号绝对值双重校验function F add_branch_flow(F, f, t, pf, pt) if pf 0 inflow -pt; % 实际流入节点 t 的功率 F(t, t) F(t, t) inflow; F(t, f) F(t, f) - inflow; else inflow -pf; % 实际流入节点 f 的功率 F(f, f) F(f, f) inflow; F(f, t) F(f, t) - inflow; end end这样封装好以后后面换算例、换数据都不用改主逻辑只改输入就行。4.3 单位搞错的结果有多离谱碳排放流里单位混用是最隐蔽的错误。我见过有人把发电出力读出来是标幺值直接当成有名值参与计算算出来的碳势比正常值差了100倍。IEEE 14节点系统基准容量正好是100MVA所以Matpower里所有功率都是标幺值形式。但 runpf 返回的 bus、gen、branch 表里功率已经是有名值MW了不需要再乘基准容量。真正要注意的是自己在计算发电碳流率时功率用MW碳势用kgCO2/kWh最终碳流率单位是kg/h。如果碳势用g/kWh算出来的结果又会差1000倍。我习惯做完一步就打印一次中间结果用两节点手算值去对照。比如先手动构造一个两节点的简单潮流结果代入代码算一遍确认结果和手推公式一致再上IEEE 14节点。4.4 多发电机节点碳强度聚合如果某个节点挂了多台不同燃料类型的发电机不能简单把每台机组的出力乘碳强度再加起来塞进C向量里面事。正确做法是先算节点总发电碳流率再反推等效碳强度或者干脆不聚合把每个发电机作为独立的碳注入源加入方程。我在代码里用的是先聚合节点发电出力再用等效碳强度 EG_bus(i) 去乘。这种情况下 EG_bus(i) 应该是该节点所有发电机碳排放流率之和除以总出力EF_node sum(gen_carbon_rate_per_gen) / PG(i);如果是做EI论文复现建议把发电机级别的碳流率保留下来算完碳势后再按节点汇总这样既能验证全网碳流守恒又能方便画不同电源贡献占比的图。5. 实操体会和扩展玩法5.1 我的几点个人体会跑通一遍碳排放流计算之后有几个感受特别深。第一碳排放流计算的难度不在于方程本身而在于潮流结果的正确解读。只要功率方向、单位这些细节处理好了核心求解就一行 F \ C。很多论文里搞得神乎其神的公式落到代码上其实很简洁。第二IEEE 14节点系统的碳流结果非常适合做“体检”。把各节点碳势排个序能明显看出清洁电源对局部碳势的拉低作用也能看到网损较大、供电距离较长的节点碳势被抬高了多少。这种结果拿来写分析报告比单纯给一个碳排放总量有说服力得多。第三参数碳强度的设置对结果影响很大。我在代码里给的0.95和0.55只是常用参考值不同文献取值差异很大。复现对比时要先确认对方论文用的排放因子否则数值没有可比性。5.2 往大系统和其他方向扩展这套流程从IEEE 14节点换到IEEE 30、39、118节点代码几乎不用改只要把case14换成对应算例名称就行。真正需要改的是发电机碳排放强度配置因为大系统的发电机类型更多、燃料更杂得逐台查资料或按区域口径统一设置。除了换算例还有几个扩展方向我觉得很有意思一是把碳排放流扩展到时序场景用96点或者8760小时的潮流序列逐时段算碳流再做日、月、年累加得到的是“碳流曲线”二是把碳排放流和最短路算法结合做特定电源的碳流溯源看某个风电场的低碳电到底送去了哪些负荷三是把节点碳势作为目标函数的一部分嵌入到最优潮流模型里做低碳经济调度。我现在手上在跑的一个项目就是把14节点这套流程扩展到一个区域实际网架难点已经不在碳流方程本身而在数据清洗和发电机碳强度口径统一上。如果你也准备往更大规模系统做建议先把F矩阵构造这部分好好封装成函数后面换算例就只是换数据的问题。各位如果在复现过程中遇到其他奇怪的坑欢迎一起交流。