做电力系统双碳方向的研究最容易卡住的不是算法本身而是找不到一块合适的试验田。真实配电网的拓扑和负荷数据动辄几千个节点拿来验证新思路成本太高不如先用标准算例把流程跑通。IEEE33节点系统就是配电网里最经典的试验田节点碳势计算又是碳排放流理论落地的核心指标再配上一套可视化整个从计算到出图的链路就完整了。这篇文章适合刚接触碳流计算的电力专业学生、做双碳方向的科研人员以及想在配电网优化里引入碳指标的工程师。我会把碳势计算的原理、潮流数据怎么衔接、可视化怎么做、踩过的坑是什么全部摊开来讲。1. 项目思路与整体设计1.1 为什么选IEEE33节点系统IEEE33节点系统来自Baran和Wu提出的33节点配电网络包含33个节点、32条支路馈线出口电压12.66kV总负荷约3715kW加2300kvar。我一直把它当成配电网的hello world原因很简单规模小手工核算每一条支路都可行但又具备真实配电网的全部要素比如多级馈线、分支负荷、联络开关、电压降落问题和分布式电源接入的可能。做碳势计算选它还有一层考量就是结果容易验证。33节点系统的潮流解是可以拿标准数据复算的碳势计算依赖潮流分布潮流对了碳势才有意义。很多论文都用这个算例做碳排放流分析这意味着你算出来的任何一组结果理论上都能和文献数值互相印证这对调试自己的代码帮助非常大。1.2 节点碳势要解决什么问题传统碳排放因子是一个宏观平均概念比如某地区电网每度电排放多少克二氧化碳它很难反映出电网内部不同节点的真实差异。同样是负荷节点有的节点接近大型火电基地接入点买到的电高碳有的节点旁边就是风电、光伏买到的电低碳。节点碳势就是在潮流分布基础上逐节点地算出一个碳强度数值单位通常是gCO2/kWh代表该节点单位用电量对应多少碳排放。这个指标的价值在双碳调度里很明显。高碳节点适合优先响应需求侧调节低碳节点可以多消纳本地新能源规划储能和分布式电源的位置也能据碳势分布来优化。说白了节点碳势把碳排放从发电侧平均概念细化成了电网空间分布概念。1.3 总体技术路线我把整个项目拆成三块潮流计算、碳流计算、可视化展示。第一步必须先跑通交流潮流拿到每条支路的有功功率和流向第二步在潮流结果上构建节点功率通量关系求解每个节点的碳势第三步把节点碳势映射到网络拓扑图上用颜色、大小、柱状图等形式直观展示。值得强调的一点是碳势计算不是独立于潮流计算的另一套系统它是潮流结果的后处理。只要支路有功和发电注入已知碳势方程组的系数就全部确定了。所以整个项目迭代最快的地方在第三步可视化前两步一旦封装成函数基本不需要再动。2. 节点碳势的核心计算方法2.1 碳排放流的基本逻辑碳排放流理论的起点是把发电机产生的碳排放看成附着在电力上流动的物质。电从发电机流向负荷碳也跟着有功功率流向负荷。一条支路输送单位电量所附带的碳排放量叫支路碳流密度一个节点单位用电量对应的碳排放强度就是节点碳势。这里有个重要的近似假设稳态情况下支路碳流密度等于该支路送端节点的碳势。可以这样理解电在输电线路上传输时碳是挂在电量上的不会自己消失也不会凭空增加到了受端节点来的这一路碳流量和电量一除碳势自然就是上游节点的碳势。实际网损会让支路末端的有功功率变小但单位电量的碳强度在经典碳排放流分析里仍近似沿用送端值工程计算精度完全够用。2.2 节点碳势的公式推导对任意节点i所有从其他节点注入它的有功功率加上节点自身的发电注入构成节点i的总注入有功。节点碳势E_i就是这些来流碳势的加权平均[ E_i\frac{\sum_{j \in U_i} P_{j \to i} E_j \sum_{g \in G_i} P_{g,i} E_{g,i}}{\sum_{j \in U_i} P_{j \to i} \sum_{g \in G_i} P_{g,i}} ]其中U_i表示向节点i输送有功的支路集合P_{j \to i}是支路j到i的有功功率E_j是送端节点j的碳势G_i是节点i上的发电机集合P_{g,i}是发电有功E_{g,i}是发电机碳势。这个式子很好理解就是把所有流入节点i的碳流量总和除以总电量得到一个混合浓度。对于辐射状网络碳势可以从送端往末端一层一层推下去只要上游节点碳势已知下游节点碳势就能直接算。但如果系统闭合了联络开关形成环网节点之间可能存在多条来流路径这时候不能直接递推必须把所有节点的方程联立成一个线性方程组。我把公式重新整理成矩阵形式[ A E b ]对角线A_{i,i}为1非对角线A_{i,j}为[ A_{i,j}-\frac{P_{j \to i}}{\sum_{j \in U_i} P_{j \to i} \sum_{g \in G_i} P_{g,i}} ]右边向量b_i为[ b_i\frac{\sum_{g \in G_i} P_{g,i} E_{g,i}}{\sum_{j \in U_i} P_{j \to i} \sum_{g \in G_i} P_{g,i}} ]这样不管是单电源辐射网还是多环网都能用线性求解器直接解出所有节点的碳势。2.3 边界条件怎么给节点碳势计算必须有一个边界条件就好比解微分方程要有初值。在IEEE33节点系统里变电站出口节点0通常作为平衡节点接外电网这个节点的碳势应该设置成外购电的碳排放因子。不同区域的电网碳排放因子差别很大比如有些地区火电占比高碳排放因子可能达到0.8tCO2/MWh以上水电和新能源占比高的地区可能只有0.2到0.4tCO2/MWh。算例里我习惯先取一个方便对照的值比如0.581tCO2/MWh对应某些区域电网的平均水平后面分析分布式电源影响时再改成别的值对比。如果节点接了光伏或风电发电机碳势取0因为新能源运行过程基本不产生直接碳排放。要注意的是如果全系统没有任何节点给定碳势方程组的解会退化成全零这就是很多新手代码跑通了结果却全为0的原因。根本问题不是公式写错而是没给松弛节点设外部碳势。3. 实操过程从潮流计算到可视化3.1 潮流计算环境准备我用的是Python下的pandapower它可以直接加载IEEE33节点算例省去手工录入线路阻抗和节点负荷的繁琐工作。如果你习惯用MatlabMatpower里的case33bw也是现成的。两种工具本质一样都是先把网络数据读进来再运行牛顿拉夫逊潮流求解。import pandapower as pp import pandapower.networks as nw # 加载IEEE33节点标准算例 net nw.case33bw() # 检查网络基本参数 print(节点数:, len(net.bus)) print(支路数:, len(net.line)) print(变压器数:, len(net.trafo)) # 运行交流潮流 pp.runpp(net, algorithmnr) # 查看节点电压结果 print(net.res_bus.loc[:, [vm_pu, va_degree]])跑通之后你应能看到各节点电压在0.9到1.0pu之间末端节点电压明显偏低这符合33节点系统的经典特性。如果末端电压低于0.9多半是基准参数不对或者联络开关状态异常先别急着算碳势回头检查网络数据。3.2 从潮流结果提取支路有功pandapower潮流跑完支路有功功率存在net.res_line里字段是p_from_mw方向是从from_bus流向to_bus。这个方向很关键碳势方程里判断谁来流全靠它。一条支路的有功功率可能是负值表示实际潮流从to_bus流向了from_bus和建模时定义的方向相反。import numpy as np line_res net.res_line[p_from_mw].values from_bus net.line[from_bus].values to_bus net.line[to_bus].values n_bus len(net.bus) # 初始化节点有功注入累积量 inflow np.zeros(n_bus) # 遍历支路按实际潮流方向累计注入 for f, t, p in zip(from_bus, to_bus, line_res): if p 0: # 实际从f流向t inflow[t] p else: # 实际从t流向f inflow[f] - p print(节点注入有功分布:, inflow)这里我踩过不小的坑。第一次算的时候没有考虑支路功率方向想当然把所有支路都当成从序号小的节点流向序号大的节点结果算出来的碳势完全乱掉高碳区域莫名其妙出现在末端。后来才意识到IEEE33系统的线路虽然有编号方向但潮流实际流向只由电压相角和线路阻抗决定。所以计算碳势前务必先做一次方向判断。3.3 碳势方程组求解实现有了节点注入有功下一步就是把第2部分的方程组变成代码。我习惯把过程中间结果全部打印出来方便和手工推导值对照。def compute_nodal_carbon(n_bus, from_bus, to_bus, line_res, gen_power, gen_carbon): # 统计各节点实际来流功率 inflow np.zeros(n_bus) for f, t, p in zip(from_bus, to_bus, line_res): if p 0: inflow[t] p else: inflow[f] - p # 构建线性方程组 A E b A np.eye(n_bus) b np.zeros(n_bus) for i in range(n_bus): total_in inflow[i] gen_power[i] if total_in 1e-9: # 孤立节点或纯无源节点跳过 continue # 统计所有流入节点i的支路 for f, t, p in zip(from_bus, to_bus, line_res): if p 0 and t i: # 功率从f流向i A[i, f] - p / total_in elif p 0 and f i: # 功率从t流向i A[i, t] - (-p) / total_in # 发电碳势项 b[i] gen_power[i] * gen_carbon[i] / total_in E np.linalg.solve(A, b) return E在这个函数里gen_power是每个节点的发电有功功率数组gen_carbon是每个节点发电机对应的碳排放因子数组。对于IEEE33节点系统我把外电网等效成节点0上的一个发电机功率设为平衡节点注入功率碳因子设为0.581tCO2/MWh其他节点如果有光伏就在对应节点设置gen_power为该光伏出力gen_carbon设为0。求解完成后E数组就是33个节点的碳势值单位是tCO2/MWh乘1000就是gCO2/kWh。我通常会加一句校验把所有节点的负荷乘碳势再求和应该等于所有发电碳流量之和忽略网损时近似相等。如果不满足这个守恒关系说明某个支路方向判断错了或者某个节点注入功率漏统计了。3.4 拓扑可视化实现碳势算出来是一堆数字不画图根本看不出空间规律。可视化我分两层做第一层是拓扑图把节点连接关系画出来节点颜色映射碳势大小第二层是柱状图直接比较33个节点的高低。先用networkx构建拓扑import networkx as nx import matplotlib.pyplot as plt G nx.Graph() for i in range(n_bus): G.add_node(i) for f, t in zip(from_bus, to_bus): G.add_edge(f, t) # 使用相对布局便于看清辐射状结构 pos nx.spring_layout(G, seed42, k0.8) nodes nx.draw_networkx_nodes( G, pos, node_colorE, node_size400, cmapRdYlGn_r ) nx.draw_networkx_edges(G, pos, alpha0.3, width1.2) nx.draw_networkx_labels(G, pos, font_size8) plt.colorbar(nodes, label节点碳势 (tCO2/MWh)) plt.title(IEEE33节点碳势分布图) plt.axis(off) plt.tight_layout() plt.savefig(ieee33_carbon_map.png, dpi300) plt.show()这个图上能很直观看到碳势从馈线出口往末端的变化。如果所有支路都不存在分布式电源外电网来的碳势在全网基本保持一致这是正常现象因为馈线上没有额外低碳源注入。一旦在某个末端节点接入光伏光伏下游节点的碳势就会明显下降而光伏上游节点碳势保持不变这个分界效果在图上非常清晰。柱状图我用Matplotlib的bar函数按碳势从高到低排序方便快速找出最高和最低节点order np.argsort(E)[::-1] labels [f节点{i} for i in order] values E[order] fig, ax plt.subplots(figsize(12, 5)) ax.bar(labels, values * 1000, color#2E86AB) ax.set_ylabel(节点碳势 (gCO2/kWh)) ax.set_xticks(range(n_bus)) ax.set_xticklabels(labels, rotation90) plt.title(IEEE33节点碳势排序) plt.tight_layout() plt.savefig(ieee33_carbon_bar.png, dpi300) plt.show()如果想把结果做成可交互的网页推荐用pyecharts把节点坐标和碳势值转成graph或map类型的图表鼠标悬停就能看到节点编号、碳势数值和对应负荷。交互版本对汇报和演示的加分很大但分析阶段我还是先用静态图因为信息密度高、出图快。3.5 分布式电源接入对比实验只算一次碳势项目价值有限我建议至少做一个对比实验在重载末端节点18接入一台光伏或风电机组观察全网碳势变化。做法很简单在pandapower网络里添加一个gen到目标节点给定出力重新跑潮流再用同一个碳势函数计算。pp.create_gen(net, bus18, p_mw0.3, namePV_18) pp.runpp(net, algorithmnr)这时再看节点碳势节点18下游或者附近节点的碳势会明显下降离光伏越近的负荷节点受益越大。这就是碳排放流理论的实际意义它能把新能源降低碳排放这个宏观结论量化到具体节点、具体用户。你也可以反过来做在节点18接入一个燃气轮机碳势给定0.2tCO2/MWh看它和外部电网0.581的排放因子混合后本地下游节点的碳势变化。这种对比实验很好写进论文或项目汇报。4. 常见问题与排查技巧实录4.1 潮流不收敛怎么办IEEE33节点系统本身很成熟如果pandapower或Matpower标准数据直接跑通常一次就能收敛。但如果手动录入线路参数最常见的问题是把电抗和电阻的单位搞错或者把联络开关全部闭合导致弱环网出现潮流反转。遇到不收敛先把所有联络开关打开恢复成标准辐射状网络再逐步闭合测试。如果还是收敛不了把算法换成runpp的iwamoto_nr或者在潮流前调整电压初始值。碳势计算对潮流收敛质量的要求比较高支路功率如果有较大误差碳势结果可能出现负值或者超过1的异常值。我做验证的时候会把潮流结果的节点功率不平衡量打出来确保误差在1e-6量级再继续。4.2 支路方向判断和碳势方向不一致这是碳势计算独有的坑。支路建模方向是从from_bus到to_bus但实际有功功率可能是反向流动的。在环形网络中同一条支路在不同运行方式下方向还可能反转。如果你的碳势代码没有按实际方向动态判断结果就会错。我的排查技巧是把每条支路的方向和功率打印出来画在有向图上肉眼扫一遍是否合理。比如IEEE33系统根节点0到下游的功率都应该是正值一旦出现负值先看是不是联络开关改变了潮流路径。另外还要注意碳势方程里的来流支路是针对碳势流向说的。如果一个节点既是某条支路的受端又是另一条支路的送端那它只统计流入自身的支路不统计流出支路。用矩阵法写代码时这个逻辑非常容易混建议在函数里单独加一个inflow统计再基于它构建方程组。4.3 所有节点碳势都为0出现这个问题的概率相当高原因基本只有一个没有给外部电网设定碳排放因子。IEEE33节点系统默认只有一个平衡节点作为电源其他都是负荷和线路。如果平衡节点上的发电机碳势没有赋值程序默认值就是0整个系统的碳势自然全是0。解决方式就是把松弛节点看成一台火电或区域电网等效电源赋上区域平均碳因子。我的习惯是在数据准备阶段单独建一个carbon_map字典显式记录每个电源节点的碳势比散落在代码里更不容易漏。4.4 碳势出现负值或数值跳变如果碳势出现负值先查发电碳势是否设置了负值或者支路功率方向判断里有负负得正的逻辑错误。另一种情况是高碳和低碳电源通过环网相连两种碳流混合过程如果功率很小数值计算上可能出现条件数很大的矩阵。这时候别急着调公式先用numpy的cond函数检查方程组系数矩阵的条件数。如果条件数超过1e12优先考虑是不是网络中存在非常小的闭环潮流把对应小功率支路做阈值过滤比如低于1e-6MW的支路功率直接视为0往往能解决问题。4.5 可视化图太乱看不出规律33节点全画在一个力引导布局里节点一多标签重叠是必然的。我的经验是第一节点编号不用全显示只显示碳势排名前五和后五的节点第二节点颜色用连续色条配合图例就能表达全部信息第三在图上叠加一个文字标注记录运行工况比如联络开关闭合或节点18接入300kW光伏避免图做好了却忘了是什么工况下跑的。如果你需要发布到网页pyecharts的graph类型可以支持缩放和拖拽比静态图友好太多。5. 项目落地的一些体会整套流程跑下来最让我觉得有用的不是碳势公式本身而是把潮流计算、网络建模、数据可视化三个模块串成了一个可以反复使用的框架。现在再拿到IEEE123节点甚至更复杂的馈线模型我只需要换掉pandapower的网络加载函数碳势计算和后续可视化代码几乎不需要改动。这种核心算法与网络数据解耦的写法才是这个项目最有迁移价值的部分。有个细节想单独提一下。最开始我把碳势计算函数写成了递推版只适用辐射状网络联络开关一闭合结果就乱了。改成线性方程组求解之后不管是辐射网还是环网代码结构都一样鲁棒性提升很多。所以如果你要做这个方向我建议一开始就按方程组来写别贪图递推实现简单而给自己埋坑。计算速度上33节点系统用numpy线性求解器基本上是毫秒级完全不需要担心性能问题。后续如果你想在这个基础上做更深的内容可以对标节点碳势时序分析这个方向。给每个节点挂24小时或全年8760小时的负荷曲线光伏风电也按时序出力逐时段计算碳势最后得到一张节点碳势随时间变化的热力图。哪个节点在哪个时段碳势最高一目了然这就能给需求侧响应、储能充放电策略和低碳调度提供很实际的参考。我在做时序扩展时最大的体会是不要试图把碳势计算嵌套在潮流迭代里而是先批量跑完潮流把支路功率统一存成数组再逐时段碳势计算。这样逻辑清晰后期加数据也方便。最后给新手一个建议第一次跑通后先别急着接分布式电源和环网工况就把标准的单电源辐射状IEEE33节点跑一遍打印出每个节点的碳势手工拿公式验证其中两三个节点的计算过程。这一步校验过关了再把代码往复杂场景推否则错误藏在混合工况里排查会非常痛苦。这套流程做到位节点碳势计算和可视化这项技能基本就能直接沉淀到你的工具箱里了。