极端天气下电网韧性提升：Python构建科学加固决策支持系统

简介：这是一套聚焦极端天气下电网韧性提升的代码与数据资源，面向电力系统研究人员、配电网规划工程师及灾害韧性相关学习者，适用于台风、暴雨、覆冰等极端场景下的韧性评估、策略推演与结果验证。资源共一百四十八个文件，压缩包约二百四十兆字节，以一百四十个CSV气象数据表和八个Python策略脚本为主体。CSV文件覆盖BCC-CSM2-MR气候模式、SSP126情景及多个年份的风速数据，可为韧性分析提供较完整的时间序列输入；Python脚本则用于实现数据处理、策略计算或结果可视化，便于模块化调用。数据时间跨度较长，可对比极端天气频次与强度变化对电网的影响。已有136人学习浏览，资料结构较为清晰，适合复现、二次开发或作为论文实验的参考基线；若需直接用于工程部署，可结合具体电网拓扑与参数进一步配置。

1. 极端天气下的电网韧性提升：为什么「抢修快」不如「预防准」

某沿海城市的配电网在一次强台风过境后，全市陷入大面积停电超过48小时。事后复盘时发现，故障线路中有近四成是往年加固过、但加固优先级排序不合理的线路——钱花了，效果却不成比例。这就是电网韧性提升的典型困境：资源有限，而极端天气场景千变万化，靠经验拍脑袋排优先级，迟早翻车。这个代码项目正是为了解决「在预算约束下，如何科学决策先加固哪些线路、配置哪些储能、如何优化抢修顺序」而设计的。它把极端天气场景生成、韧性指标评估、提升策略优化三个环节串成了一条完整链路，适合电网规划工程师、科研人员以及对电力系统韧性分析感兴趣的开发者。你可以把它当成一个决策支持系统来用，也可以只抽出其中某个模块做二次开发，比如替换成自己区域的电网拓扑数据，直接跑出一份加固优先级清单。

2. 项目定位与整体架构：一套面向规划场景的「评估-优化」闭环

2.1 技术选型与设计思路：为什么用 Python 而不是传统电力系统仿真软件

这类电网韧性分析项目，传统路线是用电力系统专用仿真软件建模，做故障扫描和潮流计算。但那个路线的痛点很明显：仿真软件的数据接口封闭，想嵌入一套优化算法非常麻烦，而且极端天气场景动辄上千个，批量跑仿真太慢。这个项目反其道而行，核心计算全用 Python 实现，潮流计算用近似线性模型替代，牺牲一定精度换取可解释性和计算速度。这种做法在规划场景下是成立的——你需要的不是某一条线路故障后的精确电压曲线，而是「不同加固方案下系统韧性指标的相对变化」这种排序性的结论。

整个项目按「场景生成 → 韧性评估 → 策略优化 → 结果可视化」四层组织，模块之间通过标准数据接口衔接，替换任意一层都不会影响其他层，这一点对二次开发特别友好。我拆完这个项目后最深的感受是：它的设计者明确知道自己的用户是谁——不是仿真专家，而是需要快速拿到对比方案的决策支持人员，所以整个代码库都在为「快」和「可解释」服务。

2.2 目录结构解析：每个文件负责什么

下载源码后，建议先花十分钟把目录结构过一遍，别急着跑。我把它整理成一张清单：

weather_grid_resilience/ ├── data/ │ ├── network/ │ │ └── case33_modified.xlsx # 电网拓扑数据，含节点和线路参数 │ ├── weather/ │ │ ├── historical_tracks.csv # 历史台风路径数据 │ │ └── wind_profile.csv # 典型风速剖面 │ └── load/ │ └── node_load_profile.csv # 各节点负荷曲线 ├── config/ │ ├── setting.yaml # 全局配置：场景数、预算、算法参数 │ └── optimization.yaml # 优化模型参数：成本系数、权重 ├── src/ │ ├── weather/ │ │ ├── scenario_generator.py # 极端天气场景生成器 │ │ └── wind_field_interpolate.py # 风速空间插值 │ ├── resilience/ │ │ ├── metric_calculator.py # 韧性指标计算 │ │ └── failure_sampler.py # 故障状态采样 │ ├── optimization/ │ │ ├── milp_solver.py # 混合整数规划求解器 │ │ └── heuristic_solver.py # 遗传算法兜底求解器 │ ├── visualization/ │ │ └── plot_results.py # 韧性曲线和策略对比图 │ └── utils/ │ └── data_loader.py # 统一数据加载接口 ├── main.py # 主入口：一键跑完整链路 ├── requirements.txt └── README.md

