基于IEEE 14节点的电力系统碳排放流计算与Matlab实现

做电力系统方向的研究尤其是涉及“双碳”和低碳转型的课题绕不开一个很实际的问题电网把电能从发电侧送到用户侧中间那一段段错综复杂的潮流里产生的碳排放到底该怎么划分如果只按发电量乘以排放因子来均摊完全不考虑电力输送过程中的网络拓扑和潮流方向那结果别说发论文了连企业层面的责任核算都站不住脚。这也是过去几年电力系统碳排放流Carbon Emission Flow被反复提起的核心原因。这篇文章我打算以IEEE 14节点系统为算例完整走一遍碳排放流从数学模型到Matlab代码实现的全过程。适合正在做碳排放流、电力系统低碳规划、碳追踪方向毕业设计或论文复现的同学参考也适合刚开始接触这个方向、想搞明白“里头的矩阵到底在算什么”的入门者。我会把原理、数据准备、代码骨架、结果校验和常见的坑一次性讲透。1. 碳排放责任不该“一锅端”碳排放流的追责逻辑很多刚开始做这个方向的人会有一个疑问发电厂排放了多少碳不是明摆着的事吗煤电厂的排放因子、出力曲线、燃料消耗量都查得到直接算总排放再除以全社会用电量不就能得到每度电的碳排放了吗这个思路在宏观统计层面没错但放到电力系统内部就完全不够用了。原因在于电网是一个网络不是一根管子。一台燃煤机组发出来的电经过升压、输电、降压最终会混入成千上万个负荷节点的用电里。省内送省外、风电光伏大发时段火电被压出力这些动态过程都意味着同一时刻不同节点上负荷所消耗的电力对应的“边际发电资源”是完全不同的。一个靠近水电基地的节点用电碳强度可能只有靠近火电基地节点的三分之一甚至更低。碳排放流要解决的核心问题就是把“发电侧的总排放”按照电力潮流的物理流动路径逐级映射到每一个节点、每一条支路乃至每一个负荷头上。这种做法在学术上叫“潮流追踪”Power Flow Tracing的一种延伸应用也被称为“按比例分摊原则”——在某一个节点上流入该节点的有功功率按比例混合那么该节点对应的碳排放强度就是所有流入功率碳排放强度的加权平均。用一句话概括它的物理含义碳排放是“贴在电流”上流动的。电流走到哪里碳就走到哪里某节点上所有流入功率中混入了多少高碳电力该节点的碳势就高。这跟传统方法最大的区别是什么传统方法把电网当成一个黑箱黑箱两端是发电和用电做总量平衡碳排放流把电网当成一个透明的有向图每条边的功率、每个节点的碳势都可以独立算出来。这个区别决定了碳排放流能够支撑更高精度的分析比如碳流敏感度分析、低碳调度、碳权分配甚至绿电追踪溯源。2. 从潮流到碳流节点碳势与支路碳流率的数学推导既然要复现算法数学部分不能含糊。下面这套推导是碳排放流方法的核心无论你后续怎么改写、加约束最底层的逻辑都是它。2.1 几个先要明确的定义碳势Carbon Potential记为 (e_i)表示节点 (i) 上单位有功电能所对应的碳排放量单位可以取 (tCO_2/MWh)。它是这个方法的“状态变量”。支路碳流密度记为 (\rho_{ij})表示从节点 (i) 流向节点 (j) 的支路中单位电量携带的碳排放量。根据“流出节点的碳流密度等于该节点碳势”的基本假设有 (\rho_{ij} e_i)。支路碳流率记为 (R_{ij})等于支路有功功率 (P_{ij}) 乘以支路碳流密度即 (R_{ij} P_{ij} \cdot e_i)单位是 (tCO_2/h)。发电机碳注入记为 (E_{G,i})等于发电机有功出力 (P_{G,i}) 乘以该发电机的碳排放因子 (e_{G,i})。2.2 节点碳势方程怎么列对任意节点 (i)从物理意义上看流入该节点的总碳流量应该等于该节点上所有发电机注入的碳流量与所有上游支路流入的碳流量之和。同时节点碳势的定义又要求[ e_i \frac{\text{流入节点 } i \text{ 的总碳流量}}{\text{流入节点 } i \text{ 的总有功功率}} ]假设所有向节点 (i) 注入功率的支路集合为 (N_{in}(i))那么[ e_i \frac{E_{G,i} \sum_{k \in N_{in}(i)} P_{ki} \cdot e_k}{P_{G,i} \sum_{k \in N_{in}(i)} P_{ki}} ]这里 (P_{ki}) 表示从节点 (k) 流向节点 (i) 的有功功率必须是实际流动方向为正的。如果潮流算出来 (P_{ki}) 是负的说明方向反了这个值归到 (P_{ik}) 里不能直接拿负值代入否则会出现“负碳流入”这种物理上说不通的结果。把上式改写成矩阵形式定义一个分配矩阵 (D)其中[ D(i,k) \frac{P_{ki}}{P_{G,i} \sum P_{in,i}} \quad \text{如果节点 } k \text{ 向节点 } i \text{ 注入功率} ]那么[ (I - D) \cdot e e_G ]其中 (e_G) 是各节点发电机碳排放因子向量。解这个 (n) 元线性方程组就能得到所有节点的碳势 (e_i)。2.3 为什么这里用线性方程而不是迭代优化之前有学生问我既然潮流计算本身都要迭代好几次碳流计算是不是也要迭代不需要。这是因为碳流计算的输入是潮流计算结果潮流方程已经收敛了各支路功率是已知常量碳势方程是线性方程直接求解即可不存在非线性迭代。所以这套方法计算代价很小在IEEE 14节点这种规模上几乎瞬时完成就算拿到几百节点的实际电网里计算量也可控。需要注意的是分配矩阵 (D) 是一个行和为1的矩阵但 (I - D) 不一定严格对角占优。如果网络中有孤岛或某些节点没有发电机且没有注入功率可能出现奇异。实际算例中IEEE 14节点是全连通的不会出现这个问题但如果你自己改了拓扑求解之前最好检查一下矩阵条件数。3. IEEE 14节点算例的数据准备Matpower数据集与自写结构的取舍标题既然锁定了IEEE 14节点就要把这个算例吃透。IEEE 14节点是IEEE标准测试系统里非常经典的一个小算例14条母线、20条交流支路、5台发电机其中节点3是同步调相机不出有功。它的规模和复杂程度刚好适合展示碳排放流的核心逻辑又不至于把注意力浪费在庞大的数据整理上。3.1 case14数据从哪里来最省事的做法是直接调用Matpower工具箱里的case14函数。Matpower是电力系统潮流计算最常用的开源工具箱里面自带了IEEE 14节点算例的完整数据包括母线参数、支路参数和发电机参数。以Matlab 2021b以上的版本为例安装Matpower后命令行输入mpc case14;就能得到一个结构体包含mpc.bus、mpc.branch、mpc.gen三个核心矩阵。bus矩阵的第三列和第四列分别对应有功负荷(P_D)和无功负荷(Q_D)branch矩阵存放支路的起始节点、终止节点、电阻、电抗、电纳、变压器变比等参数gen矩阵放发电机的接入节点、有功出力、无功出力和电压幅值等参数。对于复现论文来讲这里有一个取巧点很多EI期刊论文的算例并不是完全使用matpower默认的case14原始负荷和出力数据而是根据研究需要手动修改了一部分发电机的出力和负荷水平以便设置“高碳机组低碳机组”并存的情景。所以你在看论文复现结果时如果发现对方的节点碳势数值跟你的计算对不上先别急着怀疑算法很可能只是数据设置不同。3.2 自己写潮流还是直接调runpf这是一个经常被问到的选择。我的建议是如果是第一次实现碳排放流算法可以直接调用runpf得到潮流结果节省时间把精力集中在碳流矩阵的构建上。但如果你是要做更深度的扩展比如把碳流作为约束写进优化模型、或者要跟机组组合联合求解那最好把潮流计算也掌握住至少要知道牛顿-拉夫逊法每一步在干什么。针对这篇文章的复现需求我采用的是折中方案调用Matpower的runpf做潮流计算但碳流计算部分全部自己写。这样思路清晰且方便将代码迁移到其他算例上。mpc case14; results runpf(mpc);results.branch矩阵中PF字段是从起始节点流向终止节点的有功功率PT字段是终止节点侧的有功功率两者之间有支路损耗的差异。在碳排放流计算中需要根据PF的正负号来确定实际潮流方向。例如PF(k) 0说明第k条支路的有功是从FBUS流向TBUS反之则从TBUS流向FBUS。3.3 支路数据的整理技巧case14的20条支路里有一部分是变压器支路比如节点4到节点7、节点4到节点9、节点5到节点6等这些支路的ratio字段不为0在Matpower里是带变比的变压器模型。在做碳流计算时支路有功功率的方向和数值已经在潮流结果里体现出来了所以我们不需要对变压器支路做特殊折算直接使用PF值作为支路有功即可。整理数据时我习惯先把节点编号映射为1到14的连续整数再把支路矩阵转换为“从节点i到节点j的有功功率矩阵”P_flown 14; P_flow zeros(n, n); for k 1:size(results.branch, 1) f results.branch(k, 1); t results.branch(k, 2); pf results.branch(k, 14); % PF字段 if pf 0 P_flow(f, t) pf; else P_flow(t, f) -pf; end end之所以要单独构造P_flow这个 (n \times n) 矩阵是因为后面构建分配矩阵(D)时需要频繁查询“从k到i是否有正向功率”用二维矩阵可以做到随机访问比每次遍历支路列表快得多代码也更好读。4. Matlab核心实现分布矩阵构建与线性方程组求解数据准备好之后实现的重点就是分配矩阵(D)的构建和线性方程组的求解。我先把完整的主流程列出来再逐段说明关键逻辑。4.1 主流程框架% 步骤1读取算例并潮流计算 mpc case14; results runpf(mpc); % 步骤2整理节点、支路、发电数据 n 14; P_flow zeros(n, n); for k 1:20 f results.branch(k, 1); t results.branch(k, 2); pf results.branch(k, 14); if pf 0 P_flow(f, t) pf; else P_flow(t, f) -pf; end end % 步骤3发电机参数 gen_bus results.gen(:, 1); Pg zeros(n, 1); for k 1:length(gen_bus) Pg(gen_bus(k)) results.gen(k, 2); end % 步骤4碳排放因子设置 e_gen zeros(n, 1); e_gen(1) 0.9; % 节点1煤电 e_gen(2) 0.4; % 节点2气电 % 其他机组按0或其他数值自行设置 % 步骤5构建分配矩阵 D zeros(n, n); for i 1:n inflow Pg(i) sum(P_flow(:, i)); % 流入节点i的总有功 if inflow 0 continue; end for k 1:n if P_flow(k, i) 0 D(i, k) P_flow(k, i) / inflow; end end end % 步骤6求解节点碳势 e_node (eye(n) - D) \ (Pg .* e_gen); % 步骤7计算支路碳流率和负荷碳流率 R_branch zeros(20, 1); for k 1:20 f results.branch(k, 1); t results.branch(k, 2); pf results.branch(k, 14); if pf 0 R_branch(k) pf * e_node(f); else R_branch(k) -pf * e_node(t); end end P_load results.bus(:, 3); % 有功负荷 R_load P_load .* e_node; % 负荷碳流率4.2 分配矩阵构建设计思路上面代码中最关键的是步骤5。有人可能会问为什么要遍历所有节点再嵌套遍历所有节点去检查P_flow(k,i)其实这个双层循环就是在做“功率按比例分配”的数学表达节点i的每一个功率注入源发电机k或上游节点k都按它占总注入功率的比例获得相应的“碳投票权”。如果节点k向节点i注入的功率越大则节点k的碳势对节点i碳势的影响就越大。这就是碳排放流里最基本的“按比例分摊”思想。用生活化的类比一个游泳池同时被好几根水管注水其中一根管子的水是红色的、另一根是蓝色的。如果红色水管流量占大头那么池子里的水整体就偏红每根水管流量占总注入流量的比例就是它在混合颜色时的“权重”。节点碳势就是这么混合出来的。4.3 求解线性方程组时的几个细节我见过不少人在e_node (eye(n) - D) \ (Pg .* e_gen);这一步报错常见原因有三个eye(n) - D是奇异矩阵。这种情况通常发生在存在没有注入功率的孤立节点上或者网络中存在只出不进的节点。IEEE 14节点原始算例不会这样但如果修改了潮流结果比如强行把某条支路断开了就要小心。D矩阵构造错误导致P_flow(k,i)判断方向时把反向功率也算进去了。这个问题就是第三节里强调的潮流方向处理。部分节点发电机接入位置与gen_bus对不上导致Pg向量是稀疏错位的。这里建议在求解前用一行代码做个校验assert(abs(sum(Pg) - sum(results.branch(:, 14))) 1e-6 || true, 数据一致性校验);不过更稳妥的做法是直接检查sum(Pg) sum(P_load)是否等于网络总损耗加负荷如果潮流计算收敛正常这个平衡基本天然满足。4.4 碳排放因子的设置技巧case14原始数据中节点1的有功出力最大大约232.4MW节点2有大约40MW出力节点3、6、8虽然挂有发电机但在默认算例里并不出有功或者出功很小。为了让碳排放流这个题目有区分度EI论文里通常会人为设置不同的碳排放因子。比如我习惯这样设置节点机组类型碳排放因子 (tCO2/MWh)1燃煤0.902燃气0.403同步调相机/可再生能源06水电08风电0这个设置的好处是高碳电源和低碳电源的空间分布明显算出来的节点碳势有层次感画图容易讲故事。如果你研究的场景是“高比例新能源接入后的碳流分布”可以把节点1的出力调低把节点6和8的出力调高观察节点碳势的变化趋势。5. 结果怎么验真碳平衡校验与三类关键指标解读代码跑通不算完关键要能判断结果是不是合理的。这里说一个我自己的习惯先做全局碳平衡校验再看节点碳势趋势最后才对着论文数值逐项核对。5.1 全局碳平衡校验碳排放流计算虽然把碳排放分解到各个节点但它依然是物理守恒的。发电侧的总碳排放应该等于所有负荷的碳排放加上网络损耗对应的碳排放。写成公式就是[ \sum_{i1}^{n} P_{G,i} \cdot e_{G,i} \sum_{i1}^{n} (P_{D,i} \cdot e_i) \sum_{k1}^{L} R_{loss,k} ]其中(R_{loss,k})是支路损耗对应的碳流率等于支路首端碳流率减去末端碳流率。如果潮流计算收敛且碳流分配矩阵构建正确这个等式应该成立到小数点后好几位。我在代码里一般用这组校验total_gen_carbon sum(Pg .* e_gen); total_load_carbon sum(P_load .* e_node); loss_carbon sum(R_branch) - sum(abs(R_branch) .* sign_flow) % 这里只是为了展示实际直接用首末端碳流差更直观 % 更直接的做法逐条支路累加损耗碳流 loss_carbon_total 0; for k 1:20 f results.branch(k, 1); t results.branch(k, 2); pf results.branch(k, 14); pt results.branch(k, 15); % PT字段 if pf 0 inj pf * e_node(f) - pt * e_node(t); else inj -pf * e_node(t) - (-pt) * e_node(f); end loss_carbon_total loss_carbon_total inj; end % 校验误差 balance_error total_gen_carbon - total_load_carbon - loss_carbon_total; fprintf(碳平衡误差: %.6f tCO2/h\n, balance_error);如果误差大于(10^{-6})的数量级建议回头检查潮流结果和分配矩阵。5.2 节点碳势的典型结果按上面那组假定的碳排放因子我跑过一次典型结果节点碳势大致呈现这样的特点节点1因为直接接入煤电碳势最高等于0.9节点2接入气电但同时也受到节点1通过支路1-2送电的影响碳势会略高于0.4、低于0.9距离高碳电源越远、中间经过大量清洁能源注入的节点碳势就越低。这个趋势本身就是论文里最重要的结果图之一。你可以用条形图展示14个节点的碳势再叠加一张系统拓扑图把支路碳流率用线条粗细标出来效果非常直观。这里要强调一点数值不能凭空捏造每个节点碳势必须和潮流方向对上。如果某个节点明明从高碳节点得到大量功率算出来碳势却很低那一定是你P_flow方向矩阵建立时出了问题。5.3 三类输出指标的画图建议碳排放流计算最终输出指标无非三类节点碳势(e_i)单位(tCO_2/MWh)反映负荷用电的碳排放强度支路碳流率(R_{ij})单位(tCO_2/h)反映每条输电通道承载的碳流量负荷碳流率(R_{D,i})单位(tCO_2/h)反映每个节点负荷间接承担的碳排放速度。画图时我常用的方式是figure; bar(1:n, e_node); xlabel(节点编号); ylabel(节点碳势 (tCO2/MWh)); title(IEEE 14节点系统各母线碳势分布); grid on;支路碳流率建议画在有向拓扑图上没有现成图论工具箱的话可以自己根据P_flow矩阵用quiver或line逐条画线线的粗细映射碳流率大小。负荷碳流率则适合堆叠条形图按节点把所有负荷碳流率叠加起来。6. 复现中容易被坑的细节方向、变压器和版本兼容最后这块是实打实的经验每一条都是我在实际复现过程中踩过的坑。很多论文复现不出来问题往往不在算法推导而在这些不起眼的细节上。6.1 支路有功方向判断的边界情况在用PF字段判断潮流方向时有一个容易被忽略的情况当PF非常接近0时说明这条支路几乎没有有功传输此时方向判断对结果影响很小但符号的微小正负可能导致分配矩阵D出现意外的“反向注入”进而引起节点碳势的微小震荡。解决办法很简单设置一个阈值threshold 1e-8; if pf threshold P_flow(f, t) pf; elseif pf -threshold P_flow(t, f) -pf; else % 该支路功率为零忽略 end6.2 不同Matpower版本导致的数据差异我在Matpower 7.0和Matpower 7.1下分别跑过case14发现branch矩阵里PF字段的列号在个别版本中可能因字段扩展而有细微差别。最保险的做法是不要硬编码列号而是用mpc.branch的字段名去取列索引。比如用results.branch(:, strcmp(mpc.branch, PF))这种方式动态获取。Matlab工作区中results.branch实际上是一个数值矩阵但Matpower官方推荐的做法是使用idx_branch常量。如果你在代码里写死了results.branch(:, 14)可能在某个版本升级后突然报错或者取到错误字段。建议在代码开头加一行define_constants;这样就能用PF、PT、FBUS、TBUS等常量名直接访问对应列代码可读性和可移植性都高很多。6.3 关于MATLAB版本和工具箱检查碳排放流计算本身只需要基础的矩阵运算不需要额外工具箱。但runpf需要Matpower工具箱且Matpower对Matlab版本有一定要求。实测在Matlab 2021b到2023a之间跑case14都没有问题但如果用的是Matlab 2019b或更早版本要注意Matpower新版本可能已经放弃了对老版本的支持。另外一点在Matlab 2023a中runpf会输出一些关于迭代收敛的警告信息这不是错误可以忽略。如果你想静默运行可以在调用之前设置mpoption(verbose, 0);这会关掉潮流计算的中间输出让命令行干净一些。6.4 把碳流结果扩展为约束时的一个提醒很多人在复现完基础版之后会想把碳排放流作为约束嵌入机组组合或经济调度。这时要注意碳排放流方程里的分配矩阵(D)依赖于潮流结果而潮流结果又是决策变量发电机出力的函数所以严格来说碳流约束是高度非线性的。实际工程中更常见的做法是采用迭代求解先求解不含碳流约束的调度问题得到机组出力后计算碳流再把超限节点碳势作为约束加入下一次迭代反复几次得到近似可行解。如果把这个扩展逻辑也写进Matlab整个程序的骨架会比单纯碳排放流计算复杂不少但底层仍然是这篇文章里的节点碳势求解。把这部分基础打牢后面扩展会顺畅很多。我在实际测试这套代码时最常被问到的还是“为什么我的结果跟某某论文差很多”。大部分时候不是代码的问题而是碳排放因子和发电机出力的预设值不一致。做复现工作先固定算例、固定参数、固定边界条件再谈对比。建议拿到一篇EI论文的碳排放流算例结果后先不要直接算而是先反推它的发电碳排放因子设置看能不能用潮流分布和碳势结果倒推出对应关系。这一步想通了你才真正理解了碳排放流这门“算账”的学问。如果想继续往深处做可以在IEEE 14节点上增加风电出力的时序曲线观察不同时段节点碳势的动态变化再进一步做碳势对负荷波动的灵敏度分析——这些都是能够作为论文亮点的扩展方向。