做了这么多年配电系统规划，我越来越确信一件事：经济性和可靠性这对矛盾，躲是躲不掉的。这篇文章要聊的，就是一套基于经济与可靠性双目标的混合配电系统规划及可靠性评估方法，以及配套的Python代码实现思路。我会把建模逻辑、算法细节、代码骨架、算例结果和踩坑经历都摊开讲，适合正在做配电网规划、微电网规划，或者急着准备毕业设计、课题结题需要双目标优化加可靠性评估代码的同行参考。
先说清楚一个前提：很多刚接触这个方向的读者，拿到题目第一时间就去搜“NSGA-II代码”“蒙特卡洛Python库”，结果跑出来的结果自己都不敢信。问题往往不出在算法上，而是出在“混合配电系统”这个词的理解上。系统模型没搭对，再高级的算法都是在给错误模型做精装修。所以这篇文章先用一整个章节把规划对象拆干净，再谈目标函数、评估引擎、求解器、代码骨架和实战坑点。
1. 混合配电系统的规划对象拆解：电源、网架与负荷侧都在混合什么
1.1 三种“混合”的工程背景
“混合配电系统”在学术文献里有两种主流用法。一种指交直流混合的配电网络，交流馈线通过电压源换流器（VSC）或者智能软开关（SOP）与直流母线连接；另一种指多类型分布式电源、储能和传统配电网混合运行。现实中做规划项目时，这两件事经常同时出现——馈线上既有PV又有风机，还接了电池储能，局部直流区域通过换流器接入交流母线。
所以规划对象的“混合”至少要拆成三个维度，如下表所示。
| 混合维度 | 典型构成 | 规划阶段新增决策 | 对可靠性的影响方式 |
|---|---|---|---|
| 电源侧 | 分布式光伏、风电、微型燃气轮机、储能、EV充放电接口 | DG选址定容、储能位置与容量 | 故障后能否孤岛运行、出力波动带来的供电不确定性 |
| 网架侧 | 交流馈线、直流线路、VSC/SOP、联络开关 | 线路升级、联络位置、SOP容量 | N-1转供能力、故障隔离范围 |
| 负荷侧 | 柔性负荷、可削减负荷、EV充电桩 | 需求响应策略参数 | 紧急情况下可切负荷的灵活性 |
规划时不可能把所有变量都放进优化模型里。解空间爆炸不是一句玩笑，决策变量每多10个，进化算法的收敛难度上一个量级。我通常的做法是：把DG选址定容、储能容量、联络线改造作为主决策变量，SOP容量和需求响应参数先固定成场景输入，等选出主方案后再做校核调整。这样模型的物理意义清晰，计算量也可控。
1.2 为什么单目标规划做到最后总差口气
单目标优化不是不能做，问题是“只算经济账”和“只算可靠性账”都会得出工程上很难接受的方案。
只算经济账时，优化算法会倾向于少装DG、不建联络线、变压器和线路全部卡着容量上限走。结果是网损和购电成本确实低，但第一条馈线故障就可能让一大片负荷停电。反过来，只算可靠性账，模型会疯狂加联络线、冗余设备，把所有元件可靠度往高了抬，造价高到项目根本批不下来。
我在实际项目中见过一个很典型的例子：某个园区配网，经济最优方案和可靠性最优方案的年费用能差到1.8倍，但两者对应的EENS差了6倍以上。这种矛盾完全是结构性的，不是调调参数就能消掉的。所以“经济与可靠性双目标”不是论文里的形容词，而是工程决策的刚需。
1.3 本文规划问题的决策变量定义
为了让后面的模型和代码有明确抓手，我把规划问题定成如下形式：
- 候选DG节点集合内，为每个节点选择是否安装DG以及安装容量等级（离散档位）。
- 候选储能节点集合内，为每个节点选择储能容量（离散档位）。
- 待升级线路集合内，选择是否升级导线截面或改为电缆。
- 待选联络位置集合内，选择是否新增联络开关。
编码形式上，一个个体可以表示成整数数组，比如[DG节点1档位, DG节点2档位, ..., 储能节点1档位, ..., 线路1是否升级, 线路2是否升级, ...]。这种编码用NSGA-II处理起来非常顺，交叉变异都能对应到具体工程动作。
有一点要提前打招呼：DG容量档位不是越大越好。渗透率过高会带来反向潮流、电压抬升和继电保护配合问题。所以后面目标函数和约束里，我会专门加渗透率约束，避免算法把全系统都铺满光伏板。
2. 经济目标与可靠性目标如何在同一个优化框架里共存
2.1 先避开一个常见误区：不要急着把可靠性折算成钱
很多论文做“单目标”规划时，会把停电损失乘以单位缺电成本并入年费用，然后说自己是经济性最优。这种做法理论没问题，但工程上非常容易被质疑，因为单位缺电成本（VOLL）在不同用户类型、不同时段、不同地区差异极大，从几块钱到几十块钱每千瓦时都有。你取个中间值，决策结果往往就跟着变了。
双目标的好处在于：经济指标就是经济指标，可靠性指标就是可靠性指标，两者不强行换算。最后给决策者一条Pareto前沿，让领导在“多花钱买多少可靠性”这个层面上做选择，比在模型里偷偷塞一个电价参数要讲得清楚。
2.2 经济性目标：等年值费用模型
经济性目标F1采用等年值法，把寿命周期内的一次性投资折算到每年，再叠加运行维护、网损和购电成本。
F1 = CRF × C_inv + C_om + C_loss + C_buy
其中CRF是资金回收系数：
CRF = r × (1 + r)^n / ((1 + r)^n - 1)
r为贴现率，n为设备经济寿命。设备投资C_inv包含DG、储能、线路改造、联络开关等所有一次性费用；C_om按投资比例取一个年维护系数，光伏、储能、线路的维护比例各不相同；C_loss按典型日仿真得到的全年网损电量乘以购电价；C_buy是系统向上级电网购入的电量费用。
这个公式里最容易出错的是C_loss和C_buy。规划方案变了，DG装多了，网损下降了，但倒送功率和购电结构也变了。所以严格来说，F1的每一项都要在给定方案下跑一次潮流或时序仿真才能算准。后面代码章节会说明怎么用简化方法把这个过程控制在可接受的计算量内。
2.3 可靠性目标：用EENS作为核心指标
可靠性目标F2我推荐直接用年缺供电量EENS（兆瓦时/年），而不是只用SAIDI或SAIFI。SAIDI和SAIFI是用户视角的指标，需要用户数加权，数据不全时容易失真；EENS是能量视角的指标，跟负荷功率直接挂钩，规划优化时更方便计算。
负荷点i的年停电时间Ui由该点供电路径上所有串联元件的不可用率累加得到。元件不可用率的简化表达式为：
U ≈ λ × r / 8760
λ是故障率（次/年），r是平均修复时间（小时/次）。这是非序贯计算的根基。如果考虑转供和孤岛，Ui还要减去故障后通过联络线恢复供电的那部分时间，这部分必须靠网络分析算出来，没法用单纯串并联公式代替。
F2 = EENS = Σ (负荷点有功 × Ui)
派生指标SAIDI、SAIFI、CAIDI可以在算完Ui后很轻松地算出来，作为报告输出。
2.4 约束条件的工程清单
优化模型里要满足的约束，我一般按下面这个清单来列。
- 潮流约束：交流或线性化潮流方程约束，保证方案在稳态下可运行。
- 节点电压约束：各节点电压幅值在0.93~1.07 pu内，典型配电网工程要求。
- 支路容量约束：所有馈线、变压器、换流器不过载。
- DG渗透率约束：总DG出力占年峰值负荷的比例不超过设定值（比如40%）。
- 储能SOC约束：充放电功率和荷电状态在物理允许范围。
- N-1转供约束：任意一条支路或变压器退出后，失负荷量不超过允许值。
N-1转供约束和F2目标有耦合。如果约束卡得太严，比如要求所有负荷都能转供，那联络线投资会被推得非常高；如果放得太松，EENS指标又很难看。实际操作中我会把“允许最大失负荷比例”作为一个可调参数放进算例，而不是写死。这也是双目标规划比单目标更接近工程决策的地方。
3. 可靠性评估的蒙特卡洛引擎：元件停运模型与故障后果分析
3.1 元件可靠性参数怎么给
可靠性评估的精度上限，其实在参数输入那一刻就决定了。进化算法再聪明，也没法把垃圾参数变成可信结果。
配电网常用元件的可靠性参数大致在一个很稳定的范围，我列一组工程上比较常用的参考值：
| 元件 | 故障率λ | 平均修复时间r | 说明 |
|---|---|---|---|
| 架空线路（每公里） | 0.06~0.15次/年·km | 3~6小时 | 随地区雷暴、树障情况浮动 |
| 电缆线路（每公里） | 0.01~0.04次/年·km | 6~24小时 | 故障难定位，修复时间较长 |
| 配电变压器 | 0.005~0.02次/年 | 4~8小时 | 油变、干变差异较大 |
| 分布式光伏模块 | 0.05~0.15次/年 | 5~10小时 | 逆变器是主要薄弱点 |
| 储能系统（含PCS） | 0.2~0.5次/年 | 10~20小时 | 年检与维修周期较长 |
注意线路故障率一定要乘线路长度，单位是“次/年·km”。如果直接把0.1当作整条2公里线路的故障率，可靠性指标会被低估一半左右，这是新手最容易翻车的第一个点。
3.2 序贯蒙特卡洛与非序贯蒙特卡洛怎么选
蒙特卡洛模拟在可靠性评估里分两派。
- 非序贯蒙特卡洛（状态抽样）：对每个元件按可用率抽样一个0/1状态，得到一个系统状态快照，然后对这个快照做故障后果分析。优点是快，适合在进化算法内部做大规模候选方案的粗筛。
- 序贯蒙特卡洛（时序抽样）：对每个元件抽样“运行持续时间—修复持续时间”交替序列，合并成系统状态时间线，逐时间段做后果分析。优点是可以精确统计每次停电持续时间和缺供电量，适合最终方案的精细化验证。
工程上，我会在NSGA-II进化过程中用非序贯版本，只计算每个方案的期望缺供电量趋势；等Pareto前沿收敛之后，再对前沿上的十几个解跑一遍序贯蒙特卡洛，作为最终报告数据。这样做计算成本和结果精度都能兼顾。
非序贯状态抽样在Python里就是这么一句：
import numpy as np rng = np.random.default_rng(42) # availability = 可用率，例如0.99 state = rng.random(n_elements) > (1 - availability)序贯蒙特卡洛的核心逻辑是先给每个元件抽样第一次故障前的运行时间，然后模拟时间推进：
def sequential_mc(system, years=50, seed=42): rng = np.random.default_rng(seed) for br in system.branches: br.ttf = rng.exponential(1.0 / br.lam) # 运行时间 br.ttr = rng.weibull(1.2, system.rreat ) # 修复时间，weibull更贴近实际 current_time = 0.0 horizon = years * 8760.0 while current_time < horizon: # 找到下一个发生事件的元件 next_br = min(system.branches, key=lambda b: b.ttf) current_time += next_br.ttf # 执行故障后果分析，记录负荷点停电信息 analyze_outage(system, next_br, current_time) # 采样下一次故障/修复循环 next_br.ttf = rng.exponential(1.0 / next_br.lam) next_br.ttr = rng.weibull(1.2, system.rreat) current_time += next_br.ttr代码里用weibull替代指数分布做修复时间，是实测后我觉得更贴近真实的一个改动。后面的坑点章节我会单独展开讲。
3.3 故障后果分析：转供、孤岛与切负荷
蒙特卡洛抽样只是告诉你“这个元件坏了”，真正决定停电多严重的是故障后果分析。这部分要做三件事。
第一，把故障元件隔离，重新分析网络连通性。用深度优先搜索或广度优先搜索找出所有失电节点。第二，检查失电区域能不能通过联络开关合闸转供，如果联络线容量够，失电负荷可以恢复。这里要同时校验转供路径上的支路容量和电压约束。第三，如果不能完全转供，就按负荷优先级切负荷，直到满足约束为止。
如果失电区域内恰好有储能和分布式电源，可以进入孤岛运行模式。规划阶段做可靠性评估时，我一般只做有功平衡的简化校验——孤岛内DG和储能的总有功大于等于孤岛负荷，就认为可以维持供电；更精细的频率跌落动态校验放在后续设计阶段做，否则蒙特卡洛的每一轮都做电磁暂态仿真，计算量完全不可接受。
3.4 收敛判据：别等到天了才发现没收敛
蒙特卡洛模拟要知道什么时候停。工程上用方差系数β来衡量EENS估计的稳定性：
β = σ / (μ × √N)
σ是样本标准差，μ是均值，N是抽样次数。一般要求β小于2%~5%，指标越重要的模拟，阈值卡得越严。实际操作中每评估完一次，就把历年数据重新算一遍β，如果没达标就继续追加抽样。要注意不同方案之间β收敛速度差异很大，规划算例里不能固定写成“模拟1万次”，而要写成“模拟到β小于阈值为止”。
4. 双目标求解器的工程化细节：NSGA-II之外你必须处理的四件事
4.1 为什么是NSGA-II而不是其他算法
双目标优化算法可选的很多：MOEA/D、SPEA2、NSGA-III，甚至粒子群改一改也能跑。选择NSGA-II在这个场景里理由非常务实：两个目标用不到NSGA-III的高维参考点机制，MOEA/D的权重向量在离散整数编码下容易失效，SPEA2写起来比NSGA-II啰嗦。NSGA-II的快速非支配排序加拥挤距离，两个目标下表现稳定，工程案例和文献都能对上号，出问题也好查。
4.2 编码、修复与约束处理
前面说过个体用整数数组编码。这里有个关键处理：交叉变异后子代可能不满足工程约束，比如某个节点装了DG但又选了升级线路导致投资重复计算，或者总渗透率超限。两种处理思路，一种是罚函数，把约束违反对应一个很大的惩罚项塞进目标；另一种是修复，直接把非法基因改成合法值。
我更推荐修复为主、罚函数为辅。修复能保持种群的搜索效率，罚函数会给进化过程带入大量“名义上可行实际却偏离真实目标的个体”，极大的浪费评估次数。比如渗透率超限时，直接把某些DG档位降档，直到满足约束为止，这个操作本身就能保证后续评估都是在真实可行方案上进行的。
4.3 SBX交叉与多项式变异的参数选择
NSGA-II里最常用的实数算子是模拟二进制交叉（SBX）和多项式变异。虽然是实数算子，但在离散编码上仍可使用，只是交叉变异后要取整。
| 参数 | 常用范围 | 我的建议 |
|---|---|---|
| 交叉概率 Pc | 0.8~0.95 | 0.9 |
| 变异概率 Pm | 1/chromosome_length | 1/编码长度，别拍脑袋设0.1 |
| SBX分布指数 ηc | 10~20 | 15，越大子代越接近父代 |
| 多项式变异分布指数 ηm | 20~30 | 20，越小越容易产生远距扰动 |
还有一个容易被忽略的点：整数编码的变异不能简单地对实数加高斯噪声再取四舍五入。比如DG容量档位只有0/1/2/3四档，多项式变异产生的新值可能在-0.2附近，取整后变成0，导致变异失效。处理办法是变异后做一个映射，把数值夹到合法档位区间内再取整。
4.4 评估预算分配：粗评估和细评估怎么配合
NSGA-II一代种群200个个体，进化200代，每个个体做一次0.5秒的可靠性评估，总时间就是200×200×0.5，约5.6小时。这还只是理想情况。很多系统单次可靠性评估要跑几秒甚至几十秒，纯嵌套优化根本跑不完。
我的做法是把评估分两档跑。进化过程中用非序贯蒙特卡洛，抽样次数压到2000到5000次，方差系数可能到10%左右，但足以比较个体之间的优劣。等收敛后，取Pareto前沿上30到50个个体，换序贯蒙特卡洛加5万次抽样重新评估，得到最终报告数据。这个方法在多个算例里都被验证过：前沿个体排序和粗评估排序基本一致，但计算时间能省掉一个数量级。
5. Python实现骨架：从拓扑建模到Pareto迭代的代码结构
5.1 系统拓扑与参数的数据表示
我习惯用networkx存拓扑，用pandas存参数。networkx的Graph对象天然支持节点和边的属性字典，做连通性分析、寻找失电孤岛都很顺手。
import networkx as nx import pandas as pd G = nx.Graph() nodes_df = pd.DataFrame([ {"id": 0, "load_kw": 0, "type": "substation"}, {"id": 1, "load_kw": 100, "type": "load"}, {"id": 2, "load_kw": 120, "type": "load"}, ]) branches_df = pd.DataFrame([ {"from": 0, "to": 1, "r": 0.092, "x": 0.047, "cap_kw": 400, "lam_km": 0.1, "len_km": 0.6, "r_repair": 5.0}, {"from": 1, "to": 2, "r": 0.110, "x": 0.060, "cap_kw": 350, "lam_km": 0.1, "len_km": 0.4, "r_repair": 5.0}, ]) for _, nd in nodes_df.iterrows(): G.add_node(nd["id"], load_kw=nd["load_kw"], type=nd["type"]) for _, br in branches_df.iterrows(): G.add_edge(br["from"], br["to"], r=br["r"], x=br["x"], cap_kw=br["cap_kw"], lam=br["lam_km"] * br["len_km"], r_repair=br["r_repair"])线路的lam在存图之前就乘好了长度，这是我想强调的第二个关键细节：故障率参数从“每公里”换算成“整条线路”，否则后续所有计算都会出错。
5.2 非序贯蒙特卡洛可靠性评估骨架
import numpy as np def evaluate_reliability(G, load_nodes, mc_samples=5000, seed=42): rng = np.random.default_rng(seed) branches = list(G.edges) total_failure_loss = 0.0 for _ in range(mc_samples): # 抽样所有支路状态：1=正常，0=故障 avail = np.array([1 - G[u][v]["lam"] * G[u][v]["r_repair"] / 8760.0 for u, v in branches]) state = rng.random(len(branches)) < avail failed = [branches[i] for i in range(len(branches)) if not state[i]] if not failed: continue # 临时移除故障支路，分析失电负荷 G_tmp = G.copy() G_tmp.remove_edges_from(failed) loss = calculate_loss_of_load(G_tmp, load_nodes) total_failure_loss += loss return total_failure_loss / mc_samplescalculate_loss_of_load是核心函数，需要根据实际网络判断哪些负荷节点失电。我常用的做法是先找出所有从变电站母线能到达的节点，不在可达集合里的负荷节点就是失电节点；然后对有联络线的失电孤岛，逐一尝试合闸转供，转供成功的负荷从失电集合中扣除。
5.3 NSGA-II主循环骨架
快速非支配排序是所有多目标进化算法的地基。核心思路是：对每个个体p，维护两个集合——被p支配的个体集合S_p，支配p的个体数量n_p；n_p为0的个体属于第一前沿；然后逐层剥洋葱。
def fast_non_dominated_sort(fitness_values): pop_size = len(fitness_values) S = [[] for _ in range(pop_size)] n = [0] * pop_size front = [[]] for p in range(pop_size): for q in range(pop_size): if dominates(fitness_values[p], fitness_values[q]): S[p].append(q) elif dominates(fitness_values[q], fitness_values[p]): n[p] += 1 if n[p] == 0: front[0].append(p) i = 0 while front[i]: next_front = [] for p in front[i]: for q in S[p]: n[q] -= 1 if n[q] == 0: next_front.append(q) i += 1 front.append(next_front) return front[:-1]主循环就是标准结构：初始化种群、评估、非支配排序、拥挤距离计算、锦标赛选择、SBX交叉、多项式变异、父子合并、环境选择，迭代到最大代数。这段逻辑网上有很多开源实现，但要注意把前面说的“修复策略”插在交叉变异之后、评估之前，否则算法搜索空间会被大量非法个体污染。
5.4 并行评估与随机数隔离
蒙特卡洛加进化算法最大的痛点就是算力。Python里最简单的提速手段是用multiprocessing.Pool对种群个体做并行评估。这里有一个非常隐蔽的坑：如果所有worker共用同一个np.random.default_rng(seed)，在并行环境下很容易产生重复的抽样序列，导致不同个体的可靠性评估结果高度相关，Pareto前沿会失真。
正确的做法是按个体索引生成独立随机数流：
def eval_worker(individual, idx, total_plan_params): rng_seed = total_plan_params["base_seed"] + idx return evaluate_individual(individual, rng_seed, total_plan_params)每个个体拿到独立的种子，评估结果之间互不干扰，并且重新运行时可以复现。这个设计在后面排查“为什么跑两次结果不一样”的时候能救你一命。
6. 算例实验与结果解读：Pareto前沿上的方案怎么选
6.1 算例设置与参数清单
我用一个改进的IEEE 33节点配电系统做演示算例。基准网络节点0为变电站母线，系统峰值负荷约3.7兆瓦加2.3兆乏。候选方案设置如下表。
| 规划要素 | 候选设置 |
|---|---|
| 分布式光伏候选位置 | 节点11、18、23、29 |
| 光伏容量档位 | 0 / 0.2 / 0.4 / 0.6 MW |
| 储能候选位置 | 节点18、33 |
| 储能容量档位 | 0 / 0.2 / 0.4 MWh |
| 待升级线路 | 1-2、3-4、5-6、23-24 |
| 新增联络线候选 | 节点8-21、18-33、25-29 |
经济性参数：贴现率8%，光伏投资约150万元/兆瓦，储能投资约120万元/兆瓦时，线路升级每公里约40万元，运行维护费率光伏2%、储能3%，购电价0.55元/千瓦时。可靠性参数采用第3章表格里的参考值。仿真年限20年。
6.2 Pareto前沿长什么样
经过200代进化，跑出来的Pareto前沿大致呈现一个右上角向左下角弯曲的形态。横轴是年费用F1，纵轴是EENS。
| 方案 | 年费用（万元/年） | EENS（兆瓦时/年） | SAIDI（小时/用户·年） | 特征 |
|---|---|---|---|---|
| A 经济最优 | 238 | 18.6 | 1.82 | 少DG、无联络线，投资最省 |
| B 可靠性最优 | 372 | 3.1 | 0.31 | 联络线全覆盖，DG渗透率高 |
| C 拐点折中 | 301 | 6.2 | 0.62 | 曲线最陡处附近，边际效益最高 |
这三组数是演示用的示意结果，标准IEEE 33系统参数和实际负荷曲线不同，数值会变动，但形态规律是通用的。可以看出方案B比方案A每年多花134万元，换来EENS降低15.5兆瓦时/年；而方案C比方案A多花63万元，EENS降低12.4兆瓦时/年。也就是说拐点之后再多投钱，可靠性改善的效率会明显下降。
6.3 前沿解怎么落到最终方案
Pareto前沿给出的是一组互不支配的方案，最终选哪个取决于决策者的偏好。工程上我常用两步走：先根据项目预算上限把前沿截一段，排除那些年费用超预算的方案；然后对剩余方案做多准则决策，比如用TOPSIS或熵权法，把年费用、EENS、SAIDI、SAIFI这几个维度综合打分。
一个容易犯的错是拿一个“权重”把双目标直接压成单目标再排序。如果预算没有硬约束，权重不同选出来的点可能完全不同，报告里说不清楚。正确的做法是先把前沿摆出来，让决策者看到“多花钱买可靠性”的边际曲线，再结合预算和停电损失心里价位来选点。
7. 实测中踩过的坑：可靠性评估结果“复现不一致”的五个根源
7.1 随机数种子与并行重复抽样
这是我在多个项目里见过最多的问题。代码写在Jupyter Notebook里串行跑，结果没问题；一上服务器用multiprocessing.Pool，前后两次结果就对不上。原因几乎都是worker里的随机数生成器共享了同一个初始种子，或者每个worker在复制父进程全局变量时继承了相同的np.random状态。
解决办法就是前面说过的：按个体索引生成独立种子。这个改动的效果立竿见影，不用怀疑。
7.2 可靠性指标的口径不一致
同一个EENS，有人按负荷点年电量统计，有人按用户数加权，有人把故障后的瞬态欠压也算进去，结果能差出好几倍。同行之间对结论，第一件事不是看代码，而是看指标定义。我建议在报告里明确写出：EENS按负荷点有功乘以年停电时间累加，SAIDI按用户数加权，SAIFI按用户数加权，停电时间不含电压暂降。
7.3 负荷与DG时序曲线缺失导致的高估
如果只用年峰值负荷跑非序贯蒙特卡洛，EENS会被严重高估。道理很简单：故障发生在一个深夜低负荷时段，实际失电量比峰值负荷下小得多，但恒定峰值负荷模型会把这段停电按峰值计算。工程上至少要做季节典型日或者等值小时曲线，把8760小时压缩成几十个代表性时段，再和非序贯抽样交叉计算。算例里的数字差别动辄30%到50%，这个坑不填，后面的方案比选结论可能直接反转。
7.4 修复时间分布不能一律用指数分布
指数分布的特点是无记忆性，均值是1/λ，但修复时间实测数据往往是右偏的，偶尔出现极长的故障修复时间，对EENS的高尾段影响很大。实测中用威布尔分布拟合修复时间，尾部更厚重，更贴近真实停电事件。如果你接手过工程现场的历史故障台账，你会对这个深有感触：大部分的故障是小问题快速修复，但少数严重故障会拖几天，这少数事件往往贡献了大部分EENS。
7.5 进化算法嵌套蒙特卡洛带来的“伪收敛”
NSGA-II每一代都在用随机评估结果比较个体优劣。如果评估的蒙特卡洛抽样次数太少，个体之间的排序噪声会非常大，算法看起来收敛了，实际上只是“噪声被优化掉了”而不是“真实目标被优化了”。要验证是不是伪收敛，可以对最终前沿的每个解，用5万到10万次抽样重新评估，再看前沿是否明显移位。如果移位很大，说明进化过程中评估精度不够，需要加大抽样次数或者用公共随机数（同一组随机数评估所有解）来降低比较噪声。
最后再分享一个我自己的调试习惯。每次写可靠性评估模块，我不会直接拿IEEE 33节点全系统跑，而是先构造一个三节点两线路的微型系统，手动枚举所有单故障场景，计算出理论EENS，再拿蒙特卡洛程序去对答案。对上了再扩到完整算例。这个“故障注入式单元测试”看似多花了半小时，但能省下后面优化跑完后才发现评估代码写错的数天返工时间。双目标规划项目里，算法永远不是瓶颈，把可靠性评估写对才是真正的分水岭。