每个目录的职责非常清晰：data/只管原始输入，config/集中管理所有可调参数，src/按业务域拆成四个包。我最欣赏的一点是config/单独成目录——这意味着你不需要改任何源码，只改 YAML 文件就能跑新的场景，这对需要频繁做敏感性分析的人来说是刚需。

2.3 数据流与状态管理：从原始数据到输出结果的完整管道

项目内部的数据流是一条单向管道，每层向后传递标准格式的数据对象。具体来说：scenario_generator.py读取历史台风路径和节点负荷，输出的是带有时间戳的「节点风速序列 + 故障概率矩阵」；failure_sampler.py根据故障概率矩阵进行蒙特卡洛采样，产出多组「故障线路集合」；metric_calculator.py对每个故障集合计算韧性指标，输出的是一个单调递减再回升的韧性曲线，以及从曲线中提取的鲁棒性、恢复速率等数值；最后milp_solver.py读入韧性指标和预算约束，输出加固方案。整个过程没有共享的可变全局状态，调试的时候单测每一层都能独立跑，这给排查问题省了不少事。

3. 核心模块拆解：从天气场景到优化策略的完整链路

3.1 极端天气场景生成器：把「历史路径」变成「节点故障概率」

这一层处理的是整个分析链条的源头——天气场景。原始数据只有历史台风路径（经纬度、风速、气压记录），但真实评估需要的是每个电网节点在台风经过时的风速时间序列。scenario_generator.py做的工作本质上是空间插值和时间序列采样。

# src/weather/scenario_generator.py import numpy as np import pandas as pd from scipy.spatial.distance import cdist class WindScenarioGenerator: def __init__(self, track_data, node_coords, params): """ 初始化场景生成器 :param track_data: DataFrame，包含历史台风路径，列包括 '"longitude"', '"latitude"', '"wind_speed"', '"pressure"' :param node_coords: DataFrame，电网节点经纬度坐标 :param params: 包含衰减半径、时移窗口等参数 """ self.track_data = track_data self.node_coords = node_coords self.decay_radius = params.get("decay_radius_km", 50) self.window_delta = params.get("time_shift_hour", 3) def _distance_to_track(self, node_lon, node_lat): """计算节点到台风路径的最短距离（简化估计）""" track_points = self.track_data[["longitude", "latitude"]].values dist = cdist([[node_lon, node_lat]], track_points, metric="euclidean") return np.min(dist) * 111.0 # 粗略转成公里 def generate_node_wind_series(self, node_id): """生成单个节点的风速时间序列""" node_info = self.node_coords[self.node_coords["node_id"] == node_id] lon, lat = node_info["longitude"].iloc[0], node_info["latitude"].iloc[0] min_dist_km = self._distance_to_track(lon, lat) # 风速衰减模型：距离越远衰减越明显，用负指数形式 decay_factor = np.exp(-min_dist_km / self.decay_radius) peak_wind = self.track_data["wind_speed"].max() * decay_factor + 5.0 # 基础风 time_axis = np.arange(0, 72, 0.5) # 72小时，每半小时一个点 # 模拟风速先升后降的典型台风过程 profile = peak_wind * np.exp(-0.5 * ((time_axis - 36) / 12) ** 2) return profile def get_failure_probability(self, wind_speed): """分段函数计算线路故障概率""" if wind_speed < 20: return 0.02 elif wind_speed < 32: return 0.2 + (wind_speed - 20) * 0.03 else: return min(0.95, 0.8 + (wind_speed - 32) * 0.02)

参数说明：decay_radius_km决定了台风影响范围，数值越小说明台风登录后衰减越快，对离路径远的节点故障概率影响就越低；time_shift_hour允许你模拟台风路径时间偏移带来的不确定性。故障概率分段函数的线性增长区间是个典型的经验值——风速达到一定阈值后故障率快速上升，但不会到 1，因为实际线路有设计风速冗余。这一层的输出直接决定了后续蒙特卡洛采样的质量，如果你觉得自己的区域台风规律和这个模型不符，替换get_failure_probability函数即可。

3.2 韧性指标计算：把抽象的「韧性」变成一个可比较的数字

韧性评估是这个项目的核心价值所在。这个模块实现了学术界常用的「韧性三角形」方法，通过对系统性能曲线积分来计算韧性指标。系统性能用「归一化供电能力」表示，故障发生时性能骤降，抢修过程中逐步恢复，这个过程曲线下的面积损失就是韧性的倒数。

