简介这份资料是面向电力系统规划与运行研究者的论文复现资源针对配电网面临台风等极端灾害时的全过程韧性提升问题给出基于Wasserstein距离模糊集与分布鲁棒机会约束的完整建模与求解实现适合关注韧性提升、分布式鲁棒优化的学者与工程师。压缩包仅含1个docx文档大小56KB文档内整合了论文概要、数学模型推导和可运行的Python代码及逐段注释覆盖“预防-应急-抢修”双层三阶段规划、CVaR转化、混合整数二阶锥规划求解等内容。资源发布以来已有132人学习浏览。通过这份文档读者可以快速理解如何构建不确定性模糊集掌握台风风场建模、线路故障概率建模及韧性指标量化的具体方法并能对照代码复现仿真实验用于自身课题或工程项目中的配电网防灾规划参考。1. 配电网全过程韧性和分布鲁棒机会约束规划一次模型选择上的止损在10kV配电网的短路电流仿真里一旦叠入台风场景负荷和分布式电源的历史分布会瞬间失效之前的规划结果短时间内出现系统性偏差。常见的机会约束规划把经验分布当真实分布算出来的韧性方案往往在极端灾害下“翻车”。分布鲁棒机会约束规划不这么做它以经验分布为中心用一个概率球罩住所有可能的真实分布要求即使在球内最差的那个分布下支撑韧性目标的可靠性约束也必须满足。这使得灾前预布局、灾中响应和灾后恢复全过程都能拿到可证的下界。这套方法适合负责配网韧性规划、微网调度和并网方案评估的工程师。下面按一个能落地的求解链路走模型原理、多类型资源建模、完整代码以及调参避坑。2. 分布鲁棒机会约束规划Wasserstein球半径与机会约束的等价转化2.1 为什么普通机会约束会在灾害场景失效普通机会约束的条件写成 P( g(x, ξ) ≤ 0 ) ≥ 1 − ε。这里的P是历史样本拟合出来的经验分布ε是允许的风险概率。这个表达式的吸引力在于直观它把不确定性从“必须逐个场景满足”的确定性鲁棒里释放出来允许少量极端场景不满足。在配电网规划里典型的g(x, ξ)写成功率供给与负荷需求的差x是新建储能和柴发容量ξ是净负荷。这个约束可以算出容量配置只要把某个区间内的缺口控制在5%以内就算过关。问题出在“灾害场景下的P和建模时用的P不是同一个P”。台风、冰灾这类小概率高影响事件在正常月份几乎没有样本而在事件发生的当天负荷和光伏出力的条件分布会整体偏移。用日常样本估计的均值、方差放到灾害窗口里常常是低估计的。这时候普通机会约束允许的5%缺口在真实分布下可能放大到15%甚至更高。配电网的灾前预布局、灾中响应和灾后恢复这三个阶段共用同一套容量一旦灾中的缺口超标后期的恢复曲线也会整个右移。这也是为什么“全过程韧性”必须把分布偏差算进去而不是只关心期望成本。用黑匣子式的置信度拍脑袋一定会在第一轮极端工况下翻车。解决思路是把机会约束里的“那个P”从一个固定分布放宽成一个分布集合。设经验分布为P_hat真实分布虽然未知但可以假定它落在以P_hat为中心、半径为r的分布球里。球内的任意一个分布都可以是真实分布机会约束必须对球内最差的那个分布也成立。这个“最差分布”的构造正是分布鲁棒机会约束规划和普通机会约束的区别前者先和最难缠的分布算账再让你决定配置。2.2 Wasserstein距离把“分布不确定”收进一个球分布球用Wasserstein距离来定义。它回答“把一个分布移动成另一个分布最小成本是多少”在配电场景里直观理解是移动负荷和光伏样本值让一个概率密度变换成另一个需要付出的期望“路程”。数学上对于两个离散场景分布q和pWasserstein距离可以用线性规划求解而在工程上我习惯直接调用scipy的wasserstein_distance做一维粗略估计。半径r是这个球的核心参数。r趋向0时分布集合收缩到经验分布模型退化成普通机会约束r增大模型被迫考虑越来越偏的分布解会变得保守。半径不是越小越好也不是越大越好。常见的经验估式是 r ≈ C × sqrt( log(1/β) / (2N) )其中N是场景样本数β是希望构造的置信度比如0.05C是样本空间直径的估计量——在10kV配电网里C可以用短路电流仿真时得到的净负荷极差来标定。我记得一次台风算例中净负荷波动从0.2到2.8兆瓦C就取了2.6兆瓦这样r的量级大约是0.2兆瓦看起来不多但足以把5%的违背率顶到8%。实际做工程时我不建议直接把这个r定死而是把它放到验证集上做一维扫描r从小到大观察最坏违背率的变化曲线选在曲线拐角处再结合投资预算做权衡。这一步没什么捷径可走跟调PID参数一样得靠几次重跑找感觉。提示Wasserstein球半径的单位和约束变量一致。如果你把负荷写成标幺值r也必须落在标幺值系里否则后面所有对偶转化和数值求解都会漂移。2.3 半径怎么拍置信度、场景数和碳足迹之间的博弈场景数N对半径的影响遵循 sqrt(1/N)。N从50加到200r大致缩小到一半但场景数增多会让每次评估的计算成本线性上涨。在阶段规划中我一般先跑N50的快速扫描快速找到容量下限再用N200到500做精算。样本抽法这里也要注意直接蒙特卡洛采样会让尾部稀疏我会先用拉丁超立方生成若干个分位场景再把历史极值与合成场景混合形成带“厚尾”的样本集。下面这段代码用来估算两个场景分布之间的Wasserstein距离并给出r的初始候选值。这里的base是历史净负荷标幺序列storm是叠加台风扰动后的再采样序列。import numpy as np from scipy.stats import wasserstein_distance rng np.random.default_rng(7) base rng.normal(1.0, 0.15, 200) # 历史净负荷标幺 storm base * rng.uniform(0.8, 1.3, 200) # 台风扰动后的净负荷 wd wasserstein_distance(base, storm) N len(base) beta 0.05 C max(base.max() - base.min(), storm.max() - storm.min()) r_guess C * np.sqrt(np.log(1.0 / beta) / (2 * N)) print(wasserstein dist:, wd) print(r_guess:, r_guess)这里beta取0.05表示希望以95%的信心认为真实分布落在半径球内C是样本极差这两个参数在配网工程里都会输入到后面的规划模型中。wd本身是经验分布间的最优输运距离如果wd超过r_guess说明台风扰动已经把分布推到规划考虑范围外就要上调C或者补充样本。这段代码并不直接求解规划却决定后续所有机会约束的“口径”踩坑最多的地方也在这里。对偶转化是分布鲁棒机会约束的核心引擎。当g(x,ξ)关于ξ是仿射函数时球内最坏分布下的机会约束可以被改写成一组带对偶变量的线性约束复杂度跟原来的场景规模同阶。这个转化在代码里看起来是一段求解器的约束构建但背后涉及的条件风险值近似并不是每个工程师都会手动推一遍。所以我的经验是先用重采样法确认方案趋势再把关键的仿射机会约束交给成熟求解器去处理不要把对偶推导硬编码在自己的脚本里。处理不好数值波动会比想象的大很多。3. 多类型韧性资源建模分布式电源、储能与柔性负荷的约束写法3.1 韧性阶段资源匹配谁在灾前、灾中、灾后起作用“全过程韧性”不是单点指标它至少包含四个阶段灾前防御、灾中响应、灾后恢复以及长期修复。不同阶段对应的资源动作完全不同。灾前主要是储能充电和柴发热备把SOC拉到上限灾中储能和柴发共同承担缺额柔性负荷开始按优先级切载灾后恢复阶段光伏若复归还要考虑储能再充电和负荷逐步并入。一个约束模型如果不区分阶段把所有资源塞进一个静态不等式解出来往往是“全能但单薄”的答案容量不少关键时刻哪个都顶不上。所以建模的第一步是给资源打上阶段标签。分布式电源里光伏在台风期间出力断崖式下跌只能当“后援”柴发可以作为灾中主力但有启动时间和爬坡限制。储能是少数能同时覆盖灾前和灾中的资源因为它的功率双向可调柔性负荷则提供可短时中断的备用本质上是把缺电从“系统硬缺口”转移成“用户软降额”。电动汽车充放电也可以承担负荷侧调节但调度优先级通常低于专用储能。下面这张表格给出我常用的一组资源-约束对照关系。注意参数只是示意具体取值要按10kV馈线容量和可靠性要求调整。| 资源类型 | 参与阶段 | 主要决策变量 | 核心约束 | 成本类型 | | 储能 | 灾前、灾中、灾后 | 充放电功率P_ch/P_dis容量E | SOC递推、C-rate限制、容量上下限 | 容量年化成本 | | 柴发 | 灾中、灾后 | 出力P_gen、启停状态y | 爬坡率、最小启停时间、出力上下限 | 燃料启停成本 | | 光伏 | 灾后复归为辅 | 输出P_pv(不确定) | 时序出力上界、逆变器容量 | 固定投资 | | 柔性负荷 | 灾中 | 切负荷量shed | 切负荷比例上限、最小恢复时间 | 缺电补偿 | | 电动汽车 | 灾中 | 放电功率P_ev | 电池约束、并网点容量 | 放电补偿 |这张表在后面代码里都会转化成结构化约束而不是手感式的惩罚项。把每个资源独立建模再在目标函数里合成是避免出现“储能容量很大但功率上不去”这类尴尬局面的前提。3.2 储能和分布式电源的时序约束储能的时序约束最关键因为它的当前电量是上一时刻电量和充放电过程的积分结果。写成离散时间步时就是状态变量在24个时段的递推式。常见写法如下import cvxpy as cp T 24 E cp.Variable(T) # SOC单位kWh p_dis cp.Variable(T) # 放电功率单位kW p_ch cp.Variable(T) # 充电功率单位kW E_max, E_min 800.0, 200.0 P_dis_max 400.0 P_ch_max 300.0 eta_c, eta_d 0.95, 0.90 soc_constr [ E[0] 600.0, # 灾前初始电量 E[t1] E[t] eta_c * p_ch[t] - p_dis[t] / eta_d for t in range(T-1) ] power_constr [ E E_min, E E_max, p_dis 0, p_ch 0, p_dis P_dis_max, p_ch P_ch_max ] constr soc_constr power_constr代码里的E_t1表达式把充放电效率放到能量递推里是一种简化处理。严格的模型还会把“同时充放”的二进制变量加进去避免同一时段既充电又放电的虚增收益在投资规划阶段我会先不加二进制变量用线性递推快速跑出容量区间等进入调度验证时再加上。P_dis_max和P_ch_max分别限定了储能功率上限E_max限定了容量两者之间的比值就是C-rate比如800kWh配400kW放电刚好是0.5C对铅炭和磷酸铁锂都算合理。柴发这类DG的约束不止出力上界还有爬坡率和最小运行时间。在灾后恢复阶段如果系统负荷快速爬升柴发爬坡跟不上就会产生二次缺电。表述为相邻时段出力差限幅运行状态变量y_t如果从0变1需要持续运行至少M个时段才能再停。这类整数约束投资规划常常省略但如果你做的是小时级全过程仿真还是保留为好否则恢复阶段会出现“开了关、关了开”的奇怪解。3.3 柔性负荷与网络拓扑的耦合约束柔性负荷的实质是切负荷决策。它的约束通常会写成一个可切比例上限同时为了保证基本供电还要设置一个最低可切总量。这里有个工程取舍切负荷成本往往只有可靠供电压电价的一半优化器会倾向于把它用满所以在目标函数里要给切负荷一个较高的机会成本否则韧性约束形同虚设。网络耦合是另一个容易被忽略的部分。单节点模型把配电网看成无限大母线但10kV馈线末端电压降和热稳定容量会限制资源实际输出。常见做法是把线性潮流方程加进来对每条支路用输入功率和线路阻抗估算电压偏移ΔV ≈ (P·r Q·x)/V0。这个表达式没有迭代求解精度在规划阶段足够用也可以和短路电流仿真结果一起作为校核条件。# 用简化线性潮流校核三条馈线段末端电压 V0 10.0 # kV母线电压 Z [(0.12, 0.06), (0.15, 0.08), (0.10, 0.05)] # (r, x) 欧姆按段 P_branch [0.8, 0.6, 0.4] # 各段有功传输MW Q_branch [0.2, 0.15, 0.1] # 各段无功传输MVar dv_total 0.0 for (r, x), p, q in zip(Z, P_branch, Q_branch): dv (p * r q * x) / V0 dv_total dv print(末端电压偏移(kV):, dv_total) print(末端电压(pu):, (V0 - dv_total) / V0)这段代码不是优化模型的内部约束而是多类型资源容量方案确定后的一次快速校核。如果末端电压低于0.95说明该段馈线需要加装无功补偿或增大资源就近接入而不是无限放大储能。4. 用Python在10kV配网算例上跑通全过程韧性规划核心代码与参数解释4.1 生成台风场景的负荷与光伏序列先准备一份能代表“全过程韧性”的场景输入。我用24小时作为灾前到灾后恢复的观察窗口双维度数组net表示200个场景、24个时段的净负荷单位是标幺值。历史负荷是均值为1.0、标准差0.15的正态分布光伏在台风期出力下调到平均0.3波动加大这样净负荷序列会留出明显的“缺口窗口”。import numpy as np rng np.random.default_rng(42) N_scen, T_horizon 200, 24 load rng.normal(1.0, 0.15, (N_scen, T_horizon)) load np.clip(load, 0.0, None) pv rng.normal(0.3, 0.25, (N_scen, T_horizon)) pv np.clip(pv, 0.0, 1.0) net np.maximum(load - pv, 0.0)normal取随机数后clip到0以上是为了避免负负荷这种物理上不存在的值。光伏的均值和标准差故意调得比正常天气低且离散代表灾害窗口中的“分布漂移”。这条net数组就是后续分布鲁棒机会约束评估的样本集。4.2 基于重采样的分布鲁棒可靠性评估接下来写一个给定储能容量E_cap和柴发额定功率P_rated时评估200个场景下切负荷概率和平均缺电量的函数。评估逻辑按小时推进储能从60%SOC开始优先用柴发出力供给当净负荷超过柴发功率时放电补差仍然不够的部分记为切负荷。这个逻辑和灾中响应阶段运行方式一致每一步都在回答“如果当时只有这些韧性资源这个小时够不够”。def resilience_power(net, P_rated, E_cap): N, T net.shape soe np.full(N, E_cap * 0.6) # 起始电量60%容量 outage_p 0.0 # 至少出现一次切负荷的场景比例 shortage_acc 0.0 # 平均累计缺电量 eta_dis 0.9 for t in range(T): draw net[:, t] # 柴发先出力仍不足则储能放电 use_gen np.minimum(draw, P_rated) deficit draw - use_gen use_energy np.minimum(deficit, soe * eta_dis) soe np.maximum(soe - use_energy / eta_dis, 0.0) deficit_final deficit - use_energy outage_p (deficit_final 0.01).astype(float) shortage_acc deficit_final.mean() return outage_p / T, shortage_acc / T这里outage_p除以T的含义是“任意小时发生切负荷的平均比例”不是“灾中全过程至少缺电一次的概率”所以它更严厉要求的是每个小时都基本自给这符合全过程韧性里对持续供电能力的强调。柴油机的P_rated是常数模拟的是容量配置下的额定功率暂态爬坡和启停被略去留到验证阶段再补。eta_dis在放电能量换算里出现用来把储能释放的能量对应到SOC消耗和前面模型的效率方向保持一致。4.3 最小成本容量搜索与结果输出有了评估函数接下来在容量参数空间上做网格搜索。这是工程上最笨但最可靠的起步做法不依赖求解器内部算法每一步都能回放。搜索成本由三部分相加——储能单位容量年化成本、柴发单位功率年化成本、缺电补偿单位成本。机会约束在这里体现为硬过滤条件要求每小时平均切负荷比率不超过5%平均缺电量不超过0.05不满足的方案即使总投资更低也不进入候选。E_grid np.arange(0.0, 8.0, 0.5) # 储能容量单位MWh P_grid np.arange(0.0, 3.0, 0.25) # 柴发额定功率单位MW c_e, c_p, c_s 0.05, 0.5, 0.1 # 成本系数单位万元 best_cost np.inf best_E, best_P None, None for E_cap in E_grid: for P_rated in P_grid: v, sh resilience_power(net, P_rated, E_cap) if v 0.05 and sh 0.05: total c_e * E_cap c_p * P_rated c_s * sh if total best_cost: best_cost, best_E, best_P total, E_cap, P_rated print(best E_cap(MWh):, best_E) print(best P_rated(MW):, best_P) print(best total cost:, best_cost)网格步长是精确度和计算量的折中。0.5MWh和0.25MW的网格对10kV馈线规模够用但如果你面对的是一整座园区级微网可以先用大步长扫全域再把最有希望的区域加密换成0.1MWh和0.05MW的步长重跑。搜索结束后不要直接采信这组数字我通常还会再做一次相邻网格点的局部精调。比如best解是E2.5、P1.25就在2.0~3.0和1.0~1.5的范围内用0.1和0.05步长重跑一遍。这一步能碰到网格漏掉的更优解也能判断当前最优解是否在边界上。如果最优解压在探索边界上需要扩大网格范围否则说明可行性区域搜索不完整。这个搜索外层看起来像确定性优化内部却套着分布鲁棒可靠性评估所以结果天然带上了“概率球内最坏姿态”的保守性。要验证这一点下一章会专门讲怎么通过调参踩出这些坑。5. 配置与调参避坑分布鲁棒机会约束规划在10kV配电网落地常见问题5.1 置信度epsilon拍脑袋导致实际风险超标现象规划时机会约束设成95%也就是epsilon0.05自认为已经很保守。第一场强台风过去实际供电缺口天数接近12%可靠性指标明显低于预期。检查代码逻辑没有错约束也确实生效但结果就是不行。原因epsilon是被当成统计风险概率来用的而样本分布是正常天气下的历史分布。灾害窗口的真实分布已经漂移到经验分布半径球外部球内的最坏分布根本覆盖不到换句话说分布鲁棒的机会约束即使严格满足也只是在您自己划定的球内严格球本身划小了。解决用验证集做“最坏违背率”回测。把历史台风案例按时间切在训练集外对每个候选容量跑一遍评估函数记录所有验证场景下违背率的极大值。如果最坏违背率超过epsilon说明球半径r要上调通常的做法是把r从0.5倍初始值逐步拉到2倍观察违背率曲线拐点。我现在一般会把epsilon当作“常规工况余量”另设一个比epsilon更严格的“灾害工况门槛”比如要求5%常规违背率、2%极端违背率两个约束同时满足才允许方案进入下一轮。5.2 储能容量和功率解配比失衡现象优化结果里储能容量接近上限但灾中缺电比例仍然超标。仔细看小时级数据发现很多时刻净负荷缺口不大不小储能放电功率上限被卡死容量再多也放不出来。原因独立设置了储能容量的上下界和充放电功率的上下界但没有把两者通过倍率关系约束起来。优化器可以自由选出一个“大容量、小功率”的储能投资成本和实际缓解能力完全不成比例。解决在储能约束里补上一组倍率限制比如P_dis_max E_cap * 0.5C并把它写成硬约束而不是事后提醒。我还会在目标函数里加入0.2倍的功率成本因为配电网储能的PCS变流器往往比电池本体的单位造价更高不设功率成本解会系统性地偏向大容量方案。仿真时用C-rate0.5和1.0各跑一次观察容量搜索结果有没有明显分化这常常能暴露数据源里的储能参数是否合理。5.3 单节点模型忽略网络耦合约束现象在主变容量和总负荷平衡的规划结果里配置都很合理可一放到10kV短路电流仿真和潮流校核里就露馅馈线末端电压跌破0.9部分线路过载保护定值动作顺序错乱。原因模型把配电网当成了无穷大母线只保证“总供给总需求”没有考虑每条馈线的阻抗、变压器容量和电压降。资源在电气距离上离缺电点太远功率送不过去韧性资源再多也帮不上忙。解决在第一轮规划里至少加入支路容量和末端电压的简化线性约束。把每条馈线的有功输送上限写成P_line 0.8 * S_rated电压约束用( P_r Q_x ) / V0算偏移小于0.05pu。如果规划阶段拿不到无功数据宁可把线路有功限压到0.7也不放开。后面再用综合短路电流仿真做全模型校验把校验发现的越限线路回填到约束里。5.4 场景数太少导致分布鲁棒解抖动现象同一个历史数据集场景数N取20时连续两次运行得到的最优容量相差很大储能容量甚至出现1MWh以上的跳变。原因Wasserstein半径与sqrt(1/N)成正比N小的时候概率球半径大最坏分布对个别极端样本极其敏感N20时随机采样出来的“最坏”场景几乎每换一次种子就换一副面孔优化目标在噪声里跳动。这个现象在分布鲁棒里比普通抽样更明显因为最后还要取一次“球内最坏”。解决把N从20提到至少200优先使用拉丁超立方采样而不是纯蒙特卡洛在概率空间上铺满分位点然后做固定种子的稳定性测试换三组不同随机种子各跑一次如果容量解变化超过20%就继续加样本。更稳的做法是在重复采样时对每个场景都做一次极值扰动让尾部概率密度更厚避免少数场景主导球边界。5.5 量纲不统一导致数值病态现象求解器反复提示“problem might be unbounded”或者某些变量一直停在边界但模型本身没有写错。原因负荷用兆瓦投资成本用万元储能SOC用千瓦时三者数量级混在一起目标函数系数和约束矩阵的尺度差异达到10的6次方以上数值优化器在这种矩阵下会丧失精度。解决第一步统一标幺。在规划模型里选择一个基底容量比如100kW把负荷、光伏、储能出力全都换算成标幺值投资成本换算成“每标幺功率单位”的万元数再和缺电成本加权。这样所有变量保持在0到10的量级目标函数系数在0.001到1之间。做完这一步再回头看求解器的警告绝大多数都消失了。这也是为什么我在第4章的代码里特意先做负荷归一化不是为了好看而是为了少叫一次“求教”。6. 验证分布鲁棒方案的进阶操作回撤测试、最坏违背率表与短路电流校核一个方案能不能上线光看规划阶段的目标函数不够。我会把验证分成三档第一档是回撤式样本外验证第二档是打印最坏违背率表第三档是拿10kV短路电流仿真做最终的设备级校核。三者都过方案才能签字。先说回撤验证。做法是把历史时间序列按年份或按天气过程划分成训练集和验证集训练集用来生成经验分布和球半径验证集只用来做可靠性评估不参与规划。下面的函数记录了训练和验证两个切负荷率并返回其中更大的值def backtest(net_train, net_val, P_rated, E_cap): v_train, _ resilience_power(net_train, P_rated, E_cap) v_val, _ resilience_power(net_val, P_rated, E_cap) return v_train, v_val, max(v_train, v_val) v_tr, v_va, worst backtest(net[:150], net[150:], best_P, best_E) print(train outage:, v_tr) print(validate outage:, v_va) print(worst outage:, worst)如果worst比规划时用的0.05高出不少我不急着降容量而是先去检查验证集里有没有一个特别恶劣的时段把整体数值拉起来。这个时段经常对应光伏出力接近0的连续阴雨台风窗口正好能暴露资源组合的盲区。第二档是最坏违背率表我习惯把以下几个值排成一行打出来规划约束值、训练集违背率、验证集违背率、以及Wasserstein球半径r对应的理论违背率下界。这张表贴在方案书里调度队才能直观理解“我们预留了多少余量”。表格的行不用多五列就够。| 规划约束 | 训练集实际 | 验证集实际 | 最坏违背率 | 球半径r | | 0.05 | 0.032 | 0.058 | 0.071 | 0.20 | | 0.02 | 0.015 | 0.029 | 0.034 | 0.45 |第三档才是短路电流仿真校核。规划模型里的线性潮流只用于筛选方案实际设备能不能承受灾后恢复时的电流冲击要放回基于10kV系统阻抗参数做短路计算。校核点一般是馈线出口断路器、分段开关和用户接入点重点看三相短路电流是否超过设备遮断容量。如果超了就要反过来约束资源接入位置和并网逆变器容量再用新的容量上限重新回第4章的搜索流程。提示短路电流校核一定要用带金属性短路的三相电流值而不是单相接地电流。10kV中性点经消弧线圈接地的系统单相电流可能反而不触发保护容易误导规划决策。我的习惯是每出一版容量方案都把这三档验证结果存在同一个工程目录里半年后复盘时能直接看到当初的“最坏违背率”和实际运行数据的差距。这套习惯帮我躲过了不止一次“看起来很美”的规划方案。希望帮到你。本文还有配套的精品资源点击获取