搞电力系统的人对潮流计算都不陌生但如果你突然问我一句“一条110kV线路上面这几十兆瓦里到底对应多少吨二氧化碳”传统潮流就答不上来了。这正是电力系统碳排放流方法要解决的核心问题。简单说碳排放流就是给电网里的有功功率“染色”把发电机侧的碳排放属性跟着潮流一路往下游传播最后落到负荷侧。这篇博文我要复现的就是用Matlab在IEEE 14节点系统上实现完整的碳排放流计算输出节点碳势、支路碳流和负荷碳流三类结果。适合正在做碳排放核算、电力市场碳责任分摊、或者写相关方向论文的学生和工程师参考代码思路可以直接迁移到更大规模的节点系统上。1. 到底要算什么碳流的三类核心指标1.1 碳流不是“烟囱排放量”而是电网里的“碳排放浓度”很多人第一次接触碳排放流时容易绕进去觉得“发电厂排放多少二氧化碳直接加总不就完了吗”对单一发电厂来说确实是这样但电网是一个多电源共同供电的网络。同一个负荷节点它的电可能来自煤电、气电也可能来自水电、风电不同电源的碳排放强度差别很大。我们没法把物理上已经混合在一起的电再拆开告诉大家“你这度电是哪个电厂发的”但通过碳排放流方法可以在能量流跟踪的意义上把碳排放合理地分配到每条支路和每个负荷头上。这里必须区分三个核心量节点碳势单位是kg/MWh或者g/kWh表示某个节点每消耗1MWh电能所对应的碳排放量。可以理解成电网中某个节点的“碳排放浓度”。支路碳流单位是kg/h或者t/h表示某条输电线路上单位时间流过的碳排放量等于该支路有功功率乘以首端节点碳势。负荷碳流单位同样是kg/h或t/h表示某个负荷节点单位时间用电量对应的碳排放等于该节点负荷功率乘以节点碳势。如果把电网想象成一条水管系统有功功率是水流碳排放就是水里的“盐浓度”。发电厂相当于不同浓度的盐水源节点相当于水池支路相当于管道负荷相当于用水端。碳排放流计算做的事情就是根据每个水池的进水盐浓度和水量算出每个水池里的盐浓度再算每根管道流向下游的盐总量。这里有个关键前提碳排放流是基于有功潮流的无功功率不承担碳排放传递功能。这是物理上说得通的——真正做功的是有功功率无功功率只是维持电压的一种交换功率把它纳入碳流计算只会让结果失真。1.2 为什么选IEEE 14节点系统和Matlab做复现IEEE 14节点测试系统是电力系统领域最经典的算例之一14个节点、20条支路、5台发电机、3台变压器。规模不大不小手动检查结果很方便但又不是三节点那种过于理想化的案例。它覆盖了两种电压等级138kV和69kV、多个负荷节点和多条环网支路用来验证碳流算法再合适不过。Matlab的优势不用多说矩阵运算是它的看家本领而碳流计算的核心正好就是解一个线性方程组。搭配MATPOWER工具箱一行代码就能跑完IEEE 14节点的交流潮流然后我们只需要在潮流结果基础上做矩阵组装和求解即可。整个复现过程的核心代码量其实非常少去掉注释也就五六十行。选择这个案例还有一层原因IEEE 14节点系统的潮流结果公开可得任何人跑出来的节点电压、支路功率都可以互相校核。这意味着你复现的碳流结果如果和别人论文对不上问题一定出在代码逻辑上不太可能是因为原始数据不同。这一点对初学者非常友好。2. 数学原理与计算流程拆解2.1 节点碳势方程的推导碳排放流最核心的方程是节点碳势方程。先看任意节点k它在某一时刻同时存在多种功率注入来源上游支路注入功率、本地发电机出力同时它又往外送出功率并供应本地负荷。在稳态潮流下流入功率等于流出功率但碳排放不像功率那样“无条件守恒”它需要乘上每个来源对应的碳强度。设节点k的碳势为e_k单位为kg/MWh。那么节点k单位时间“带走”的碳排放可以分为两部分本地负荷消耗的部分P_load,k * e_k以及从节点k流向相邻节点的支路功率所携带的碳排放Σ P_k→j * e_k。节点k的碳排放来源则包括上游节点i流入支路功率携带的碳流Σ P_i→k * e_i以及本地发电机出力对应的碳排放P_gen,k * c_gen。于是可以写出碳平衡方程e_k * (P_load,k Σ P_k→j) Σ (P_i→k * e_i) P_gen,k * c_gen写成矩阵形式就是A * e b其中矩阵A的对角元素是节点k的总流出功率支路流出功率 本地负荷非对角元素是支路流入节点的功率取负号向量b是各节点本地发电机的碳排放注入量。求解这个线性方程组就得到了所有节点的碳势。这里需要特别强调的是支路功率方向必须根据潮流计算结果判断。不能只看支路的参考方向。如果某条支路的实际潮流方向和参考方向相反就必须把流入节点和流出节点对调。这一点在做代码时非常容易出错。2.2 支路碳流、负荷碳流和网损的处理节点碳势求解完成后其他两类结果就是简单的乘法支路碳流F_l |P_flow| * e_from其中P_flow是该支路实际输送功率e_from是实际送端节点的碳势。负荷碳流F_load,k P_load,k * e_k。但有一个绕不开的问题网损。交流潮流的线路上送端功率和受端功率不相等差值就是线损。线损这部分功率怎么算碳目前工程上最常见的处理方法有两种一种是把线损按比例分摊到送端和受端另一种是忽略线损直接用支路首端功率乘首端碳势作为支路碳流。后者简单直观在网损率不高的情况下误差很小IEEE 14节点系统的网损率只有百分之几完全在可接受范围内。我在复现时采用的是首端功率法即直接使用潮流结果中支路首端的有功功率P_fr乘首端节点碳势得到该支路碳流。这样碳流在全网存在一个很小规模的“损耗缺口”可以单独作为网损碳流项在做全网碳平衡校核时把它考虑进去。2.3 整体计算流程整个碳排放流计算可以拆成六个步骤用MATPOWER或其他潮流工具求解IEEE 14节点系统的交流潮流。从潮流结果中提取每个节点的负荷功率、每台发电机的出力和对应节点编号。根据潮流方向判断每条支路的实际送端和受端。组装矩阵A和向量b求解线性方程组得到节点碳势e。根据节点碳势计算支路碳流和负荷碳流。进行全网碳排放平衡校核检查发电侧总排放是否等于负荷碳流加网损碳流之和。这套流程不依赖任何专用碳流工具箱纯手工搭建优点是逻辑透明改到其他节点系统时只需要换掉case数据不需要改代码框架。3. Matlab代码实现手把手复现3.1 准备工作数据与潮流求解假设你已经安装好Matlab并且把MATPOWER的路径加进来了。整个复现的第一步是加载IEEE 14节点系统并求解潮流define_constants; mpc loadcase(case14); res runpf(mpc); bus res.bus; branch res.branch; gen res.gen; n size(bus, 1);这里define_constants会定义MATPOWER里常用的列索引常量比如PG、PF、F_BUS方便后面按名字取列。潮流求解完成后res.gen里存的是各台发电机的出力res.branch里存的是各支路首端和末端的功率单位都是MW。下一步把这些数据提取出来gen_bus gen(:, GEN_BUS); Pg gen(:, PG); % MW Pd bus(:, PD); % MW有一点要提醒MATPOWER潮流结果中的GEN_BUS是节点编号IEEE 14节点系统的节点编号正好是1到14所以后面可以直接把这个编号当作矩阵索引用。如果你换到其他系统节点编号不连续时最好先做一个编号到矩阵索引的映射表。3.2 构建碳流方程矩阵现在进入核心部分构建矩阵A和向量b。先初始化稀疏矩阵和几个辅助变量A sparse(n, n); b zeros(n, 1); outflow zeros(n, 1);然后设置各台发电机的碳排放强度。这里的取值可以根据研究场景调整我先用一组典型值c_gen zeros(n, 1); c_gen(1) 400; % 节点1火电机组单位kg/MWh对应400g/kWh c_gen(2) 300; % 节点2气电机组对应300g/kWh % 节点3、6、8的机组在默认潮流中只发无功有功为0碳排放注入为0接下来遍历所有支路判断每条支路的实际潮流方向更新outflow向量和矩阵A的非对角元素for l 1:size(branch, 1) i branch(l, F_BUS); j branch(l, T_BUS); P_fr res.branch(l, PF); if P_fr 0 % 实际潮流从i流向j outflow(i) outflow(i) P_fr; A(j, i) A(j, i) - P_fr; else % 实际潮流从j流向i P_back -P_fr; outflow(j) outflow(j) P_back; A(i, j) A(i, j) - P_back; end end这里的关键点在于非对角元素取负号对角元素是总流出功率加本地负荷。等到所有支路遍历完最后给A补上对角元素for k 1:n A(k, k) outflow(k) Pd(k); end发电机碳排放注入向量b的组装非常简单for k 1:size(gen, 1) b(gen_bus(k)) b(gen_bus(k)) Pg(k) * c_gen(gen_bus(k)); end注意这里的单位运算Pg单位是MWc_gen单位是kg/MWh乘积就是kg/h正好是碳排放流量单位。3.3 求解并输出三类碳流结果矩阵A和向量b组装完成节点碳势就是解线性方程组e A \ b; % 节点碳势单位kg/MWhMatlab的稀疏矩阵左除用的是非常成熟的求解器14节点这种规模几乎瞬间完成。拿到节点碳势后支路碳流和负荷碳流就是逐项计算branch_carbon zeros(size(branch, 1), 1); for l 1:size(branch, 1) i branch(l, F_BUS); if res.branch(l, PF) 0 from_node i; else from_node branch(l, T_BUS); end branch_carbon(l) abs(res.branch(l, PF)) * e(from_node); % kg/h end load_carbon Pd .* e; % kg/h为了便于查看把结果汇总成表格T_node table((1:n), e, Pd, load_carbon, ... VariableNames, {节点, 碳势_kg每MWh, 负荷_MW, 负荷碳流_kg每h}); disp(T_node);支路碳流表可以按碳流从大到小排序看看哪条支路是“碳流主干道”T_branch table(branch(:, F_BUS), branch(:, T_BUS), ... abs(branch(:, PF)), branch_carbon, ... VariableNames, {首端, 末端, 有功_MW, 碳流_kg每h}); T_branch sortrows(T_branch, 碳流_kg每h, descend); disp(T_branch);3.4 结果可视化只出一堆数字不够直观我习惯再加两个图。节点碳势用柱状图一眼扫过去就能看出全网碳势的高低区间figure; bar(1:n, e); xlabel(节点编号); ylabel(节点碳势 (kg/MWh)); title(IEEE 14节点系统节点碳势分布); grid on;支路碳流分布可以用线宽来映射碳流越大的支路画得越粗这样可以直观看出碳流从电源侧往负荷侧的“流动骨架”figure; hold on; for l 1:size(branch, 1) i branch(l, F_BUS); j branch(l, T_BUS); x [bus(i, 6), bus(j, 6)]; % MATPOWER的bus矩阵第6列不是坐标 y [bus(i, 7), bus(j, 7)]; % 实际case14的bus矩阵没有坐标列建议手动定义坐标数组 plot(x, y, LineWidth, 1 10 * branch_carbon(l) / max(branch_carbon)); end hold off;这里有个小坑MATPOWER的case14格式里bus矩阵第6列是区域编号第7列是电压基准值并不是节点坐标。要画网络拓扑图需要另外定义一份节点的可视化坐标。IEEE 14节点系统的经典坐标网上很容易找到也可以根据单线图自行设置不影响计算逻辑。4. IEEE 14节点算例结果与分析4.1 节点碳势分布特征用默认的case14数据和上面这组碳排放强度跑完后节点碳势会呈现一个非常清晰的分布规律。节点1的碳势就是它本地发电机的碳强度400kg/MWh因为节点1只有一台发电机直接接入没有上游支路给它注入其他碳源。节点2的碳势会略低于节点1但高于节点2本地发电机的300kg/MWh原因是节点2既有本地300的气电注入又从节点1方向吸收了碳势更高的功率两者混合后节点碳势大约在330到350kg/MWh之间。这个“混合后介于两个来源之间的数值”正是碳流方法最核心的表达能力。再往负荷密集的区域看节点3、节点4、节点5这类纯负荷或轻负荷节点的碳势是由多个上游节点的碳势加权平均得到的。因为IEEE 14节点系统中环网支路不少比如节点2-4-5这一圈碳流会同时从多个方向汇入节点碳势体现的是“多源供电路径下的综合碳排放浓度”这是单看发电机碳排放强度完全得不到的信息。节点8的情况值得单独说。节点8的本地机组在默认潮流中有功出力为0只发无功但它通过变压器支路7-8从节点7吸收有功功率所以节点8的碳势不为0而是等于节点7的碳势。这说明碳流方法天然把“虚拟发电机”和“同步调相机”这类不发出有功的装置排除在碳排放贡献之外非常符合物理直觉。4.2 支路碳流和负荷碳排放核算从支路碳流的排序结果看最大的几条支路通常集中在两个区域一是发电机节点向主网送电的通道比如节点1到节点2、节点1到节点5二是连接高压网络和低压负荷区的关键断面。支路碳流的大小区分于两个因素一是功率大小二是首端碳势。一条功率巨大但首端碳势只有100的支路碳流未必比一条功率适中但首端碳势400的支路大。这也是碳排放流方法比单纯看潮流更有价值的地方。负荷侧的结果同样有看点。IEEE 14节点系统中负荷最大的是节点394.2MW但它的碳势在全网处于中等水平所以它的负荷碳流不是全网最大。节点9、节点14这类低电压等级的负荷节点虽然单个负荷不大但碳势往往会受到上游多路径碳流汇聚的影响算出来的单位用电碳排放反而偏高。这种空间差异正是后续做“谁用电、谁排碳”责任分摊的基础。4.3 全网碳排放平衡校核算完之后一定要做一步校验把全网发电侧的总碳排放加起来和负荷碳流加网损碳流加起来对比。发电侧总排放等于各台发电机有功出力乘自身碳强度之和。负荷碳流之和就是sum(load_carbon)。两者之间的差值就是网损对应的碳流也就是支路首端碳流总和与末端碳流总和的差。在默认case14的潮流下节点1出力约232.4MW碳强度400贡献92.96吨每小时节点2出力约40MW碳强度300贡献12吨每小时其他几台发电机有功出力为0不贡献碳排放。全网发电侧总碳排放大约104.96吨每小时。全网负荷总量约259MW平均碳势大约在360kg/MWh上下负荷碳流总量约93到95吨每小时。中间的差值是线损对应的碳排放大约10吨每小时约占发电侧总排放的9%到10%和14节点系统的网损率量级基本吻合。如果这一步校核差得太多就要回头检查方向判断和矩阵组装了。5. 常见问题与调试经验5.1 潮流不收敛或结果异常碳流计算完全依赖潮流结果潮流先发散后面全白搭。用MATPOWER跑IEEE 14节点基本不存在收敛问题但如果你把代码改成其他算例很可能遇到潮流不收敛的情况。这时不要急着怀疑碳流代码先用runpf看看潮流残差检查发电机有功、无功上下限、负荷基准值是否合理。还有一种常见情况是修改了负荷水平后初始潮流点离解太远MATPOWER默认用平启动也就是电压幅值1、相角0大多数场景没问题个别重载场景可以换runpf(mpc, mpoption(pf.alg, NR))手动切换算法。5.2 矩阵奇异与孤立节点构建A矩阵最怕遇到对角元素为0的节点。如果某个节点既没有负荷也没有支路流出功率同时本地发电机出力为0那么这个节点在A中对应对角元素就是0整个矩阵奇异求解直接给出NaN。IEEE 14节点系统设计得比较合理几乎所有节点都有负荷或有流出支路所以默认不会出问题。但如果你把某个负荷改成0就可能踩雷。遇到奇异时先检查该节点是不是“悬空节点”。如果确实是单点孤立可以把矩阵中对应行和向量b中对应行删掉再求解或者用e pinv(A) * b强行求解。不过最稳妥的办法是在建模阶段就给所有节点至少保留一个有效耦合关系避免节点脱离主网。5.3 方向、单位和损耗的坑方向上最大的坑是支路功率超过0的判断。MATPOWER的PF列指的是“从首端流向末端的功率”如果PF为负说明实际功率从末端流向首端。很多人在写代码时直接拿PF乘首端节点碳势一旦系统存在反向潮流结果就全错了。我的处理方式是显式判断方向反向时将送端节点换成末端节点再乘也就是代码里from_node那一行的作用。单位方面最容易搞混的是kg/MWh和g/kWh。记住这两个单位在数值上是相等的也就是说400kg/MWh就是400g/kWh。如果潮流功率用MW、碳强度用g/kWh乘积出来的单位是“克每小时”数值会偏大1000倍。我建议在代码开头统一定义碳强度为kg/MWh计算出来的碳流单位就是kg/h不需要额外换算。损耗方面采用首端功率法后全网会存在一个“碳流损耗缺口”这是正常现象。你要是非要把全网碳流闭环需要额外做损耗分摊常见做法是按各支路负载率把网损碳分摊到下游节点。但对于技术复现和论文验证来说标明“损耗未分摊”完全能接受。6. 从复现到扩展的几个方向6.1 扩展到实时碳流计算这套代码本质上是时断面计算也就是给定一组运行工况求一个碳流分布。如果想做全天24小时的动态碳流只需要把每个时刻的机组出力、负荷水平依次代入潮流求解再循环调用上面的碳流计算函数即可。注意每15分钟或1小时采集一次状态数据时机组的组合和出力都会变化每个时刻都要重新判定支路潮流方向不能复用上一个时刻的方向判断结果。6.2 与机组组合和电力市场结合碳排放流模型从IEEE 14节点往更大系统扩展时瓶颈不在碳流求解本身而在潮流数据的获取。实际电网的拓扑、线路参数、机组出力和负荷曲线往往分布在多个部门先把数据清洗成MATPOWER的case格式才是最大的工作量。另外碳流结果可以直接作为机组碳排放责任分摊的依据也可以结合电力市场价格分析不同节点的“碳排放足迹”会不会影响边际电价这是目前论文里比较热门的方向。我在实际使用这套方法时最满意的一点是它足够透明每一步都能手工验算不像某些所谓的碳流分析平台黑箱跑出一个数你根本不知道中间发生了什么。把Matlab代码吃透之后换成任何节点系统无非是改case数据、改碳强度设置核心框架完全不动。最后提醒一句碳流结果的使用一定要标注计算方法和边界条件不同文献里支路损耗处理方式不同直接对比数值没有意义。