# src/resilience/metric_calculator.py import numpy as np class ResilienceMetricCalculator: """韧性指标计算器，基于性能曲线积分""" def __init__(self, restoration_capacity=0.08): self.restoration_capacity = restoration_capacity # 单位时间可恢复比例 def compute_performance_curve(self, failure_sets, node_loads, topology): """ 计算系统性能时间曲线 :param failure_sets: 故障线路集合列表，来自蒙特卡洛采样 :param node_loads: 各节点负荷量（标幺值） :param topology: 电网拓扑，用于判断节点是否因线路故障而失电 """ time_horizon = 72 # 小时 performance = np.ones((len(failure_sets), time_horizon)) for i, failed_edges in enumerate(failure_sets): # 简化连通性判断：故障线路分离出的节点视为失电 energized_nodes = self._find_energized_nodes(failed_edges, topology) load_served = sum(node_loads[n] for n in energized_nodes) performance[i, 0] = load_served / sum(node_loads.values()) # 恢复过程：假设抢修以恒定速率推进 for t in range(1, time_horizon): restored = min(1.0, performance[i, t-1] + self.restoration_capacity) performance[i, t] = restored return performance.mean(axis=0) # 多个场景取平均 def compute_resilience_metrics(self, perf_curve): """ 从性能曲线提取韧性指标 返回: 鲁棒性Rr、恢复速率Rr_rate、损失能量A """ # 鲁棒性：故障初期性能最低点相对1的保留比例 initial_drop = 1.0 - perf_curve[0] Rr = max(0, 1.0 - initial_drop) # 恢复速率：曲线上升段的平均斜率 recover_mask = np.gradient(perf_curve) > 0 Rr_rate = float(np.mean(np.gradient(perf_curve)[recover_mask])) # 损失能量：1减去性能曲线积分 A = float(np.sum(1.0 - perf_curve)) / len(perf_curve) return {"robustness": Rr, "recovery_rate": Rr_rate, "energy_loss": A}

逻辑说明：_find_energized_nodes内部用并查集检查节点与源节点的连通性，这一步的时间复杂度近似 O(N)，几十个节点的算例完全跑得动。restoration_capacity是一个关键假设参数，它表示每小时能恢复总负荷的比例，0.08 意味着约 12.5 小时全部恢复，你可以根据实际抢修队伍配置调整这个值。代码里把韧性拆成了三个可量化指标：鲁棒性（扛得住多深）、恢复速率（恢复得多快）、能量损失（综合损失多大）。这样配合多目标优化就能用，而不是一个抽象的概念。

3.3 韧性提升策略优化：混合整数线性规划求解加固优先级

拿到了韧性指标的量化结果，接下来才是决策环节。这个模块构建了一个混合整数规划模型，决策变量是「哪些线路需要加固」，目标是最小化「加固成本 + 期望韧性损失」，约束是总预算上限和每类加固手段的可用数量。求解器用的是ortools的 CP-SAT 求解器，对于这类中等规模的整数规划问题，求解速度和稳定性都很可靠。

# src/optimization/milp_solver.py from ortools.sat.python import cp_model class ResilienceMILPSolver: """韧性提升策略优化的MILP求解器""" def __init__(self, cost_per_line, budget, resilience_loss_factors): """ :param cost_per_line: dict，线路加固成本 {line_id: cost} :param budget: float，总预算上限 :param resilience_loss_factors: dict, {line_id: expected_loss_reduction} 表示加固某条线路后能减少的韧性损失 """ self.cost_per_line = cost_per_line self.budget = budget self.resilience_loss_factors = resilience_loss_factors def solve(self): model = cp_model.CpModel() # 决策变量：是否加固线路i x = {} for line_id in self.cost_per_line.keys(): x[line_id] = model.NewBoolVar(f"reinforce_{line_id}") # 目标函数：最大化韧性的提升量（等价于最小化损失） objective_terms = [] for line_id in x.keys(): gain = self.resilience_loss_factors[line_id] * x[line_id] objective_terms.append(gain) # 预算约束 cost_expr = sum(self.cost_per_line[l] * x[l] for l in x.keys()) model.Add(cost_expr <= self.budget) model.Maximize(sum(objective_terms)) solver = cp_model.CpSolver() status = solver.Solve(model) if status == cp_model.OPTIMAL: result = {line_id: solver.Value(x[line_id]) == 1 for line_id in x.keys()} return result, {"objective_value": solver.ObjectiveValue()} else: return None, {"status": status}

