最近在做一个双碳方向的电力系统分析项目,需要把碳排放流算法落到具体算例上跑通,目标就是复现EI期刊里的那套方法。折腾了一周,把IEEE 14节点系统上的Matlab实现完整跑通了,这里把整个思路、公式推导、代码实现和踩坑过程写出来。这篇内容适合正在做碳计量、碳排放核算、电力系统低碳规划方向的研究生和工程师参考,也适合想接触电力系统碳排放流计算但对完整流程还不太熟悉的读者。看完你可以直接照着复现一套可运行的结果。
碳排放流这个概念,简单说就是把“碳排放”看成随电能一起流动的虚拟物质。发电厂发出电量,同时也“发出”碳排放;负荷用电,实际上也在“消耗”碳排放。这种视角能把发电机侧的排放责任精确分摊到每个负荷节点,是碳核算从宏观到微观的关键一步,也是“双碳”目标下电网精细化管理的核心算法之一。IEEE 14节点是电力系统分析最经典的标准测试算例,节点规模适中,既能验证算法正确性,又不至于被拓扑复杂度淹没,非常适合作为碳排放流方法的验证平台。
1. 碳排放流的核心思路与公式拆解
1.1 碳排放流到底在算什么
先说个最容易糊涂的问题:碳排放流不是真实存在的物理流,它不能像有功功率、无功功率那样用仪表去测量。它是一种“伴随流”,跟着有功功率的流动路径走。为什么是伴随有功而不是无功?因为碳排放的本质是燃料燃烧产生的,而燃料消耗和发电机有功出力直接相关,无功功率虽然影响网损和电压,却不会直接导致碳排放在电网中的“分摊路径”改变。
这个思想其实可以用生活里的例子来理解。一个小区有两台锅炉供热,一台烧天然气、一台烧煤,热水通过管道送到各家各户。我们想知道每家每户消耗的热量对应了多少煤、多少气,最简单的方法就是看热量怎么流、流量多大,然后把两个锅炉的“污染指标”按流量比例摊到每一户头上。电网里的碳排放流就是这么干的:发电机是“锅炉”,负荷是“住户”,输电线路和变压器是“热水管道”,有功功率的大小和方向决定了“碳”往哪儿流、流多少。
这个思路的核心价值在于“责任可追溯”。传统碳核算只能做到区域或电网层面,属于宏观大盘子;碳排放流算法则可以把排放责任精确到节点、到用户、到一条具体的支路。比如某个工业园区从节点5取电,它消耗的每一度电背后是什么类型的机组在发电、碳排放是多少,都能算得清清楚楚。
1.2 核心公式:节点碳势、支路碳流密度、负荷碳流率
碳排放流计算的本质,是建立一个“功率流-碳流”的映射关系。整个算法体系建立在三个核心概念上:节点碳势、支路碳流密度、负荷碳流率。
节点碳势(nodal carbon potential)是所有进入该节点的碳流总和除以进入该节点的有功功率总和。这个概念可以类比为“混合后平均浓度”:多个上游支路和本地发电机组把不同碳强度的电能注入一个节点,在节点上充分混合,形成一个统一的碳势,所有从该节点流出的电能都带有这个碳势。数学上写成:
e_n = (Σ P_G,g × e_G,g + Σ P_in,l × e_branch,l) / (Σ P_G,g + Σ P_in,l)其中P_G,g是节点上第g台发电机的有功出力,e_G,g是这台发电机的碳排放强度,P_in,l和e_branch,l分别是进入该节点的第l条支路的有功功率和该支路的碳流密度。
支路碳流密度(branch carbon flow density)的物理含义是“单位电量通过这条支路时携带的碳排放量”。根据比例共享原则,一条支路的碳流密度等于其送端节点的碳势。也就是电量从哪个节点流出,就带有哪个节点的碳势特征:
ρ_l = e_i其中e_i是支路送端节点i的碳势。这里有个关键细节:送端指的是有功功率实际流出的那一端,不是线路的“电路图上端”。如果潮流方向是从节点j流向节点i,那么送端就是j,支路碳流密度等于e_j而不是e_i。
负荷碳流率(load carbon flow rate)真正回答“这个负荷产生了多少碳排放”这个问题。负荷功率乘以节点碳势,就得到该负荷每小时消耗电量所对应的碳排放量:
R_L,n = P_L,n × e_n按照这个公式,如果节点碳势是800 kgCO2/MWh,节点负荷是50 MW,那么这个负荷的碳流率就是40000 kgCO2/h,也就是每小时产生40吨二氧化碳。
这里要特别强调一个容易被忽略的问题:计算节点碳势需要“按拓扑顺序”推进,而不能随便选节点先算。原因在于节点碳势依赖于上游支路的碳流密度,而上游支路的碳流密度又依赖于更早节点的碳势,存在明显的因果依赖关系。正确做法是从平衡节点开始,顺着有功潮流的实际流向逐级向下游推进,类似有向图的拓扑排序。如果电网中存在环路,则需要迭代计算,直到所有节点碳势不再变化。
2. IEEE 14节点系统与数据准备
2.1 为什么选IEEE 14节点作为验证平台
IEEE 14节点系统是电力系统领域最著名的标准测试系统之一,从20世纪60年代提出至今,几乎所有电力系统分析工具都内置了它的数据文件。它包含14个节点、20条支路(含输电线路和变压器支路)、5台发电机组、11个负荷节点,规模不大但拓扑特征非常丰富:既有辐射状结构,又有环网结构,还有变压器支路和多机节点。这些特征恰好能覆盖碳排放流算法的主要难点。
选它做验证平台有三个实际好处。第一,标准算例的参数是公开的,任何人都能下载到同一套数据,结果可以互相校验;第二,规模适中,用普通电脑跑Matlab几秒钟就能完成潮流计算和碳流计算,调试迭代效率很高;第三,很多EI期刊论文都选择在这个系统上展示碳排放流算例,复现时方便和已发表文献做数值对照。
相比更大的IEEE 118节点或IEEE 300节点系统,14节点的拓扑足够把算法逻辑说清楚,又不至于让调试过程淹没在数据问题里。对刚接触碳排放流的人来说,这是性价比最高的起点;对有经验的研究者来说,在小系统上先把逻辑跑通、再换到大系统,也是标准的稳妥路线。
2.2 系统参数与碳排放强度设置
IEEE 14节点系统的经典数据包含母线参数、支路参数、发电机参数三部分,Matpower的case14函数直接内置了这些数据。以下几个关键参数在处理碳流时必须心里有数。
发电机分布情况:节点1是平衡节点,节点2、3、6、8是PV节点。节点3、6、8上的发电机在经典参数中通常作为调相机运行,有功出力很小或为0,主要提供无功支撑。这一点对碳流计算影响很大,因为碳排放流只跟随有功功率流动,有功出力为0的发电机实际不对节点碳势产生贡献,在代码里可以直接按0处理。
负荷分布情况:系统总负荷在260 MW左右,主要集中在节点2、3、4和9,其中节点3负荷最大,约94 MW。负荷越大的节点,即使碳势不高,负荷碳流率也可能很高。这是后续结果分析中需要重点关注的节点。
碳排放强度的赋值是最体现研究者主观性的环节。IEEE 14节点本身不提供碳排放参数,需要自己设定。我在复现时的设定如下:
| 节点 | 机组类型 | 有功出力/MW | 碳排放强度/(kgCO2/MWh) |
|---|---|---|---|
| 1 | 燃煤机组 | 232.4 | 800 |
| 2 | 燃气机组 | 40.0 | 400 |
| 3 | 调相机 | 0 | 0 |
| 6 | 调相机 | 0 | 0 |
| 8 | 调相机 | 0 | 0 |
这个设定参考了当前国内电源结构的一般特征:燃煤机组碳排放强度在750-900 kgCO2/MWh区间,燃气机组在350-450 kgCO2/MWh区间。实际复现时,这个参数完全可以根据研究目标调整,比如模拟高比例新能源接入,可以把部分机组碳强度设为0,考察零碳电源对全系统碳势的“稀释”效果。
2.3 环境准备:Matlab与Matpower安装要点
这个项目的实现依赖Matlab和Matpower工具箱。Matpower是电力系统潮流计算最常用的开源工具箱,自带case14数据文件和runpf潮流计算函数,省去了自己写潮流程序的麻烦。
Matpower的安装比想象中简单,但有几个细节值得说。首先是版本适配问题,Matpower对Matlab版本没有特别严格的要求,R2016a到最新的R2024a都能正常使用,但建议用R2020b以上版本,避免一些语法兼容问题。安装时把解压后的文件夹放进Matlab路径即可,通过Set Path把整个文件夹加入搜索路径,然后在命令行输入test_matpower检查是否安装成功。
这里提示一个实际工作中经常遇到的坑:如果电脑上装了多个版本的Matlab,或者Matpower路径下有多个版本,会出现函数冲突。表现是调用runpf时提示找不到函数或者版本错误。解决办法是在路径设置里只保留一个Matpower文件夹,并且用rehash toolboxcache重刷工具箱缓存。
如果不想依赖Matpower,也可以自己写牛顿-拉夫逊潮流,但工作量会大很多。我的建议是:对算法本身感兴趣、目标是跑通碳排放流逻辑的,直接用Matpower,把精力集中在碳流计算层;如果是做深入研究需要改动潮流算法本身,再考虑手写。
3. Matlab实现流程与关键代码
3.1 整体代码架构设计
碳排放流的Matlab实现,我采用的是“先潮流、后碳流”的两段式结构。这样设计的好处是模块职责清晰:潮流部分负责算出有功功率的分布和方向,碳流部分专注处理碳势分配逻辑。两个模块通过一个结构体result传递数据,互不干扰,后续如果想替换潮流算法或修改碳流模型,只需要改动对应模块。
整个程序拆成三个文件:主脚本、碳流计算函数、结果可视化函数。主脚本负责加载数据、设置参数、调用潮流计算和碳流计算、展示结果;碳流计算函数实现节点碳势迭代求解和支路碳流密度分配;可视化函数则把碳势分布和碳流率以图表形式呈现。
主脚本的核心代码如下:
%% run_carbon_flow_ieee14.m % 加载IEEE 14节点标准算例 mpc = loadcase('case14'); % 设置潮流计算选项,verbose=1可以在命令行看到收敛信息 mpopt = mpoption('verbose', 2, 'out.all', 0); % 求解交流潮流 result = runpf(mpc, mpopt); if ~result.success error('潮流计算未收敛,请检查输入数据'); end % 设置发电机的碳排放强度,单位: kgCO2/MWh % 按节点1,2,3,6,8的顺序对应5台发电机组 gen_ems = [800; 400; 0; 0; 0]; % 调用碳排放流计算函数 cf = carbon_flow_calc(result, gen_ems); % 展示核心结果 disp_table(cf.bus_carbon_potential);这个结构任何人拿到都能快速看懂,替换成自己的算例时只需要换loadcase的名称和gen_ems的赋值。
3.2 潮流结果的数据提取与方向判断
潮流计算完成后,result.branch矩阵包含了所有支路的潮流结果,其中第14列是支路首端有功PF,第16列是支路末端有功PT。这里需要特别小心方向问题:PF和PT都是带符号的,正负取决于Matpower内部的节点编号方向约定。
碳流计算中,支路方向的判断直接影响送端节点碳势的取值。我处理的方法是把每条支路的潮流方向统一为“从送端到受端”的表达:
% branch矩阵列含义可查Matpower文档: [F_BUS T_BUS ... PF QT ...] for k = 1:nl f_bus = branch(k, 1); % 首端节点 t_bus = branch(k, 2); % 末端节点 pf = result.branch(k, 14); % 首端有功 if pf > 0 % 实际潮流从f_bus流向t_bus, 送端是f_bus send_bus(k) = f_bus; recv_bus(k) = t_bus; branch_flow(k) = abs(pf); else % 实际潮流反向, 送端是t_bus send_bus(k) = t_bus; recv_bus(k) = f_bus; branch_flow(k) = abs(result.branch(k, 16)); end end这一步是整个碳流计算中最容易出错的环节。如果方向判断错了,后面所有碳势计算都会跟着错,而且错误不会很明显,因为数值看起来“合理”但跟文献对照时对不上。我在开发过程中专门写过一段校验代码,统计每条支路的PF和PT是否满足功率平衡关系,确保方向处理正确后再计算碳流。
还有一个细节:变压器支路和普通输电线路在这个算法里是无差别处理的,碳排放流只关心有功功率的方向和大小,变压器变比只影响潮流计算,不影响碳流分配逻辑。
3.3 节点碳势的迭代求解与逆流剔除
节点碳势求解是这个程序的核心。由于14节点系统存在环网结构,节点间的碳势依赖关系不是严格的树状结构,直接用拓扑排序可能漏算,所以我采用了迭代法:初始化所有节点碳势为0,反复扫描计算每个节点碳势,直到变化量小于阈值。
function cf = carbon_flow_calc(result, gen_ems) bus = result.bus; branch = result.branch; gen = result.gen; nb = size(bus, 1); nl = size(branch, 1); % 节点本机注入功率和碳势 % 先统计每台发电机对应的节点位置 gen_bus = gen(:, 1); gen_power = zeros(nb, 1); gen_emission = zeros(nb, 1); for g = 1:size(gen, 1) b = gen_bus(g); gen_power(b) = gen_power(b) + gen(g, 2); % 有功出力 if gen_power(b) > 0 gen_emission(b) = gen_emission(b) + gen(g, 2) * gen_ems(g); end end % 支路潮流方向判断(简化版) send_bus = zeros(nl, 1); recv_bus = zeros(nl, 1); branch_flow = zeros(nl, 1); for k = 1:nl f = branch(k, 1); t = branch(k, 2); pf = result.branch(k, 14); if pf >= 0 send_bus(k) = f; recv_bus(k) = t; branch_flow(k) = pf; else send_bus(k) = t; recv_bus(k) = f; branch_flow(k) = -result.branch(k, 16); end end % 迭代求解节点碳势 e_node = zeros(nb, 1); max_iter = 100; tol = 1e-6; for iter = 1:max_iter e_node_old = e_node; for n = 1:nb % 节点n的所有注入支路: recv_bus == n 且潮流方向正常 inflow_power = 0; inflow_carbon = 0; for k = 1:nl if recv_bus(k) == n s = send_bus(k); rho = e_node_old(s); % 支路碳流密度=送端碳势 % 逆流剔除规则 if rho >= e_node_old(n) - tol inflow_power = inflow_power + branch_flow(k); inflow_carbon = inflow_carbon + branch_flow(k) * rho; end end end % 加上本机发电 total_power = gen_power(n) + inflow_power; if total_power > 1e-6 e_node(n) = (gen_emission(n) + inflow_carbon) / total_power; else e_node(n) = 0; end end if max(abs(e_node - e_node_old)) < tol break; end end cf.bus_carbon_potential = e_node; cf.branch_flow_direction = [send_bus, recv_bus]; cf.branch_carbon_density = arrayfun(@(k) e_node(send_bus(k)), (1:nl)'); cf.branch_carbon_flow = cf.branch_carbon_density .* branch_flow; cf.load_carbon_flow = bus(:, 3) .* e_node; % bus第3列是有功负荷 end这段代码里有几个值得琢磨的细节。
第一,逆流剔除规则。如果一条支路送端碳势低于受端碳势,也就是ρ_l < e_n,意味着高碳势节点在“吸收”低碳支路流入的碳流,这在物理上说不通。实际处理时我把这类支路排除在碳势计算的注入项之外,否则会出现高碳节点被“稀释”的反常现象。这个规则在论文里叫“逆流剔除”,我在复现时发现它直接影响环网结构的收敛性和合理性。
第二,迭代初值的选择。全部设0是最省事的做法,但如果在某些特殊拓扑下迭代次数会比较多。实测IEEE 14节点系统一般迭代10-20次就能收敛到1e-6的精度,性能完全不是问题。对于更大的系统,可以先用拓扑排序按“从平衡节点往外扩展”的方式给一个更好的初值,能明显加速。
第三,平衡节点1的碳势是“源头”。在迭代过程中,平衡节点的发电碳流直接决定了整个系统碳势的基准水平。如果平衡节点的碳势设置错误,后面所有节点的碳势都会按比例偏移,但“相对高低关系”是对的。这也是为什么碳流结果分析时既要看绝对数值,也要看节点间的相对差异。
3.4 结果输出与可视化
算完碳流之后,结果展示也很关键。我习惯用三张图和一个控制台表格来呈现:节点碳势柱状图、负荷碳流率柱状图、IEEE 14节点拓扑碳势分布着色图。
节点碳势柱状图能直观看出哪些节点“电比较脏”,哪些节点“电比较绿”。负荷碳流率柱状图则反映碳排放责任的分布,碳势高和负荷大的节点都会“名列前茅”。拓扑着色图最直观,把所有节点按碳势高低从红到绿着色,一眼就能看出碳势从电源中心向负荷末端衰减的空间分布特征。
%% 节点碳势可视化 figure; bar(cf.bus_carbon_potential); xlabel('节点编号'); ylabel('节点碳势 (kgCO2/MWh)'); title('IEEE 14节点系统节点碳势分布'); grid on;控制台表格输出则按节点编号、碳势、负荷、负荷碳流率的顺序排列,方便把结果粘贴到论文或实验报告中。
4. 结果验证与对比分析
4.1 碳势分布与负荷碳流率结果解读
按前面的参数设定跑完程序,得到的节点碳势分布有一个非常清晰的特征:以节点1为中心向外辐射衰减,经过多级功率分配后碳势逐级下降。节点1碳势等于燃煤机组的800 kgCO2/MWh,因为平衡节点直接由燃煤机组注入功率;节点2有燃气机组注入,碳势被“稀释”到500到600左右的水平;离平衡节点更远、经过变压和长距离传输的节点,碳势进一步降低。
这里有一个反直觉的现象值得解释:碳势低的节点不一定负荷碳流率低。节点3的碳势在系统中不算最高,但由于它的负荷高达94 MW,其负荷碳流率反而是全系统最大的。这个现象说明碳流分析的结论要看两个维度:碳势衡量“电的干净程度”,碳流率衡量“碳排放责任大小”。做碳减排决策时,两者要结合看。
对照已发表的EI论文,IEEE 14节点碳势分布的总体规律是一致的:碳势从注入源向外递减,支路碳流密度等于送端碳势,负荷碳流率与节点负荷和碳势的乘积呈正比。如果复现结果与文献数值有差异,优先检查两点:发电机碳强度设置是否一致、潮流运行方式是否一致(比如平衡节点出力可能因负荷模型不同而不同)。
4.2 碳流守恒校验
判断碳流计算是否正确,最有力的依据是“碳流守恒”。这个守恒关系说的是:全系统所有发电机组产生的碳排放总量,等于所有负荷消费的碳排放总量加上所有支路网损对应的碳排放量。写成公式就是:
Σ P_G,g × e_G,g = Σ P_L,n × e_n + Σ P_loss,l × ρ_l如果这个等式成立,说明碳流分配过程既没有“凭空产生碳”,也没有“凭空消灭碳”,计算逻辑是自洽的。这是我在复现过程中最依赖的验证手段。
在代码里,可以用一句简单的命令完成校验:
total_gen_carbon = sum(gen(:,2) .* gen_ems); % 发电总碳流 total_load_carbon = sum(bus(:,3) .* e_node); % 负荷总碳流 total_loss_carbon = sum(cf.branch_carbon_density .* abs(result.branch(:,14) - result.branch(:,16)));在我的运行结果中,发电总碳流、负荷总碳流加网损碳流三者的误差在0.1%以内,主要误差来源是潮流计算的数值精度和支路功率损耗的近似处理。如果这个误差明显偏大,说明碳流分配逻辑有bug,要回头检查支路方向或逆流剔除规则。
4.3 与EI原始文献的数值对比策略
“完美复现”的检验标准不是数值一字不差,而是在合理误差范围内复现原文献的核心规律和关键数值。我通常采用以下几个层次的对比策略。
第一层,比对节点碳势的相对大小关系。即使文献设置的发电机碳强度和我的不完全相同,各节点碳势从高到低的“排序”应该基本一致。如果排序乱了,说明碳流分配逻辑有问题。
第二层,比对负荷碳流率的分布特征。重点看负荷最大的几个节点是否如预期承担了最大的碳流率,以及碳流率的占比和负荷占比的差异是否与文献描述一致。
第三层,比对守恒关系。文献中如果给了发电总碳流和负荷总碳流的数据,可以直接核对绝对值是否吻合。如果文献只给相对值,就要换算后再比。
复现时不要盲目追求“一模一样”。不同文献对发电机碳排放强度的设置不同,潮流运行点也可能因版本差异而略有不同,这些都会导致具体的碳势数值有出入。关键是抓住“规律正确”这个核心,再根据文献的具体参数调整设置后做精确对比。
5. 常见问题与排坑实录
5.1 潮流不收敛怎么办
碳排放流计算的前提是潮流正确收敛,如果runpf报错,后面全白搭。我遇到的最常见情况是Matpower版本更新后,case14数据格式或默认选项发生变化,导致原本收敛的算例报错。
排查思路:先看报错信息是在数据检查阶段还是迭代求解阶段。如果在数据检查阶段,多半是mpc结构里某些字段格式不对,或者当前Matpower版本要求额外的字段。如果在迭代求解阶段,可以调整mpoption里的迭代次数上限和收敛精度:
mpopt = mpoption('verbose', 2, 'out.all', 0, 'pf.max_it', 50, 'pf.tol', 1e-8);还有一个容易踩的坑:如果修改了case14数据(比如改了负荷或发电机出力),可能导致潮流无解。这时候要检查修改后的总负荷与总发电是否匹配,平衡节点出力是否有足够调节空间。
5.2 平衡节点的碳势设定和特殊处理
平衡节点的碳势设置是整个计算里最容易困惑的地方。因为平衡节点的有功出力是潮流计算自动算出来的,不是预先给定的,所以它的发电碳流和碳势是“结果”而不是“输入”。
调用runpf后,平衡节点的有功出力在result.gen矩阵中已经包含了最终值。比如经典case14的平衡节点出力大约是232.4 MW,这个值会随着负荷水平变化。在你写碳流计算函数时,用的是result.gen里的值,而不是自己预先给定的值,否则会造成碳流不守恒。
还有一种特殊情况:如果某个节点上有多台发电机,比如研究场景中某个节点同时接了燃煤和燃气机组,那么这个节点的发电碳流应该是各台机组碳流之和,节点等效发电碳势等于加权平均值。代码中我已经用累加方式处理了这种情况。
5.3 环网结构导致的计算顺序问题
IEEE 14节点系统存在环网,节点之间的碳势依赖关系不是严格的上下游关系。如果完全按节点编号顺序计算,可能出现“下游节点已经算完了,但上游节点碳势还没更新”的问题。
解决办法有两种:迭代法和拓扑排序法。迭代法简单稳健,代码写起来快,适合14节点这种小规模系统;拓扑排序法需要先对潮流方向做图分析,把节点排成严格的有向无环图顺序,然后再单遍计算,效率更高,适合上百节点的大系统。
如果你的系统规模在14节点左右,我的建议是直接上迭代法,代码简洁,不容易出错。等以后换到IEEE 118节点再做性能优化也不迟。
5.4 网损碳流怎么处理
网损对应的碳排放经常被人忽略。每条支路上的功率损耗也是由发电侧供给的,也会产生碳排放。在碳流守恒校验里,如果不把网损碳流算进去,等式一定不成立。
网损碳流的算法很简单:支路首端功率减去末端功率,差值就是网损;再乘以该支路的碳流密度,就是网损碳流。在代码里用abs(result.branch(:,14) - result.branch(:,16))就能算出来,注意要用绝对值,因为首末端有功功率的符号定义在反向潮流时会变。
5.5 与Matlab环境相关的几个实际问题
这个项目对Matlab环境并不挑剔,但有几个实际问题值得提醒。
第一,如果电脑上同时装了Matlab和Octave,不要用Octave运行Matpower,Matpower的某些函数依赖Matlab特有的工具箱,虽然基础版本能在Octave上跑,但结果可能不稳定。
第二,如果想把结果导出为高清图片用于论文,用exportgraphics(gcf, 'output.png', 'Resolution', 300)而不是老的print命令,清晰度高而且不会出现字体或者线条畸变。
第三,如果系统是中文字体缺失导致的绘图乱码,可以在绘图前加一句set(0, 'DefaultAxesFontName', 'Times New Roman'),曲线和标注就不会出现方框。
最后再分享一个我在复现过程中总结的小技巧:调试碳流程序时,不要一上来就跑完整的14节点系统。先构造一个3节点或者4节点的极小系统,手动把潮流结果和碳流结果算出来,然后用代码跑同一组数据对比。这样能快速定位问题是在方向判断、节点碳势公式还是逆流剔除规则上。等小系统完全对上了,再换回IEEE 14节点,基本就不会有原则性错误了。我靠这个办法省下了至少两天的调试时间,效率提升非常明显。