半年前我接手一个县级配电网的扩建规划,用的还是老一套最小投资模型,结果方案被业主方退回来三次。头一次只压投资,电网公司问“停电时间多少”;第二次加了个N-1约束,把所有走廊都按最大型号扩容,预算超了百分之三十;直到我把经济性和可靠性拆成两个目标,在Pareto前沿上挑出一组方案,才开始被认可。以我自己的经验,混合配电系统规划从来不是一个“算出最优解”的问题,而是一个“在风险和成本之间找出可接受解集”的问题。这篇文章就聊一聊基于经济与可靠性双目标的混合配电系统规划及可靠性评估思路,以及我整个Python代码实现过程中踩过的坑、验证过的方法、积累的数据上的体会。适合正在做配网规划项目的工程师、写相关论文的研究生,以及想用Python把电网拓扑优化和可靠性仿真串起来的人。
1. 混合配电系统规划:为什么“省钱”和“可靠”必须放在同一个模型里
1.1 所谓“混合”到底混了什么
很多同行一看到“混合配电系统”这个说法,第一反应是“交流+直流”。实际工程中,这个词的覆盖面宽得多。我接触的项目里至少见过三种形态:传统辐射网接入分布式光伏、风电和储能的多馈线网络;含联络开关、能被上游变电站转供电的“辐射+多联络”复合网架;还有通过电压源型换流器把直流线路嵌入交流配电网的交直流混合系统。这几种形态的共同点是:供电拓扑不再单一,电源不再只来自变电站母线,运行方式存在大量不确定性。因此,规划的约束和评估维度都要比传统单辐射配电网复杂很多。论文或项目标题里写“混合配电系统规划”,大多数情况下就是指这类含分布式电源、储能和联络通道的复合配电网规划。
1.2 单目标规划的盲区:最便宜的方案往往最贵
我刚做这类项目的时候,第一版模型目标函数只有“年度综合费用最小”,包括线路建设、DG投资、运维和网损。优化出来的方案确实便宜,但它把新增负荷全部压在了少数几条主干线上,末端电压贴近下限,任意一条中压馈线故障,下游一大片负荷停电。整个方案在“正常运行”工况下看起来没问题,一旦遇到故障或检修,风险就暴露得干干净净。
这让我意识到一个关键点:配电系统规划和输电网规划有本质区别。配网故障离用户更近,供电可靠性直接对应到用户侧的停电时户数。省下来的一段电缆钱,在停电损失里可能翻几倍还回去。这就是单目标规划最典型的盲区——它把“未来可能发生的故障事件”默认成了零成本。实际上,从全寿命周期角度看,一个可靠性的微小改善,往往比一次重大停电的经济价值更可观。
1.3 双目标优化的实质:不是找唯一最优,而是给出决策序列
为什么双目标比单目标更贴合工程实际?因为经济性和可靠性之间天然存在对抗关系。你多花一千万新建联络线,可靠性指标确实会变好,但经济性变差;你不花这笔钱,指标反过来。这种问题不存在一个解在所有维度上都碾压其他解,只存在一个Pareto最优解集:在解集中,改善经济性必然牺牲可靠性,改善可靠性必然增加成本。
工程师真正的价值恰恰在于把这一整条“成本—可靠性”曲线交给决策者,让决策者根据投资上限、可靠性考核指标和负荷重要性去选点。这个思路几乎适用于所有配电网扩建规划项目,也特别适合写成论文——目标清楚、方法标准、结果可以可视化。它的输出不再是一个武断的“最优解”,而是一组有工程解释的候选方案,决策过程也变得透明。
2. 双目标数学模型:决策变量、目标函数与约束条件的工程化设计
要把问题交给Python,第一步就是把工程语言转成数学语言。这里有两个原则:一是决策变量必须对应真实工程操作,二是目标函数能算、可比、可解释。
2.1 决策变量怎么定义
规划问题里,决策变量就是规划人员要“决定”的东西。我一般把这类项目的决策变量分成三类:
| 决策变量 | 含义 | 类型 |
|---|---|---|
| x_l ∈ {0,1} | 是否新建/升级馈线线路 l,含路径选择、导线型号 | 0-1整数 |
| P_PV,i、Q_PV,i | 节点 i 接入光伏的额定容量 | 连续 |
| P_ESS,i、E_ESS,i | 节点 i 储能额定功率和额定容量 | 连续 |
| y_DG,i ∈ {0,1} | 是否在节点 i 建设分布式电源 | 0-1整数 |
这里有一个经验:分布式电源和储能的“是否建设”与“建设多少”要分开建模。原因很简单,混合整数非线性规划的求解难度与0-1变量数量高度相关。如果每个节点只能装固定容量,那直接一个0-1决策就够了;如果需要选容量等级,就得再加一组离散档位变量,解空间瞬间膨胀。我通常先做灵敏度分析,确定哪些节点值得装,再把这些节点的容量设成连续变量,能省掉大量无效解。
另外还要注意变量之间的耦合关系。比如储能只能建在有DG接入的节点时,就需要加逻辑约束 y_ESS,i ≤ y_DG,i,否则优化器会在一个无电源节点凭空放一套储能,结果无法工程实现。
2.2 经济性目标函数:把钱算清楚
经济目标建议用“全寿命周期等年值费用”而不是一次性总投资。原因在于DG和线路的寿命不同:光伏25年、电缆30年、储能可能只有10年。如果直接把总投资相加,模型会天然偏好短期便宜、寿命短的设备。我采用等年值法:
f_econ = r_a × (C_inv + C_OM) + C_loss
其中C_inv是折现到建设年的总投资,C_OM是年运维成本,C_loss是网损年费用,r_a是资金回收系数。网损年费用按下式计算:
C_loss = 8760 × c_elec × P_loss_avg
这里的P_loss_avg需要用典型年曲线的时序潮流算出来,不能用一个峰值负荷代替。我第一版建模时就偷懒用了最大负荷时段的网损,结果经济目标被严重高估,优化方案里网损排序全乱,后来换成典型日曲线取加权平均才正常。
关于停电损失C_ENS要不要放进经济目标,这里有个容易踩的坑。如果把它放进经济目标,实际上就是一种“经济化可靠性”的单目标思路,把可靠性按停电损失单价折算成钱;如果同时又把ENS单独当作第二目标,两个目标就不是真正独立的,等于暗中给两个目标绑定了线性权重,Pareto前沿会扭曲。我的做法是:优化阶段经济目标只含投资/运维/网损,可靠性目标用ENS,最后在方案比选阶段再用停电损失费用做一次整体评估。这样模型干净,解释也清楚。
2.3 可靠性目标函数:为什么我选ENS
可靠性指标有SAIFI、SAIDI、CAIDI、ENS等,都能刻画“可靠程度”。在多目标优化里,我建议把ENS(期望缺供电量)作为第二目标函数。原因有三点:第一,ENS是能量量纲,单位是MWh/年,在数值上与投资费用可比,画Pareto图时两个坐标轴量纲不会差出好几个数量级;第二,ENS天然与停电损失费用线性相关,后续做经济评估时衔接很顺;第三,ENS能直接反映“故障期间有多少负荷被削减”,对电网公司来说,这个数字比停电次数更好理解。
ENS并不能完全替代SAIFI和SAIDI,因为它们刻画维度不同。我的做法是:优化阶段以ENS为第二目标函数,方案输出后,再用SAIFI、SAIDI做一次离线核算,用来填电网公司的可靠性考核表。这样在优化器里模型简洁,在工程报告里指标全面。
2.4 约束条件里最容易写漏的三个东西
建模时约束条件决定优化结果是否可信。我吃过亏的地方有三个:
第一是潮流平衡约束。配电网常用DistFlow模型描述有功、无功注入与支路流量的关系,以及电压降方程。在优化阶段我坚决用线性化DistFlow,不要直接上完整交流潮流,否则求解时间完全失控。线性化后的精度在规划阶段完全够用,最终方案的精确校验再交给完整潮流工具。
第二是辐射状结构约束。正常运行状态配电网要求开环运行,拓扑必须是无环连通图。数学上常用 |E| = |V| - 1 加连通性约束表达。但注意,这个约束拆开写才有效:只写边数关系不写连通,优化器会构造出“有环但有孤岛”的图;只写连通不写边数,又可能自动闭合联络线。两点必须同时成立。
第三是N-1安全约束。任意一条线路故障退出后,系统仍能通过转供保证不失负荷,或者切除负荷在允许范围内。这个约束最严格,也最消耗计算资源,通常不会原样写成显式约束,而是在评估阶段对每个候选方案做N-1校验,违反的用罚函数或舍弃处理。
还有一个细节经常被忽略:电压约束不能只看稳态下限。DG出力高峰且负荷低谷时,可能出现功率倒送,电压不是偏低而是偏高。我初版模型只加了电压≥0.95 p.u.的下限约束,跑出来的方案在仿真中被巡检出大量过电压节点。现在我的模型一律写作 0.95 ≤ V_i ≤ 1.05 p.u.,并且用风光出力典型场景分场景校验。
3. 可靠性评估:从指标到序贯蒙特卡洛模拟
模型建好之后,下一步就是怎么算出ENS。这一步是整个项目技术含量最高的地方,也最考验工程判断。
3.1 四个核心指标先厘清
| 指标 | 定义 | 公式含义 |
|---|---|---|
| SAIFI | 系统平均停电频率 | 总停电用户次数 / 总用户数 |
| SAIDI | 系统平均停电持续时间 | 总停电用户小时 / 总用户数 |
| CAIDI | 用户平均停电持续时间 | SAIDI / SAIFI |
| ENS | 期望缺供电量 | 所有故障场景下损失电量期望 |
规划阶段我主要盯ENS,校验阶段再核算SAIFI/SAIDI。这里要特别提醒:SAIFI和SAIDI的计算需要知道每条馈线上的用户数。有些项目资料只有负荷功率没有用户数,此时只能用ENS或AENS做目标,避免口径不一致。
3.2 解析法为什么在混合配电系统里不够用
最早的可靠性评估常用故障枚举法:把每条线路、每台变压器的故障事件列出来,算出故障概率、负荷削减量,再求和。这个方法在单一电源、单一辐射网络里非常好用,因为故障影响范围清晰。可一旦网架含DG、储能、联络开关,情况就变了:储能充放电行为与故障发生时刻强相关,光伏出力随时间变化,联络转供需要设备动作时间,这些已经不再是“事件概率相乘”能准确刻画的。解析法状态数量随元件数指数膨胀,几十条线路就够呛,更别提时序相关性。
3.3 序贯蒙特卡洛模拟的正确打开方式
我实际采用的是序贯蒙特卡洛模拟。核心思想是:在时间轴上模拟每个元件的“正常—故障—修复—正常”循环,让系统逐年、逐小时演化,统计停电事件,最后用多年的平均结果估计可靠性指标。
具体步骤如下:
- 设置仿真时钟,按小时步长推进,一般按8760小时/年。
- 对每个元件(馈线段、变压器、断路器)抽取首次故障时刻TTF。
- 找出最早故障元件,将时钟推进到该时刻。
- 判断故障影响:隔离故障区段后,哪些负荷可以被转供,哪些负荷必须切除。
- 若负荷被切除,累计停电时户数和缺供电量;同时抽样该元件修复时间TTR。
- 该元件修复后网络恢复,循环回到步骤2,直到达到仿真年限。
这里最花时间的是步骤4的“影响分析”。在只有十几条支路的小算例里,每步甚至可以用networkx的连通性判断完成;到了真实规模,就必须对每个时间步做潮流和电压约束校验,性能会急剧下降。我在第5节会专门讲提速方案。
Python侧的采样代码本身很简单,核心就几行:
import numpy as np def sample_ttf(failure_rate_per_year, t_now): # 指数分布:元件平均无故障时间 = 1/λ(年) # 输入λ单位是次/年,输出为小时 return t_now + np.random.exponential(365 * 24 / failure_rate_per_year) def sample_ttr(mttr_hours): # 修复时间分布,工程上常用指数分布近似 return np.random.exponential(mttr_hours)真正复杂的是“怎么把故障事件连成时间流”和“故障后的负荷削减逻辑”,这两部分必须在网络拓扑数据结构上做文章,而不是靠简单循环。
设备切换动作时间在月级长周期模拟中通常小于1小时,所以中压配网年可靠性分析按小时步长已经足够。只有做秒级短时可靠性分析时,才需要细化联络开关动作时间。
3.4 负荷削减判断的两种策略
故障后的负荷削减量,取决于有没有冗余通道。我的判断分两层:
- 快速判断:先用networkx判断故障后网络是否依旧连通。如果故障元件两端存在备用通道,整个网络仍然连通,则可以初步认为不需要切负荷(在容量允许前提下)。这一步是O(|V|+|E|)级别,非常快。
- 精确判断:如果网络虽连通但某些支路过载或电压越限,需要做一次线性DistFlow,找到需要削减的量。这一步较慢,只对少数严重场景执行。
这个分层判断能筛掉80%以上的非故障影响场景,把蒙特卡洛时间从“不可用”降到“可用”。
4. Python实现全链路:从拓扑建模到NSGA-II求解与Pareto前沿分析
下面进入代码实现的主线。我把整个项目拆成四步:拓扑建模、评估函数、优化求解、结果输出。环境依赖很简单,pip install numpy networkx matplotlib pymoo就够了。
4.1 用NetworkX构建配电网拓扑
我最开始用自定义字典存节点和线路,代码写起来灵活,但后期可靠性模拟处处要查连通性、算路径、找邻居,把字典从头实现一遍太痛苦。后来全部改用NetworkX,工程效率高一个数量级。
构建拓扑的骨架代码:
import networkx as nx import numpy as np G = nx.Graph() # 节点:0为变电站,其余为负荷或分布式电源节点 G.add_node(0, node_type='source', v_rated=10.0) # 10 kV中压配网 G.add_node(1, node_type='load', peak_p=2.0) # 峰值负荷2MW G.add_node(2, node_type='load', peak_p=1.5) G.add_node(3, node_type='pv', cap_p=0.8) # 光伏节点 # 线路边:属性包括长度、单位阻抗、容量、故障率 G.add_edge(0, 1, length=2.1, r=0.13, x=0.08, cap=6.0, fail_rate=0.08, mttr=5.0) G.add_edge(1, 2, length=1.5, r=0.13, x=0.08, cap=4.5, fail_rate=0.06, mttr=5.0) G.add_edge(2, 3, length=1.8, r=0.13, x=0.08, cap=4.0, fail_rate=0.07, mttr=5.0)这里fail_rate的单位到底是“次/年·公里”还是“次/年·条”,非常容易搞混。我的经验是:每条线路的故障率=单位长度故障率×线路长度(km),在代码注释里写清楚,不然换算例后指标会整体漂移,且很难察觉。
邻接矩阵在论文里有需求时可用nx.to_numpy_array(G)一行生成。但注意,这个矩阵只保留拓扑连通信息,不含电气参数。真正的电气计算我建议用边列表加属性的DataFrame,处理起来比矩阵直观得多。
4.2 优化里不能天天跑全潮流:线性化DistFlow
配网潮流精确计算有Pandapower这类工具,但放在NSGA-II里每个个体跑一年8760小时潮流根本不现实。我在优化阶段用的是线性化DistFlow,忽略损耗和电压降的高阶项:
P_j = P_i - p_load,j + p_dg,j - p_ess,j
V_j = V_i - (r_ij × P_ij + x_ij × Q_ij) / V_0
线性化之后,全网电压可以直接用一次矩阵运算序贯推出来:
def linear_distflow(G, p_inject, q_inject, v0=1.0): # p_inject:节点净注入有功(负荷取负、DG取正) # 仅适用于单源辐射网,通过后序节点逐个推算 from networkx.algorithms.traversal import dfs_postorder_nodes order = list(dfs_postorder_nodes(G, source=0)) v = {0: v0} for j in order: if j == 0: continue # 因为是辐射网,j的上一级节点唯一 i = next(iter(G.neighbors(j))) v[j] = v[i] - (G[i][j]['r'] * p_inject[j] + G[i][j]['x'] * q_inject[j]) / v0 return v这只是一个简化示意。真正实现要注意后序节点的负荷聚合,而且对于含联络开关、存在多馈线的情况,辐射状假设不成立,这时要退一步用直流潮流近似做N-1校验,不必解完整交流潮流。
4.3 NSGA-II求解:pymoo库的接入方式与骨架代码
双目标优化的求解,我推荐用pymoo里的NSGA2,这是当前Python生态里最省事的多目标进化算法库,Pareto前沿的可视化也完整。核心思路是自定义一个Problem类,把前两节定义的评估函数放进去:
from pymoo.core.problem import Problem from pymoo.algorithms.moo.nsga2 import NSGA2 from pymoo.optimize import minimize class DNPPlanning(Problem): def __init__(self, n_var): super().__init__(n_var=n_var, n_obj=2, xl=0, xu=1) # 决策变量归一化 def _evaluate(self, X, out, *args, **kwargs): f1 = np.array([self.calc_economic(x) for x in X]) f2 = np.array([self.calc_ens(x) for x in X]) out["F"] = np.column_stack([f1, f2]) algorithm = NSGA2(pop_size=40, n_offsprings=20) res = minimize(DNPPlanning(n_var=12), algorithm, ("n_gen", 50), seed=42, save_history=True)calc_economic和calc_ens就是前两节定义的评估函数。这里有一个性能关键点:_evaluate里不要一个一个循环调用蒙特卡洛,尽量向量化,或者用多进程并行。我把评估函数设计成传入整个种群,内部用multiprocessing.Pool分摊,提速非常明显。
4.4 Pareto前沿绘制与方案优选
优化完成后,把结果画出来。横轴是经济性f1(万元/年),纵轴是可靠性ENS(MWh/年),每个点是一个规划方案。正常情况下散点会呈现一条左低右高或左高右低的下降带,这就是Pareto前沿:
import matplotlib.pyplot as plt F = res.F plt.scatter(F[:, 0], F[:, 1], c='steelblue', alpha=0.7, s=40) plt.xlabel('Annual total cost (10k CNY/yr)') plt.ylabel('ENS (MWh/yr)') plt.grid(alpha=0.3) plt.tight_layout() plt.savefig('pareto_front.png', dpi=200)方案优选我一般用“拐点法”:找前沿斜率变化最剧烈的位置,拐点附近的解通常性价比最高——再增加一单位投资,可靠性改善已经明显放缓;再往后走,性价比持续下降。如果项目方有硬性投资上限,直接取前沿上不超过预算的最左侧点;如果有可靠性考核指标,就取满足考核的最小成本点。
5. 实测中的坑:可靠性仿真时间失控、收敛性与结果可信度
代码能跑通只是第一步,结果能被项目组信任才算完。我在这套代码上反复调试,踩过的坑集中在三个方向。
5.1 蒙特卡洛仿真慢到怀疑人生
第一版代码,NSGA-II迭代50代、种群40,总共2000个候选解;每个候选解用500年序贯蒙特卡洛评估,年步长8760小时。算完一个实验要三天。这不叫优化,这叫测试耐心。
后来做了四个改动,时间从72小时压到4小时左右。
第一,把评估从“逐年逐小时扫描”改成“事件驱动”。序贯蒙特卡洛不需要每个小时都评估系统状态,只在元件状态变化的时间点做网络分析。大量时间步里系统拓扑根本没变,重复计算毫无意义。
第二,故障影响分析用“连通性预筛+潮流精算”两级判断,快速排除无影响的故障事件。
第三,对蒙特卡洛抽样精度做动态控制。不是固定500年,而是每模拟200年就检查一次ENS的方差系数,降到阈值就提前停止。可以用简单的收敛判据:
β = σ / (μ × √N)
当β < 2%且ENS变化率小于1%时停止,既能保证精度又省算力。
第四,多进程并行。pymoo把候选解分给多个核,每个核独立跑自己的评估,这是理想并行场景,加速比接近核心数。
5.2 辐射状约束没写死,N-1校验全是违例
第二次修改模型时,为了让求解器跑得快,我把辐射状约束里的“连通性”条件漏了,只保留了 |E| = |V| - 1。结果优化器很聪明地构造出一堆“有环但边数刚好”的图,拓扑看起来满足条件,实际存在孤岛。这些孤岛在N-1校验里全部暴雷。
解决方式:不要在目标函数里软性惩罚,要写成硬约束,或者干脆在初始化种群时生成合法的辐射状拓扑,变异操作限定在“断一边、连一边”的环交换范围内。这样种群从初始阶段就保证连通,后期不用补救。生成辐射状连通拓扑的经典方法很朴素:先随机生成一个连通图,然后用“破圈法”反复删除能形成环的边,直到边数等于节点数减一。这样做出来的拓扑天然满足开环运行约束,后续优化只需在合法拓扑间做小幅变异。
5.3 指标波动大:收敛判据与参数兜底
早期遇到一个诡异现象:同一个算例跑两次,ENS结果差了8%。一开始怀疑随机种子问题,后来发现是收敛条件设置得过于宽松。改进后,我使用5.1节提到的方差系数作为停止条件,并要求ENS滑动平均变化率小于1%,结果重复运行基本稳定在2%以内。
另一个参数敏感性问题是可靠性基础数据。线路故障率、修复时间、联络开关动作时间,这些参数在论文里可能是假设值,但在工程里必须来自运行统计数据。我建议把所有可靠性参数单独放在一个config.yaml文件里,方便换算例、做敏感性分析:
reliability: line_failure_per_km_year: 0.08 transformer_failure_per_year: 0.02 mttr_hours: 5.0 switch_operation_hours: 0.5计算完成后还要做校验:把特定故障场景下的ENS与手工推演的结果比对,或者与历史停电统计对照,确认模型行为符合实际。
6. 最后对想做这类项目的人说几句心里话
写到这里,主体内容基本讲完了。我个人在这类项目里最大的体会是:模型复杂度永远要跟决策精度匹配。给电网公司做规划,不要把算法堆到他们不敢用、甚至看不懂的程度。先把单辐射网络的潮流、可靠性串起来,再逐步加DG、储能、联络转供,每加一层都做一次结果验证。Python这个工具链最大的优势不是哪一步有多高级,而是你能把“网络建模—优化求解—可靠性仿真—出图报告”全部打通在一个环境里,不用在多个软件之间导数据,改一个参数重跑全流程的成本也低得多。
还有一个容易被忽略的交付细节:最后给业主的图,一定要有原始负荷分布、方案网架结构、DG/储能位置、故障后转供路径这几张图;如果能把Pareto前沿图和一个对比表格同时给出,把经济指标、ENS、SAIFI、SAIDI全部列在一起,逻辑闭环就非常完整,评审时省很多口舌。
如果你正在写类似论文或者刚接手类似项目,建议先把基础算例跑通,再逐步加复杂度。遇到性能瓶颈,优先优化评估函数而不是优化器;遇到指标对不上,先查元件故障率单位,再查网络连通性。这个项目最终交付时,一组经过校核的可靠参数和一版可读性良好的代码,比任何花哨的算法都重要。