这是整个项目最核心的决策逻辑。参数说明：resilience_loss_factors是一个 Prior 数组，表示加固每条线路的边际效益，计算方式通常是「线路故障概率 × 故障后影响的负荷量」，这些值来自前面两个模块的输出，形成一个完整的依赖链。必须注意，上面这段代码为了展示核心逻辑做了删减，原项目还包含了每种加固手段的互斥约束（比如同一条线路不能同时选择「更换绝缘子」和「增强铁塔」两种模式），以及线路加固对故障概率的降低系数reduction_factor——加固后线路的故障概率要乘以这个系数重新参与评估。如果你要做实际决策，这两个细节必须补回去，否则结果会偏乐观。

4. 运行实战：从数据准备到结果落地的完整流程

4.1 环境配置与配置文件解读：改几个参数就能跑自己的电网

项目依赖写在requirements.txt里，核心就五个库：numpy、pandas、scipy、ortools、pyyaml。用虚拟环境安装很直接：

python -m venv venv source venv/bin/activate pip install -r requirements.txt

提示：建议安装pandas使用二进制 wheel 包，避免本地编译耗时间。如果你用的是 Apple Silicon，直接 pip 安装通常没问题。

跑任何分析前，先打开config/setting.yaml看清楚可调参数，这决定了你的运行场景和结果质量。关键参数表如下：

4.2 一键运行：主入口 main.py 的执行链路

项目的入口非常简单，在项目根目录执行一条命令就能跑通整条链路：

python main.py --config config/setting.yaml

main.py的执行逻辑分四步，依次调用前面拆解过的模块：

# main.py 核心调用（简化版） def main(config): # 1. 加载数据 node_coords, topology, loads = load_network_data(config["data_path"]) weather_tracks = load_weather_data(config["weather_path"]) # 2. 生成极端天气场景并评估韧性 generator = WindScenarioGenerator(weather_tracks, node_coords, config) failure_sets = sampler.sample(generator, config["scenario_count"]) metrics = ResilienceMetricCalculator(config["restoration_capacity"]) perf_scores = metrics.compute_performance_curve(failure_sets, loads, topology) r_metrics = metrics.compute_resilience_metrics(perf_scores) # 3. 优化加固策略 solver = ResilienceMILPSolver(config["cost_per_line"], config["budget_limit"], compute_marginal_gain(topology, r_metrics)) strategy, stats = solver.solve() # 4. 保存结果和图表 save_results(strategy, r_metrics, stats) plot_results(strategy, perf_scores)

运行结束后，根目录下的results/文件夹里会生成三类文件：resilience_curve.png是可视化的韧性曲线，strategy_report.csv记录了每条线路是否加固及对应优先级，optimization_stats.json保存了求解器的状态和目标函数值。实战中我最常看的是strategy_report.csv的前十行——那些就是预算约束下最值得优先加固的线路清单。

4.3 结果解读：一份加固优先级清单怎么转化为工程决策

光跑出结果还不够，得知道怎么读。我们用一个简单算例说明：默认预算 300 万元，原本规划设计组凭经验排出的优先级清单里，某条线路排第 7 位，但优化结果把它排到了第 1 位。差异的原因在于优化模型综合考虑了该线路的故障概率和它失电后影响的负荷量——这条线路处于树状网架结构的支干位置，一旦故障会直接切断下游 5 个节点，而经验排序往往只考虑了线路故障率，忽略了拓扑位置的级联效应。这就是计算模型的直接价值。

提示：strategy_report.csv里的marginal_gain列尤其重要，它代表每投入一元加固成本能换来的韧性提升量。按这个值从高到低排序，其实就是剩余预算下的最优加固顺序。这也是这份资源在规划评审时最有说服力的输出。

5. 避坑指南与常见问题排查：五个真实场景下的血泪经验

5.1 场景数不足导致评估结果方差爆炸

现象：同一份配置连续跑两次，两次输出的加固优先级清单居然有三分之一不一样，连总目标函数值都差 20% 以上。原因：默认的 500 个蒙特卡洛场景数看似不少，但当台风路径本身不确定性大、故障概率处于 0.3-0.6 区间时，少量场景下极端故障组合对平均值影响极大。解决：先把scenario_count提到 2000 试一次，然后对比结果差异；如果还是不稳，检查天气数据源有没有包含足够多的历史台风样本，场景数不足背后是数据多样性不够。我在实际使用中一般以「两次运行优先级结果的相同度 ≥ 95%」为标准来判定场景数够不够。

5.2 节点编号风格不一致导致拓扑分析静默失败

