这几年做双碳相关的电力系统分析有一个问题绕不开发电侧的总排放算得清但电送到用户侧之后“这一度电到底对应多少碳排放”就很难说清了。传统做法是全网平均一个值比如电网平均每度电0.5公斤二氧化碳但这忽略了不同时段、不同区域、不同电源结构带来的巨大差异。这次做的“基于IEEE33节点的节点碳势计算与可视化”项目就是把碳排放流理论落到一个标准配电网算例上精确算到每一个节点的碳势水平并可视化展示出来。对做配电网、微电网、双碳核算和电力市场方向的同学来说这套流程可以直接改数据复用也可以作为理解碳排放流理论的完整实操案例。项目本身不复杂但链路很长先要对IEEE33节点系统做潮流计算拿到支路有功功率分布再基于碳排放流理论把电源侧的碳排放强度“顺着潮流方向”分摊到每个节点最后用Python把33个节点的碳势画到拓扑图上。下面我把整个思路、算法、踩坑和代码细节都拆开讲。1. 项目思路与方案选型为什么是IEEE33 碳排放流1.1 碳排放流理论到底解决了什么问题先花点篇幅讲清楚碳排放流是个什么东西。电力系统的碳排放核算过去主要有两个层次一是发电侧的直接排放就是电厂烧了多少煤、排了多少碳这个有在线监测和盘查数据相对好算二是用户侧的间接排放也就是用户用电所对应的发电侧排放这个难就难在电是“混着走的”——一个节点可能同时接受来自多个电源的电力你没法在物理上区分某一度电来自哪台机组。碳排放流理论的核心做法是把碳排放看作一种伴随有功功率流动的“虚拟流”。有功功率从发电机流向负荷碳排放也跟着流过去。在某个节点上所有流入的碳排放量之和除以流入的总有功功率就得到这个节点的“碳势”单位是kgCO2/kWh。这个碳势就是该节点单位用电量所对应的碳排放水平。节点的下游负荷如果消耗1度电就意味着承担了“节点碳势×1kWh”的碳排放责任。这个思路的意义在于它把宏观的电力碳排放核算从“区域平均”推到了“节点级别”。在同一个配电网里靠近光伏节点的负荷碳势可能只有0.2而靠近火电供电入口的节点碳势可能高达0.9。有了节点碳势就能识别高碳负荷、评估分布式电源减排效果、甚至支撑碳电耦合的价格机制。这次项目选择这套理论除了学术热度高更关键的是它计算过程清晰、验证方便特别适合拿标准算例练手。1.2 为什么选IEEE33节点系统做载体IEEE33节点是配电网研究里最常用的标准算例学术论文里出现的频率极高跟电路课里的“经典电路”一个地位。它由33个节点、32条支路组成拓扑是典型的辐射状配电网根节点接上级电网其余节点带负荷。原始参数里还包含5条联络开关支路平时打开故障时合上实现转供可以用来研究配电网重构。选它有几个现实原因。第一参数公开且标准全网都能找到原始数据做出来的结果别人可以复现、可以对比这是研究工作最看重的事情。第二规模适中——33个节点不算大调试方便几秒钟就能跑完潮流但又比那些三节点五节点的教学算例复杂得多能够覆盖配电网计算的主要细节比如分支线、末端电压偏低、非均匀负荷分布等。第三IEEE33是辐射状网络碳排放流计算会退化成比较简单的“沿支路逐级向下游传递”的模式便于用清晰的代码实现也便于验证。1.3 整体技术路线这次项目整体分了四步走。第一步是整理IEEE33节点数据包括支路阻抗、节点负荷、根节点电压并把有名值转换成标幺值第二步是潮流计算求出各支路首端有功功率和各节点电压这一步是所有碳势计算的基础第三步是碳排放流计算按拓扑顺序从电源节点开始逐点算出每个节点的碳势和每条支路的碳流密度、碳流率第四步是结果可视化把节点碳势用颜色映射到拓扑图上输出图表和结果表。这个技术路线本身没有什么新奇的地方但它是目前把碳势计算落地的比较直接的一条路径。难点不在单点技术上而在各环节的衔接——比如潮流结果怎么传给碳势计算、支路方向怎么统一、多电源场景怎么构造。这些我都会在后面的实操部分展开讲。2. 核心算法拆解节点碳势的计算逻辑2.1 第一步潮流计算拿到功率分布没有准确的潮流分布碳势计算就是空中楼阁。碳排放流是跟随有功功率走的所以碳势计算的第一步永远是潮流计算而且必须得有稳态潮流结果里有功功率的方向和数值。IEEE33是辐射状配电网工程上最常用的方法是前推回代法。这个方法原理不复杂先假设所有节点电压都是额定值然后从网络末端向根节点逐段回推根据节点负荷和支路阻抗计算各支路流过的功率再根据支路功率和阻抗从根节点向末端逐段推算新的节点电压反复迭代直到两次迭代的电压差小于收敛判据。用前推回代法的好处是对辐射状网络几乎无条件收敛不需要形成雅可比矩阵代码写起来也直观。当然也可以用牛顿-拉夫逊法但配电网R/X比值大牛顿法有时需要更小心的初值处理前推回代反而更稳。在IEEE33这种纯辐射状网络上我不太建议用牛拉法除非你后面要做全网状态估计没必要在碳势项目里给自己加难度。收敛判据一般选节点电压幅值偏差小于1e-6 p.u.功率基准取10 MVA。收敛后拿到每个节点的电压幅值和相角同时把每条支路首端的有功功率保存下来这才是后面碳排放流计算的“输入燃料”。2.2 第二步支路碳流密度与节点碳势拿到潮流结果后碳排放流的计算就清晰了。先明确两个概念支路碳流密度和节点碳势。支路碳流密度指流过某条支路的单位有功功率所对应的碳排放量单位也是kgCO2/kWh节点碳势指流入该节点的单位有功功率所对应的碳排放量两者在概念上是对应的。辐射状网络中节点碳势的计算顺序是确定的从根节点开始沿功率流动方向逐级向下游推进。对于节点i如果它从上一级支路流入的有功功率是P_in该支路的碳流密度是ρ_in节点i本地接入的电源有功功率是P_G电源碳排放强度是e_G那么节点i的碳势ρ_i为ρ_i (P_in × ρ_in P_G × e_G) / (P_in P_G)也就是说节点碳势等于“流入碳排放总量”除以“流入有功总量”。纯负荷节点如果没有本地电源P_G就是0节点碳势直接等于上游支路的碳流密度。从该节点流向下游各支路的碳流密度都等于该节点的碳势因为节点对电能混合后流向各个方向的电都具备同样的碳强度。这里有个关键点节点负荷本身不参与碳势计算的分母——负荷是“流出”功率不是“流入”功率。它只是带走了一个等于节点碳势的碳排放强度。换句话说节点碳势描述的是节点处混合电能的属性而不是某一笔负荷的属性。2.3 三步理解碳势计算一个三节点手算示例为了把公式讲明白我拿一个三节点放射链做手算推演。节点1接上级电网碳排放强度0.9 kgCO2/kWh节点2带负荷60kW节点3带负荷40kW同时接入一个20kW的光伏电站光伏碳排放强度按0计算。先算潮流。节点1流向节点2的功率是100kW节点2流向节点3的功率是40kW因为节点3本地光伏出20kW但负荷有40kW还差20kW需要从上游取加上节点2负荷60kW所以节点1到2要送100kW。这个数字不一定精确满足电压损耗做概念示例够用了。然后算碳势。节点1是电源入口碳势就是0.9。支路1-2的碳流密度等于节点1的碳势也就是0.9。节点2没有本地电源碳势就是0.9。节点2流向节点3的功率40kW碳流密度0.9。节点3有本地光伏20kW碳强度0所以节点3的碳势为ρ_3 (20 × 0.9 20 × 0) / (20 20) 0.45 kgCO2/kWh可以看到由于光伏接入节点3的碳势从上游的0.9降到了0.45相当于这个节点上的负荷用电有一半是绿电。这就是节点碳势计算最直观的解读本地接入的低碳电源能够降低该节点及其下游所有节点的碳势。手算完这个例子再去看IEEE33代码逻辑就顺畅很多。3. 实操过程从数据准备到可视化落地3.1 数据准备IEEE33标准参数与标幺化IEEE33节点系统的基准数据需要先交代清楚。我采用的版本是基准电压12.66 kV基准功率10 MVA根节点电压1.0 p.u.总负荷有功3715 kW、无功2300 kvar。节点编号采用0基也就是根节点为节点0其余负荷节点为1到32。还有一种文献把根节点编为1、负荷节点编为2到33引用数据时要注意辨别不然代码里对不上号。支路参数是线路模型的核心。举前几段为例支路0-1的电阻0.0922Ω、电抗0.0470Ω支路1-2的电阻0.4930Ω、电抗0.2511Ω支路2-3的电阻0.3660Ω、电抗0.1864Ω。完整数据网上都能查到建议在代码里以数组或DataFrame形式维护别手敲容易错。标幺化这一步很多初学者会忽略但它是保证潮流计算不跑偏的前提。IEEE33的基准阻抗是Zbase UB² / SB 12.66² / 10 16.0276 Ω以支路0-1为例R* 0.0922 / 16.0276 ≈ 0.00575 X* 0.0470 / 16.0276 ≈ 0.00293负荷标幺化同理10 MW基准下总有功负荷为3.715 p.u.、无功负荷为2.3 p.u.。我在代码里直接用有名值计算最后再标幺也可以直接在数据表里提前标幺好看个人习惯。关键是潮流计算内部必须统一用标幺值或统一用有名值混用一定会出问题。3.2 潮流计算实现细节前推回代法的代码实现比较直接核心就两个方向相反的遍历过程。我给出一个简化版的Python实现思路实际使用时要封装成函数。import numpy as np # IEEE33节点: 支路数据 [首端节点, 末端节点, R(pu), X(pu)] # 负荷数据: [节点号, P(pu), Q(pu)] # 节点0为根节点, 电压恒定1.0pu def forward_backward_sweep(branch_data, load_data, max_iter100, tol1e-6): n 33 V np.ones(n, dtypecomplex) # 初始电压 S_load np.zeros(n, dtypecomplex) for node, p, q in load_data: S_load[node] p 1j*q # 用深度优先/广度优先确定节点层级顺序 # parent数组记录每个节点的父节点, children列表记录下游节点 for _ in range(max_iter): V_old V.copy() S_branch np.zeros(len(branch_data), dtypecomplex) # 前推: 从末端向根节点计算各支路功率 for node in reversed(order_from_root): # 本节点流向下游的总功率 本节点负荷 所有下游支路功率之和 S_total S_load[node] for child in children[node]: S_total S_branch[child] if parent[node] is not None: # 支路功率 本节点总功率 支路损耗 idx branch_index[(parent[node], node)] S_branch[idx] S_total (abs(S_total) / abs(V[node]))**2 * branch_data[idx].Z # 回代: 从根节点向末端更新电压 for node in order_from_root: if parent[node] is not None: idx branch_index[(parent[node], node)] S_end S_branch[idx] # 电压降落 ΔV (S/V)^* × Z V[node] V[parent[node]] - (np.conj(S_end) / np.conj(V[parent[node]])) * branch_data[idx].Z if np.max(np.abs(np.abs(V) - np.abs(V_old))) tol: break return V, S_branch代码里有两个细节值得注意。一是前推时计算支路电流要用节点电压模值的平方而不是恒用1.0不然精度会变差二是回代时支路功率要用前推得到的该支路末端功率而不是首端功率否则会丢支路损耗。这个版本的算法在IEEE33上通常5到10次迭代就能收敛。跑完潮流后从S_branch里取每条支路的有功分量P_branch这就是碳势计算所需的支路有功功率。记得要记录方向由于辐射状网络功率基本从根向负荷流动方向是明确的如果某个节点接了大的分布式电源可能出现倒送后文单独讲处理办法。3.3 碳势计算代码实现碳势计算的输入有三个支路有功功率、支路碳流密度、节点本地电源出力和碳强度。辐射状网络下从根节点沿拓扑顺序做一遍前推即可。下面给出核心逻辑。# 碳势计算: 输入潮流结果和电源参数 # node_gen: {节点号: (有功出力MW, 碳排放强度kgCO2/kWh)} # grid_emission: 上级电网碳排放强度, 取0.9 def carbon_intensity_calculation(branch_data, S_branch, node_gen, grid_emission0.9): n 33 node_carbon np.zeros(n) # 节点碳势 branch_density {} # 支路碳流密度 # 根节点碳势 上级电网碳势 node_carbon[0] grid_emission # 按拓扑顺序从根节点向末端推进 for node in order_from_root: # 流向下游所有支路的碳流密度 当前节点碳势 for child in children[node]: idx branch_index[(node, child)] branch_density[idx] node_carbon[node] # 下游节点碳势 本节点碳势默认无本地电源 node_carbon[child] node_carbon[node] # 若下游节点有本地电源, 做混合计算 if child in node_gen: p_gen, e_gen node_gen[child] p_in abs(S_branch[idx].real) # 从上游流入的有功 p_total p_in p_gen node_carbon[child] (p_in * node_carbon[node] p_gen * e_gen) / p_total return node_carbon, branch_density这段代码我做了简化假设每个节点最多只有一条上游支路这是辐射状网络的特征。多个分布式电源混合时只要在对应节点把P_G和e_G换成多台机组加权汇总即可。还有一点支路有功功率方向如果与拓扑方向相反abs取绝对值会掩盖反向潮流的问题更严谨的写法是判断S_branch[idx].real的符号如果为负说明功率从子节点流向父节点此时碳流方向也要反过来计算逻辑要相应调整。这个坑我在后面常见问题里详谈。计算结果输出为每个节点的碳势值。节点0碳势等于上级电网数值其余节点会因为有光伏、燃气机组的接入而出现差异化分布这也是后端可视化的数据基础。3.4 可视化让碳势分布看得见可视化部分我用了networkx做拓扑、matplotlib做颜色映射。networkx读取IEEE33的边列表非常方便节点颜色根据碳势映射色带我选的是RdYlGn_r也就是红色代表高碳、绿色代表低碳一眼就能看出哪些分支碳排放高。import networkx as nx import matplotlib.pyplot as plt from matplotlib.colors import Normalize from matplotlib import cm G nx.Graph() for f, t in branch_endpoints: G.add_edge(f, t) pos nx.spring_layout(G, seed42, k1.2) node_colors [node_carbon[node] for node in G.nodes()] norm Normalize(vminmin(node_colors), vmaxmax(node_colors)) cmap cm.get_cmap(RdYlGn_r) fig, ax plt.subplots(figsize(12, 8)) sm cm.ScalarMappable(cmapcmap, normnorm) sm.set_array([]) nx.draw_networkx_nodes( G, pos, node_colornode_colors, cmapcmap, node_size[300 30 * abs(load_p[node]) for node in G.nodes()], axax ) nx.draw_networkx_edges(G, pos, alpha0.5, axax) nx.draw_networkx_labels(G, pos, font_size8, axax) plt.colorbar(sm, labelNode Carbon Intensity (kgCO2/kWh)) plt.title(IEEE33 Nodal Carbon Intensity Visualization) plt.axis(off) plt.show()node_size这里我做了个处理节点大小随负荷有功增加这样既能看碳势高低又能看负荷轻重。如果要把碳流率也画出来可以在边上标注或调整边的粗细数值取P_branch×branch_density单位为kgCO2/h。实际显示上spring_layout随机性比较大每次运行图的位置都会变。输出到报告里最好给个固定的seed比如seed42保证图和论文文字对应得上。图出来后建议再用颜色条把数值范围标注清楚并且单独输出一张CSV表记录每个节点的碳势方便后续分析和论文引用。4. 常见问题与排查技巧实录4.1 潮流不收敛多半是参数和初值的问题前推回代法虽然稳健但不是不会卡住。我在调试过程中遇到几次“迭代500次还不收敛”的情况排查下来主要是两类原因。一类是支路数据出错尤其是联络开关支路没有正确断开。IEEE33原始数据里包含5条联络支路它们是常开开关在默认潮流计算时必须剔除否则网络变成环网前推回代法就推不回去了。网络变成环网后前推回代法没有可用的树结构结果会乱掉。第二类是单位混用比如节点负荷用了kW但支路阻抗用了标幺值换算后的欧姆两个不同量纲混在一起潮流结果自然发散。建议统一用标幺值并且在程序入口加断言检查总负荷有功应该在3.7 p.u.左右根节点电压应该在1.0 p.u.附近偏差过大说明数据有误。另外一点初值全部取1.0 p.u.是前推回代法的标准做法但如果某个节点接了很大的分布式电源功率倒送严重迭代收敛会变慢。此时可以把分布式电源有功出力先逐步加载从0开始每步加10%分10步逼近目标值也就是“连续化”方法收敛会稳很多。4.2 为什么所有节点碳势都一样这个问题我估计每个第一次做碳势计算的人都会遇到。如果你用IEEE33原始数据、不接任何本地电源直接算碳势结果大概率是33个节点的碳势全部等于上级电网碳排放强度比如都是0.9。这时候不要慌这是正确结果不是bug。原因在前面手算示例里已经说过纯负荷节点没有本地碳源碳势完全继承上游而整个网络只有一个碳源上级电网所有节点自然会等于同一个值。碳排放流理论的意义就在于“差异化”但这个差异化必须建立在多碳源或本地电源的基础上。所以在设计演示场景时至少要在两到三个节点上接入不同碳排放强度的电源比如光伏碳强度0和燃气轮机约0.57 kgCO2/kWh这样节点碳势才会出现有意义的分布。反过来说如果你的代码在单一电源场景下算出了不同的节点碳势那反而说明公式写错了需要回去检查混合计算的分母和分子。我建议先跑这个“全等验证”确认代码正确后再加多电源场景。4.3 支路碳流方向与潮流方向不一致怎么办辐射状网络在大规模光伏接入后局部功率倒送非常常见。光伏大发而本地负荷小时光伏节点可能向上一级支路送出有功这时候支路有功潮流方向与拓扑中的“父节点→子节点”方向相反。碳势计算遇到这种情况不能简单取绝对值。因为碳排放流的方向必须与实际有功潮流方向一致。当支路的有功功率为负方向时节点碳势的计算顺序也要相应改变光伏节点作为“上游”向父节点方向送出的电其碳流密度应该是光伏节点的碳势而不是父节点的碳势。具体处理时有两条路一条是算了潮流之后根据每条支路实际的有功方向重新建立“碳流拓扑”再做BFS另一条是直接构建有向图以有功功率方向为弧方向再按拓扑序计算。第一条路在常规场景下够用但要注意如果电网里同时存在多处倒送拓扑序会变得复杂最好用有向无环图排序来自动处理。这也是IEEE33这种规模算例的极限所在实际工程网络还是要上正式的状态估计和有向图计算框架。4.4 可视化配色与展示建议可视化的坑主要集中在颜色映射和布局稳定性上。颜色映射如果直接使用默认的jet色带视觉上会很强但从阅读体验看红绿语义更符合“碳势高红色、碳势低绿色”的直觉。另外如果节点碳势范围很小比如最大最小只差0.05直接用min到max拉伸色带会发现图上颜色差别不明显。这时候可以在Normalize里人为设定vmin0、vmax1把图谱整体框定在0到1的碳势区间内让分布直观可比较。布局稳定性是另个容易被忽略的问题。networkx的spring_layout受随机数影响不固定seed的话每跑一次节点的位置都不一样同一份代码两次运行结果长得完全不同不利于报告和演示。用seed42这类固定随机种子或者干脆把IEEE33节点的位置按经纬坐标和馈线走向手工摆好效果更专业。我在项目里直接取了一个IEEE33的经典坐标布局把节点坐标存在字典里这样不管怎么运行图都是同一个位置。5. 从计算到应用节点碳势还能怎么用5.1 高碳节点的识别与改造节点碳势计算出来后最直接的用途就是识别配电网中的高碳负荷节点。在多电源供电场景下碳势高的节点意味着它消耗的电能中高碳电源占比更大。在碳势计算结果的模型里负载率高、且由高碳电源供电的节点会呈现明显红色标记。针对这类节点可以优先开展节能改造、错峰用电引导或者在条件允许时建议业主加装分布式光伏。这种识别方式的价值在于“精准”。传统方案是整个配电台区用同一个平均碳排系数核算所有用户高碳的用户被平均掉、低碳的用户被平均高核算结果既不公平也不准确。节点碳势把核算单位细化到了节点甚至分支用户侧的碳账算得更明白。5.2 分布式电源低碳效益的定量评估在配电网里装光伏或风机到底能带来多大的碳减排效果如果只看“装了10MW光伏取代了10MW火电”这是一个相当粗糙的估算。实际上分布式电源接入位置不同对节点碳势的影响范围也不同——装在负荷中心就近抵消高碳负荷碳势下降幅度大装在馈线末端可能只对末端少数节点产生作用。通过节点碳势的前后对比可以算出每个低碳电源接入后下游所有节点碳势的加权下降量再乘以对应电量就是该电源的“碳减排贡献”。这个指标既反映了机组本身的零碳特性也反映了接入位置和就地消纳程度。对规划人员来说这个指标比单纯的“年度减排量”更有决策价值。5.3 动态碳势与实时碳排监测如果把潮流计算从静态扩展到时序仿真就能得到节点碳势随时间的变化曲线。配电网里的分布式光伏日出而作日落而息负荷也有早晚高峰节点碳势一天之内波动会很明显。中午光伏大发时节点碳势可能降到接近0晚高峰时又升上来。基于这个思想可以做实时碳排监测系统采集电网运行数据每隔几分钟算一次节点碳势把结果推送到可视化大屏上。相比传统的月度碳排放核算这种“动态碳势图”能实时反映电网的绿色程度也可以作为碳电联动需求响应的参考信号——用户在碳势高的时候少用电、碳势低的时候多用电本身就是一种“碳需求响应”。这个方向后续可以接时序数据、接数据库、接前端图表框架从单机脚本迭代成一个小型应用。最后再分享一个我实际使用的技巧在跑碳势计算之前一定要先对潮流结果做一遍合理性检查。如果某个节点电压低于0.9 p.u.或高于1.1 p.u.说明网络参数或负荷数据大概率有问题这时候算出来的碳势是不可信的。做工具和做项目一个很大的差别就在于工具能给出结果但好的项目会在结果前面多几道“校验闸门”。这个习惯比这33个节点的计算本身更值得养成。