基于最优孤岛划分的含DG配电网可靠性评估Matlab实现

1. 项目整体思路为什么可靠性评估要绑上最优孤岛划分1.1 DG接入改变了传统故障后果如果把配电网看成一张单向供应的树状网那么传统可靠性评估其实比较“粗暴”一条线路发生故障保护动作切除故障段之后故障点下游所有负荷都停电直到修复完成。这里几乎不需要考虑“还能不能救回来一部分”的问题因为放射状结构下没有第二个电源。分布式电源DG接入之后情况就变了。光伏、风机、储能这些小电源在故障时如果能与主网解列形成一个独立运行的孤岛就可以继续给孤岛内的负荷供电从而避免整条馈线“全黑”。这也是含DG配电网可靠性评估和传统配电网最本质的区别评估时不仅要问“哪里故障、停多久”还要回答“故障隔离后哪些负荷能由DG带起来、怎么带”。我在做这个项目时最开始也踩过一个认知误区以为只要在网络里加几个DG可靠性指标就会自动改善。实际跑仿真之后发现如果没有做合理的孤岛划分DG甚至可能成为故障扩大的因素。比如孤岛内功率严重不平衡DG过流跳闸最后整个孤岛失电又比如孤岛边界切错了开关非故障区域也跟着被甩掉。所以DG的价值到底能兑现多少完全取决于故障后怎么划分孤岛。1.2 孤岛划分在可靠性评估中的定位可靠性评估的核心是“遍历故障事件统计后果”。传统方法对每个故障事件只需计算影响范围内负荷的停电时间和停电电量。而含DG的项目中每个故障事件都对应一个“恢复策略”要不要孤岛运行哪些开关需要断开哪些负荷保留DG出力怎么分配这就是所谓的“最优孤岛划分”。它不是某一个固定时刻的静态划分而是在每个可能故障场景下以最大化恢复负荷、最小化失电量为目标求解出一个可行的供电孤岛方案。可靠性指标SAIFI、SAIDI、ENS的期望值本质上就是这个“故障事件-孤岛恢复策略-停电后果”循环反复执行后的统计结果。所以整个项目可以拆成三层底层是配电网的拓扑、DG和负荷数据中间层是一系列故障场景顶层是每个场景下的最优孤岛求解和可靠性指标累积。Matlab代码实现时这三层要分离否则后面想换求解器或换DG出力模型会非常痛苦。1.3 适用场景与解决的核心矛盾这个项目的典型使用场景有两类。一类是配电网规划在某个馈线改造方案中规划人员想知道新增一台光伏或储能到底能把用户年平均停电时间降多少值不值得投资。另一类是运行方式校核配电自动化策略已经定了但DG接入后的孤岛方案是否可靠、是否会造成恢复过度或恢复不足需要用定量指标评估。核心矛盾在于精度和计算量之间。可靠性评估要求遍历大量故障场景而“最优”孤岛本身是一个组合优化问题。每个场景都做一次严格最优求解在Matlab里会慢到让人怀疑人生。本文后面会讲到怎么用场景分类、线性化模型、启发式算法来平衡这对矛盾。如果你正准备用Matlab实现类似的项目建议先看1.1和1.2节把模型逻辑理顺再往后看代码拆解。2. 可靠性指标体系与DG出力建模先把输入搞干净2.1 选哪些指标怎么算配电网可靠性评估最常见的指标就是那么几个但每个指标背后的统计口径不能搞错。我做项目时习惯分为三类频率类、持续类、能量类。指标全称计算公式含义单位SAIFI系统平均停电频率年停电用户次数之和 / 总用户数次/用户·年SAIDI系统平均停电持续时间年停电用户时数之和 / 总用户数小时/用户·年CAIDI用户平均停电持续时间SAIDI / SAIFI小时/停电用户次ASAI供电可用率1 - SAIDI / 8760无量纲ENS年总缺供电量各故障事件缺供电量之和kWh/年AENS平均缺供电量ENS / 总用户数kWh/用户·年传统解析法的计算方式是这样的对每个线路元件假设其年故障率为λ修复时间为r则下游用户每年承受的停电时数为 λ × r × 下游用户数。把所有元件累加SAIFI的分子累加 λ × 受影响用户数SAIDI的分子累加 λ × r × 受影响用户数。但含DG后停电持续时间不再是简单的修复时间r。如果故障后某个负荷能通过孤岛快速恢复那么它只经历开关倒闸操作时间t_sw比如0.5小时而不是一直等到故障排除。所以实际计算时要区分两类负荷能转入孤岛的负荷停电时间为t_sw不能转入孤岛的负荷停电时间为r。这一点很多入门代码会忽略结果会让SAIDI明显偏乐观。2.2 光伏、风机、储能的出力时序模型可靠性评估不是只看额定容量因为故障发生时DG到底能发多少电取决于当时的光照、风速和储能荷电状态。国内工程评估中常见做法是采用时序负荷和时序出力按8760小时逐小时模拟或者将一年分成若干个典型时段。光伏出力我建议用“理论辐照度曲线 × 随机折减系数”。理想晴天曲线可以按当地纬度和日辐照峰值拟合再用Beta分布或正态分布随机扰动模拟多云天气。风机出力则先根据威布尔分布生成小时风速序列再用典型风机功率曲线折算[ P_{WT} \begin{cases} 0 v v_{in}\ P_r \frac{v^2 - v_{in}^2}{v_r^2 - v_{in}^2} v_{in} \le v v_r\ P_r v_r \le v v_{out}\ 0 v \ge v_{out} \end{cases} ]储能模型比DG复杂一些因为涉及荷电状态SOC的时间连续性。在可靠性评估里如果把储能当成一个“随时满电”的电源结论会严重失真。我自己的做法是在孤岛仿真开始前根据前一时刻的运行工况初始化SOC再限制孤岛持续时间内储能放电量不能超过可用容量上限。Matlab代码里对应的是一个简单的状态转移函数而不是固定出力点。2.3 负荷模型与故障参数负荷建模同样重要。可靠性统计要求按“用户数”加权而潮流约束要求按“功率”平衡两者不能混为一谈。通常每个10kV节点会挂一个配变配变下假设若干用户节点有功和无功负荷可以用峰值负荷并乘以时序曲线得到。另外要提前准备线路故障参数单位长度故障率、平均修复时间、开关操作时间。不同线路类型故障率差别很大电缆和架空线不要混用同一组参数。如果数据来自某地区统计注意单位是“次/百公里·年”还是“次/公里·年”换算错一个数量级可靠性指标会差得离谱。我习惯把所有输入参数集中放在一个结构体baseData里统一管理包括拓扑矩阵、支路阻抗、负荷表、DG表、用户表、故障参数表。这样后面写故障枚举和孤岛求解时输入输出都很清晰。3. 最优孤岛划分的实现策略数学建模与算法选型3.1 孤岛划分问题的约束体系最优孤岛划分本质上是一个带约束的0-1整数规划问题。决策变量至少包括每个节点是否被划入孤岛、每条支路是否投入运行、每个DG是否孤岛运行。目标函数一般是最小化失负荷量或者最大化恢复负荷的加权和权重可以设置为负荷等级或用户数。约束条件必须覆盖四个方面第一是连通性约束。孤岛内部不能出现“电气孤岛内还有孤立节点”这种病态结果。最简单的做法是引入网络流变量让孤岛内每个被选中节点都能通过投入的支路连接到一个公共的“根节点”这个根节点通常是DG接入点。第二是功率平衡约束。孤岛内DG出力加储能放电必须满足孤岛内所有保留负荷的有功和无功需求。这个约束可以写成不等式保住负荷的总功率不超过可用电源总功率。第三是线路容量与电压约束。风机、光伏接入后孤岛内的潮流可能重新分配某些线路可能过载。精确做法是求解交流潮流但可靠性评估场景多、计算量大我一般用线性化DistFlow模型忽略无功的细节损耗保留线路传输功率上限和节点电压幅值范围。第四是DG与储能运行约束。DG出力不能超过当前时刻最大可用出力储能充放电功率受逆变器容量限制孤岛持续时间要满足SOC上下限约束。3.2 混合整数线性规划与启发式算法怎么选这个项目标题里“最优”两个字意味着理论上要做“精确优化”。Matlab里最顺手的路径是用YALMIP搭建MILP模型接入商业求解器或者开源的CBC求解器。对于IEEE 33节点这种规模单次求解时间在几百毫秒到几秒不等可靠性评估中如果只枚举“N-1故障”总共也就32个线路故障场景左右是可以接受的。但如果你的网络规模上百个节点或者要做序贯蒙特卡洛模拟成千上万次每次都调用MILP就会很吃力。这时候我建议分档处理对于规划评估、年度可靠性指标计算使用严格MILP每个故障场景求解上限设20秒求解失败就用启发式兜底。对于在线运行、日内滚动评估可以使用二进制粒子群BPSO或遗传算法把目标函数稍作松弛快速给出一个可行孤岛方案。还有一种更工程化的做法先基于拓扑搜索生成“可孤岛运行区域候选集”再在每个候选集内做功率分配。这种方法类似“预生成孤岛方案表”可靠性评估时查表而不是实时优化速度非常快。我在实际项目中用的是“MILP为主、BPSO为辅”的组合。因为可靠性评估是离线计算求精确解更有说服力也能作为启发式算法的参照基准。3.3 一个算例说明“最优”的价值简单举一个例子。某馈线有12个节点DG接在节点6容量为1.5MW。故障发生在节点1到节点2之间的主干线路上。主网断开后节点2到节点12都失去主网供电。如果采用“就近恢复”的简单策略可能只把节点6、7划入孤岛因为这两个节点离DG最近却能恢复的负荷只有0.8MW。如果采用最优孤岛划分可以闭合节点8旁边的分段开关将节点8、9、10的部分负荷也纳入孤岛前提是DG容量还有富余。假设这部分负荷是1.0MW那么最优方案比就近方案多恢复0.2MWES缺供电量就少800kWh按4小时修复时间算。所以“最优”不是学术上的自嗨。它直接决定了可靠性评估结果中DG对停电损失削减的贡献值。只做静态就近孤岛的代码写起来简单但结果会保守到让人觉得DG“没什么用”。4. Matlab代码实现从数据组织到核心模块拆解4.1 文件结构与输入数据准备这个项目我建议按下面这样的文件结构组织方便后期维护和扩展文件/函数作用setup_data.m定义拓扑、负荷、DG、故障参数run_reliability.m主程序故障枚举与指标计算入口fault_island.m对单个故障场景求解最优孤岛dgp_forecast.m生成光伏/风机/储能的时序出力calc_indicators.m累加SAIFI、SAIDI、ENS等指标plot_reliability.m绘制指标对比图、孤岛拓扑图输入数据准备阶段最容易出错的是支路表。配电网支路表要按方向严格定义首端和末端否则后面做连通性检查时会乱。我通常这样组织每行记录支路编号、首端节点、末端节点、电阻、电抗、长度、容量上限、单位故障率、平均修复时间。负荷表必须同时包含有功和无功否则后面做潮流校验时会缺数据。用户数这一列不要忘记SAIFI和SAIDI的计算全靠它作为权重。DG表要写清楚类型、接入节点、额定容量、时序出力文件名。如果以后要加入电动汽车充电桩直接在DG表里新增一行可靠性主程序不用改。4.2 主循环框架故障枚举与可靠性累计可靠性评估主程序的核心逻辑其实很短。下面是我常用的骨架% run_reliability.m 主程序伪代码 baseData setup_data(); nBranch size(baseData.branch, 1); sum_freq 0; % 停电用户次数 sum_duration 0; % 停电用户时数 sum_ens 0; % 缺供电量 for k 1:nBranch % 设置第k条支路为故障状态 faultLine baseData.branch(k, 1); % 求解该故障场景下的最优孤岛 result fault_island(baseData, faultLine); % 结果累加 sum_freq sum_freq baseData.lambda(k) * result.users_affected; sum_duration sum_duration baseData.lambda(k) * result.hours_affected; sum_ens sum_ens baseData.lambda(k) * result.ens; end totalUsers sum(baseData.load.userCount); SAIFI sum_freq / totalUsers; SAIDI sum_duration / totalUsers; ASAI 1 - SAIDI / 8760; ENS sum_ens;这里有一个容易忽略的细节result.hours_affected并不是所有受影响用户都等于故障修复时间。在fault_island函数里要返回两部分用户数据切入孤岛的用户停电时间是切换时间t_sw未恢复的用户停电时间是修复时间r。主循环累加时按用户数分别累加最后得到SAIDI才是正确的。如果想要更高的精度可以把故障枚举改成序贯蒙特卡洛抽样。做法是用均匀分布随机抽样生成每年的故障时刻和故障元件然后重复多年模拟。此时主循环不再是“遍历所有支路”而是“遍历所有随机事件”。代码结构不变只是把lambda乘法逻辑改为事件计数。4.3 孤岛划分求解器的代码实现fault_island是整个项目最核心的函数。这里给出一个基于YALMIP的简化MILP实现思路function result fault_island(baseData, faultLine) n size(baseData.bus, 1); loadP baseData.load.active; dgP baseData.dg.maxActive; % 当前时刻最大可用出力 branch baseData.branch; % 决策变量 x binvar(n, 1); % 节点是否被恢复 y binvar(size(branch,1), 1); % 支路是否投入运行 % 目标最大化恢复负荷 objective - sum(loadP .* x); % 约束 constraints []; % 根节点DG所在节点必须恢复 constraints [constraints, x(findDGNode) 1]; % 支路投入后两端节点必须恢复 constraints [constraints, y x(branch(:,2))]; constraints [constraints, y x(branch(:,3))]; % 功率平衡DG总出力不小于恢复负荷总和 constraints [constraints, sum(dgP .* x) sum(loadP .* x)]; % 线路容量约束 constraints [constraints, y .* ... capacity]; % 线性化处理 % 连通性约束需要额外引入网络流变量这里省略细节 ops sdpsettings(solver, milp, verbose, 0, savesolveroutput, 1); sol optimize(constraints, objective, ops); % 输出恢复节点集合与失负荷量 result.nodes find(value(x) 0.5); result.ens (1 - sum(loadP .* value(x))) * outageHours; end注意上面的代码为了展示核心逻辑非常简化实际应用时连通性约束必须补上。一个常用的线性化办法是引入从DG点出发的单商品流变量对每条投入支路和节点建立流量平衡约束。这样才能保证解出来的节点集合能够构成一个连通的孤岛。如果不方便接MILP求解器也可以用BFS做候选孤岛生成再在候选集上用线性规划分配功率。BFS的优势是不需要额外安装求解器Matlab自带的图算法就能跑缺点是候选集选择不好时容易陷入局部最优。所以我在代码里保留了两种求解接口用options.useMILP开关切换方便对比。4.4 指标汇总与负荷点可靠性输出故障枚举完成之后除了输出系统级指标我还建议输出每个负荷节点的可靠性明细。因为规划时经常需要看“某个偏远负荷点为什么SAIDI高”是线路故障率太高还是孤岛恢复策略没覆盖到它。节点级指标累加逻辑类似对每个故障事件如果该节点在孤岛恢复集合内则它只停电t_sw小时如果不在集合内则停电r小时。把各个节点的年停电时数除以8760就能得到节点供电可用率。Matlab里用矩阵累加可以写得非常简洁node_outage_time zeros(n, 1); for k 1:nBranch node_outage_time node_outage_time lambda(k) * outageDurationVector{k}; end其中outageDurationVector{k}是长度为n的向量记录了第k个故障下每个节点的实际停电小时数。这个向量来自fault_island的返回结果。再用bar函数画节点SAIDI对比图效果很直观。5. 实操中的参数整定、求解器问题与避坑记录5.1 求解时间与收敛性控制我最初跑这个项目时用的是YALMIP加MILP求解器对IEEE 33节点做N-1故障枚举每个故障场景都严格求解整体耗时约两分钟。看起来还能接受但一旦把故障枚举扩大到含支路和变压器所有元件或者加入多故障场景计算量会指数上升。控制求解时间有几个实用招数。第一在sdpsettings里设置cputime上限比如单场景20秒超过就接受当前最优解。第二先做拓扑简化把故障段下游区域之外、且不可能被DG恢复的节点直接排除缩小优化变量规模。第三把DG的时序出力按“低、中、高”三个典型水平聚类每个故障场景只在对应的出力水平下求解一次然后用概率加权而不是每个时序点都求解一次。收敛性问题更多出现在求解器版本和YALMIP兼容性上。如果遇到“infeasible problem”但理论上明明有可行解先检查是不是根节点约束没加对或者故障支路在决策变量y中没有被强制置0。我踩过好几次这样的坑故障支路虽然断开了但优化模型以为它还能用结果孤岛把故障线路又“接回去”了。解决办法是在约束里强制y(faultLine) 0。5.2 常见问题速查表问题现象可能原因解决建议孤岛结果不连通缺少连通性约束加入单商品流约束或割集约束故障支路仍出现在恢复方案中未强制故障支路为0在模型中显式设置y(faultLine)0DG出力明明不足却恢复了更多负荷功率平衡约束写成了等式且允许负功率改为不等式并限制DG输出上限SAIDI明显偏乐观未区分孤岛切换时间与修复时间用t_sw和r分别累加停电用户时数储能型DG可靠性贡献过大储能SOC始终按100%计加入SOC时序约束限制孤岛持续供电能力不同求解器结果差异大MILP求解器对数值容差不同统一设置求解器精度并对比最优目标值部分负荷点用户数为0用户数缺失导致指标异常检查负荷表用户数至少为1这个表是我项目整理过程中真实遇到的问题不是空写尤其是“故障支路仍出现在恢复方案中”这个问题代码里如果没有强制置零结果会非常荒谬甚至出现孤岛通过故障线路向主网反送电的假象。5.3 结果合理性怎么验证做完仿真不能只打印一个SAIDI就完事。我通常用三个步骤验证结果是否可信。第一步无DG对比测试。把DG容量全部置为0跑一遍传统可靠性评估。此时SAIFI、SAIDI和ENS应该与手算结果一致因为孤岛划分退化为“无孤岛恢复”所有受影响负荷都停到修复完成。如果这一步都对不上说明故障枚举或指标统计逻辑有问题。第二步极限测试。将DG容量设得非常大比如10倍负荷理论上所有故障后的下游负荷都能被孤岛恢复那么SAIDI应当约等于切换时间t_sw而不是修复时间。如果指标明显高于t_sw说明孤岛划分的约束条件限制了DG能力需要检查是不是线路容量或连通性约束过严。第三步求解器一致性测试。对同一个故障场景分别用MILP和启发式算法求解比较恢复负荷量。MILP的目标值应该不差于启发式。如果启发式反而更好先检查是不是MILP模型里漏掉了某些恢复路径再检查求解器是否因为超时返回了次优解。6. 个人体会与后续扩展建议我在复盘这个项目时最大的体会是最优孤岛划分虽然在可靠性评估里非常有用但真实工程中还要考虑保护配合、自动化开关动作失败率、通信时延等问题。模型里的“最优”是一个理想上限实际能兑现多少取决于现场自动化设备能不能准确执行这个孤岛方案。所以我在报告里一般会同时给两种结果理想最优孤岛下的可靠性指标以及考虑开关拒动概率后的可靠性指标后者更能贴近运行实际。代码扩展方面我建议下一步可以加入分布式储能的多时间尺度协调策略因为孤岛持续期间的储能SOC变化非常关键。另一个可以扩展的方向是负荷分级把重要负荷、一般负荷分为不同权重最优孤岛划分的目标函数从“最大化恢复电量”改为“最大化重要负荷恢复比例”。这会直接改变孤岛划分结果也更符合实际保电需求。最后说一个小技巧Matlab代码里我习惯在每一轮故障枚举后都保存一份result日志包括恢复节点集合、未恢复负荷节点、DG出力水平、求解耗时。这样即使后面指标对不上也能定位到具体是哪个故障场景出了问题。这个习惯帮我在项目后期省了大量排查时间建议你也试试。