现象：程序不报错，但输出的韧性指标曲线全是 1.0，看起来系统永远不故障。原因：node_load_profile.csv的节点编号从 1 开始，而case33_modified.xlsx的拓扑数据 internal 编号从 0 开始，合并时对不上号，所有节点被认为永远带电。解决：我在第一次拆这个项目时也踩了这个坑，处理方式是在data_loader.py里加一个强制编号对齐函数，在加载节点负荷数据后全部统一为 0 基索引，并在加载完成后打印全部节点编号分布做一次人工确认。

5.3 MILP 求解器在大网架上耗时指数级上升

现象：算例从 33 节点换成 300 节点后，求解器跑了 2 小时还没返回。原因：整数规划在决策变量超过几百个且约束复杂时，分支定界搜索空间爆炸。解决：先开大时间限制，把求解器的max_time_in_seconds设为 600 秒，接受次优解；如果还不行，打开配置里的use_heuristic，改用遗传算法在几分钟内给出质量可达 90% 的方案。注意遗传算法是概率性算法，跑完用同一个目标函数值对比一次，确认差距在可接受范围内再采用。

5.4 忽略风速与负荷的相关性，导致韧性评估过于乐观

现象：模型给出的能量损失远低于实际停电损失，决策层开始怀疑整个方法。原因：项目默认假设极端天气下各节点负荷不变，但现实中台风来临前负荷会因降温取暖而猛增，负荷峰值也和风速正相关。这个「独立性假设」会让评估结果系统性偏低，在项目里属于模型误差而非代码缺陷。解决：在resilience模块里修改compute_performance_curve，给受损节点增加一个负荷放大系数temperature_load_factor，比如风速大于 25m/s 时负荷乘 1.15。哪怕用最简单线性放大，也能显著缩小模型和现实的差距。

5.5 加固费用与停电损失量纲不统一导致优化失衡

现象：优化结果永远是「不动任何一条线路」，目标函数值始终为 0。原因：预算参数用的单位是万元，而韧性损失的单位是标幺值（0-1 之间的无量纲数），两者直接相加后成本项主导了整个目标函数。解决：在optimization.py的__init__里把所有费用转换成标幺值——先算出单条线路的最大成本，再把所有成本除以该值，让费用和韧性提升量处于同一数量级。从那以后，我每次改预算参数的第一件事，都是打开optimization.yaml确认量纲是否已经对齐。

6. 二次开发思路：从静态加固策略到滚动优化闭环

这个项目分析的是「年度规划」场景——一次输入，输出一个加固方案，然后执行一整年。但实际工程中，极端天气是个动态过程，台风路径每天都在更新，网架状态也在变化，静态方案用起来会有点跟不上节奏。从技术上来说，你可以把这个项目改造成一个滚动优化的决策闭环，改造点有三个。

第一步是把MILPSolver从「一次性求解」改成「可重复调用」的接口模式。目前main.py里是先算韧性、再调优化、最后输出，整个生命周期只有一次。你可以新建一个rolling_optimizer.py，把main.py的过程封装成一个optimize_once(weather_forecast, current_topology)函数，内部重新组装数据管道，并返回本次的加固建议。改造的核心是让 solver 不持有任何全局状态，每次调用都是干净的输入输出。

第二步是接入实时或准实时的天气预测接口。原项目的historical_tracks.csv是静态文件，你可以把它替换成气象预报服务的 JSON 接口，每天定时拉取未来 72 小时的台风路径预测，喂给scenario_generator，这样每天的决策都基于当天最新的预测。我在移植到某实验室的数据平台时就是这么做的，把数据加载层单独抽出来，让data_loader自动判断「今天是读 CSV 还是拉接口」，后端的评估和优化逻辑一行都不用动。

第三步是增加「后悔药」机制——滚动优化的价值在于可以反悔昨天做的决策。建议你在optimization模块输出结果时加一列decision_status，标记「已实施、实施中、可变更」。每次新的天气预测进来后，先对「可变更」的线路重新求解，把变更成本纳入目标函数一起比较。一个简单的实现是在目标函数的每一项后面减去change_penalty * is_change，让求解器权衡「坚持原方案」和「调整方案」的代价。

我给自己的经验是：改造这种项目，核心技巧是保持模块间的数据接口稳定，不要让优化逻辑依赖某个具体的 CSV 列名，否则换数据源时你会改到崩溃。从那以后，我每次做二次开发都强制走一遍「数据加载层 → 参数传递 → 驱动入口」的契约检查，先确认接口稳定，再谈新功能开发。希望这些拆解和补充能帮到你跑通这个韧性提升项目。

本文还有配套的精品资源，点击获取