简介面向配电网规划与可靠性研究者的Matlab程序包围绕含分布式电源DG的配电网可靠性评估提供两种建模思路一是基于概率模型的最小路法二是基于时序模型的序贯蒙特卡洛模拟法参考《基于仿射最小路法的含分布式电源配电网可靠性分析》方法编写内置IEEE RBTS-BUS6 F4测试系统可直接运行得到可靠性指标。压缩包共8个文件其中5个.m脚本按主程序、节点分析、DG概率模型、蒙特卡洛时序模型等模块划分1个.mat保存光伏数据另有节点编号示意图vsdx和代码说明文档docx整体仅111KB轻量紧凑便于逐模块对照学习。代码注释清晰适合高校研究生或电力工程师快速理解最小路法与序贯蒙特卡洛模拟的编程实现并在此基础上扩展分布式电源渗透率、线路故障率等参数分析。目前已有396人学习下载对正在开展分布式电源接入可靠性研究的学生和从业者是一份可直接复用的实用参考资料。1. 分布式电源接入后传统可靠性评估为何不再可信去年做配电网规划项目时碰到一个现象某条10kV馈线接入了分布式光伏按照传统可靠性评估方法算出来的用户年平均停电时间和实际统计值差了将近三成。原因在于传统方法默认系统是单电源辐射状运行故障发生后负荷只能等主网侧修复而分布式电源改变了故障传播路径——故障隔离区间内的部分负荷原本只能停电等待现在有机会由DG继续供电。所以《基于仿射最小路法的含分布式电源配电网可靠性分析》这篇文献里作者用概率模型和时序模型分别去刻画DG出力不确定性对应了两套评估方案。下面用IEEE RBTS BUS6 F4馈线作为测试系统拆解这套MATLAB代码中最小路法和序贯蒙特卡洛模拟法的完整实现与参数设计适合正在做含DG配电网规划、可靠性指标计算的同学参考。2. 最小路法与RBTS BUS6 F4先把网络拓扑拆清楚2.1 最小路法的核心思想从负荷点反向看网络最小路法的基础假设很朴素对辐射状配电网中的某个负荷点从电源节点到该负荷点存在一条唯一的供电路径这条路径就是该负荷点的最小路。最小路上的任意元件故障都会直接导致负荷点停电路径之外的元件故障只有在其影响能通过分段开关或联络开关传导到这条路径时才会间接影响负荷点。因此计算被拆成两步最小路上的故障率直接叠加最小路外的故障率通过影响因子折算叠加到负荷点等值故障率上。DG接入后这个逻辑多了一个分支当主电源因上游故障断开时负荷点如果处在DG的孤岛运行范围内且DG当前出力能覆盖孤岛内负荷那这段时间实际上不停电。于是一个负荷点的可靠性不再只由网架结构决定还取决于DG出力概率分布和孤岛成功概率。这也是这套代码里概率模型和最小路法要放在一起的原因。2.2 RBTS BUS6 F4系统参数文件怎么组织IEEE_RBTS_BUS6_F4.m 本质上是数据文件不是计算脚本。它把F4馈线的节点、支路、负荷、用户数全部写成MATLAB可读矩阵后续所有函数直接引用这些变量不需要每次重新输入网络参数。F4馈线是RBTS BUS6测试系统的一条主馈线包含7个负荷点节点总数在20个左右规模不大但分段开关、联络开关等配置完整用来验证可靠性评估算法很合适。%% IEEE_RBTS_BUS6_F4.m —— F4馈线系统参数定义格式示例 % 节点矩阵节点编号 | 有功负荷(kW) | 无功负荷(kVar) | 用户数 NODE [ 1 0.000 0.000 0; % 电源节点BUS6母线 2 0.000 0.000 0; % 分段节点 3 96.870 60.880 134; % 负荷点LP1 4 0.000 0.000 0; % 分支节点 5 181.630 114.130 134; % 负荷点LP2 6 174.040 109.370 134; % 负荷点LP3 7 193.740 121.760 134; % 负荷点LP4 8 0.000 0.000 0; % 联络节点 9 74.060 46.540 134; % 负荷点LP5 10 173.370 108.950 134; % 负荷点LP6 11 95.410 59.960 134; % 负荷点LP7 ]; % 支路矩阵支路编号 | 首端节点 | 末端节点 | 故障率(次/年) | 修复时间(h/次) | 长度(km) BRANCH [ 1 1 2 0.1036 5.00 1.30; 2 2 3 0.1086 5.10 1.36; 3 2 4 0.0982 5.00 1.00; 4 4 5 0.1136 5.00 1.42; 5 2 6 0.1166 5.10 1.38; 6 6 7 0.1098 5.00 1.35; 7 4 8 0.1042 5.00 1.24; 8 8 9 0.1074 5.00 1.40; 9 8 10 0.1122 5.10 1.32; 10 10 11 0.1018 5.00 1.18; ];第一个代码块定义了节点数据NODE矩阵第1列是节点编号必须与BRANCH矩阵的首末端节点对应第2、3列是有功和无功负荷用于ENS指标和孤岛校验第4列用户数直接参与SAIFI、SAIDI的加权计算不能省略。BRANCH矩阵第4、5列是每个元件的故障率和平均修复时间最小路法对这两个参数非常敏感修改变量时优先检查这里。具体数值以脚本内RBTS标准数据为准上面只是格式示例。矩阵字段含义单位影响指标NODE第1列节点编号与BRANCH首末端对应无拓扑搜索NODE第2~3列有功/无功负荷kW/kVarENS、孤岛校验NODE第4列该节点用户数户SAIFI、SAIDIBRANCH第4列元件故障率lambda次/年全部可靠性指标BRANCH第5列平均修复时间r小时/次SAIDI、ENSBRANCH第6列馈线段长度km按长度折算故障率2.3 node_analyse.m节点分类与最小路搜索的前处理node_analyse.m 负责把BRANCH矩阵转换成图结构并对每个负荷点执行最小路搜索。F4馈线是开环运行的节点之间的父子关系可以通过一次深度优先搜索确定搜索过程中记录每个节点的父节点找到目标负荷点后逐级回溯还原整条路径。function path search_min_path(adj_matrix, source, target) % 基于DFS搜索最小路返回从source到target的节点序列 n size(adj_matrix, 1); visited false(1, n); stack source; parent zeros(1, n); visited(source) true; while ~isempty(stack) cur stack(end); if cur target, break; end nxt find(adj_matrix(cur, :) ~visited); if isempty(nxt) stack(end) []; else n nxt(1); visited(n) true; parent(n) cur; stack(end1) n; end end path target; while parent(path(end)) ~ 0 path(end1) parent(path(end)); end path fliplr(path); end这段代码用栈模拟DFS迭代过程parent数组记录每个节点首次被访问时的前驱节点找到target后通过回溯还原最小路。nxt(1)取第一个未访问邻居在辐射状网络中不会产生分支歧义。实际使用时先由BRANCH矩阵构造邻接矩阵再传入本函数。node_analyse.m还会顺带识别联络节点和分段节点——这两类节点决定了非最小路支路故障时的影响因子后面算等值故障率时需要用到。3. 概率模型 最小路法DG出力如何进入可靠性计算3.1 为什么DG出力要建模成概率分布分布式光伏出力不是一个固定值它受光照、云层、温度影响在一天内呈现强烈的时序变化。可靠性计算如果直接用额定功率或全天平均值会带来系统性偏差中午光伏出力高孤岛运行成功率被高估傍晚和夜间出力低孤岛供电能力被低估。正确做法是把光伏出力看作随机变量用概率分布描述其取值可能性。pv_data.mat 里存放的就是光伏出力的样本数据通常为标幺化时间序列概率模型从这些样本中统计出力在各区间的概率供最小路法使用。3.2 DG_probability_model.m出力状态枚举与离散化DG_probability_model.m 的核心工作是把连续出力曲线离散化为M个状态并计算每个状态的概率。常见做法是把出力按额定容量比例等距划分区间统计样本落在每个区间的频率。function state DG_probability_model(pv_pu, M) % pv_pu: 光伏出力标幺值序列来自pv_data.mat % M: 出力状态数 % state: 结构体含各状态出力值及概率 pi linspace(0, 1, M1); % 划分区间边界 state.p zeros(1, M); % 各状态概率 state.value zeros(1, M); % 各状态出力值取区间中点 for k 1:M state.value(k) (pi(k)pi(k1))/2; idx pv_pu pi(k) pv_pu pi(k1); state.p(k) sum(idx) / length(pv_pu); end state.p(end) state.p(end) sum(pv_pu 1) / length(pv_pu); end逻辑说明先按等间距划分M个区间统计每个区间的频率作为状态概率状态出力值取区间中点。最后一行单独处理出力等于1的边界样本否则有一部分峰值出力会被漏掉。M的取值直接决定精度与计算量的平衡工程上推荐按下面的表选择。M取值出力状态粒度适用场景5粗粒度快速估算、趋势对比10~20中等粒度RBTS级配电网可靠性评估30以上细粒度高渗透率、时序相关性显著3.3 仿射最小路法的关键公式与main_minimal_path.m有了DG出力状态概率下一步把孤岛成功概率引入可靠性指标计算。仿射最小路法和普通最小路法的区别在于前者用区间或仿射形式表示DG出力的不确定性而不是简单取期望值算出来的指标往往是一个范围能直观反映DG出力波动对可靠性指标的影响边界。其中最关键的两个量是孤岛成功概率P_iso和DG可用率avail_DG。% main_minimal_path.m 主流程节选 % 载入系统参数与负荷数据 IEEE_RBTS_BUS6_F4; load(pv_data.mat, pv_pu); lambda_DG 5.0; % DG等效故障率次/年 r_DG 20; % DG平均修复时间小时 rated_DG 800; % DG额定容量kW M 15; DG_state DG_probability_model(pv_pu, M); for i 1:length(LP_node) % 遍历所有负荷点 path_i search_min_path(adj, source, LP_node(i)); branch_in find_branch_in_path(BRANCH, path_i); branch_out setdiff(1:size(BRANCH,1), branch_in); lambda_LP sum(BRANCH(branch_in, 4)); % 最小路故障率累加 for j 1:length(branch_out) % 非最小路支路按影响因子折算与隔离开关位置有关 lambda_LP lambda_LP BRANCH(branch_out(j), 4) * ... isolation_factor(branch_out(j), LP_node(i)); end % 孤岛成功概率DG出力状态覆盖孤岛负荷的累积概率 P_iso sum(DG_state.p(DG_state.value * rated_DG ... isolated_load(i))); % 负荷点年平均停电时间小时/年 U_LP(i) (lambda_LP lambda_DG) * 8760 * ... (1 - P_iso * avail_DG) / 8760; end这段代码体现了两个关键点。第一lambda_LP由最小路支路直接累加和非最小路支路折算叠加得到isolation_factor由隔离开关分段位置决定取0到1之间的值0表示该支路故障不会影响该负荷点1表示全影响。第二孤岛成功概率P_iso是把DG出力离散状态的累积概率求和条件是DG容量乘以当前出力标幺值能覆盖孤岛内总负荷isolated_load(i)。注意U_LP用的是年停电小时数lambda_LP单位是次/年两者相乘直接得到小时/年不需要额外换算。4. 时序模型 序贯蒙特卡洛模拟把8760小时逐时跑一遍4.1 概率最小路法算不出哪天停电概率模型最小路法算的是长期期望值处理不了时序相关问题。光伏出力有很强的日内规律和季节规律——夏天午间出力高峰、冬季晚峰出力为零概率模型把这些信息全部抹平了。工程上常说概率模型适合方案比选时序模型适合精确校核。若光伏渗透率超过20%概率模型和时序模型算出的SAIDI差距会明显拉大。这时候要用序贯蒙特卡洛模拟法对每个元件建立正常运行和故障停运两状态模型通过抽样生成整个仿真周期的状态持续时间序列在系统层面按小时推进统计每个负荷点的停电次数和停电时长。4.2 状态持续时间抽样两状态模型的数学基础两状态模型假设元件的正常运行时间和故障修复时间都服从指数分布。指数分布的无记忆性使得状态持续时间可以通过逆变换法抽样U是[0,1)均匀随机数则正常运行时间TTF -1/lambda * ln(U)修复时间TTR -1/mu * ln(U)其中mu 1/MTTR。function [TTF, TTR] sample_state_transition(lambda, mu, hours) % lambda: 故障率次/年mu: 修复率次/小时 % hours: 需要生成的仿真时长小时 % TTF/TTR: 正常持续时间、修复持续时间序列小时 n ceil(hours / 200) 10; % 预留足够的事件段数 TTF zeros(1, n); TTR zeros(1, n); for k 1:n TTF(k) -log(rand()) / (lambda / 8760); % 年转小时 TTR(k) -log(rand()) / mu; end end这里特别要注意lambda的单位换算lambda是次/年抽样时间单位是小时因此lambda要除以8760转换为次/小时。修复率mu直接使用次/小时通常mu 1/MTTR_hour。如果在单位换算上出错最常见的现象是模拟出的停电次数比实际大8760倍这个坑几乎每个初学者都会踩一次。4.3 系统状态评估与指标统计生成所有元件的状态序列后需要按小时推进仿真时钟每个小时检查哪些元件处于故障状态判断故障隔离范围再根据DG出力和孤岛内负荷判断孤岛是否成立最后更新每个负荷点的供电状态。for yr 1:sim_years % 仿真年数 for h 1:8760 for k 1:n_comp % state1表示故障remain为剩余状态持续时间 remain(k) remain(k) - 1; if remain(k) 0 if state(k) 0 state(k) 1; remain(k) TTR_next(k); % 进入故障状态 else state(k) 0; remain(k) TTF_next(k); % 恢复正常 end end end % 根据当前元件状态集合确定停电负荷点 % 判断孤岛内DG出力 孤岛负荷决定是否转孤岛供电 % 累加负荷点停电次数与停电小时数 end % 统计第yr年的SAIFI、SAIDI、ENS end注意性能问题n_comp为20、sim_years为5000时内层循环共执行20500087608.76亿次纯MATLAB逐小时仿真非常慢。常见优化是改成事件驱动仿真每次跳跃到下一个状态变化时刻而不是逐小时扫描所有元件同样精度下耗时能下降一个数量级。主程序里如果默认仿真年数较大建议优先检查是否采用了事件驱动写法。仿真参数推荐取值说明sim_years2000~5000太少不收敛太多耗时剧增n_comp按系统元件数包含线路、开关、DGlambda单位次/年抽样时除以8760mu单位次/小时即1/MTTR收敛判据方差系数1%主要看SAIFI的收敛性5. 可靠性指标对比的一个实操技巧指标差异定位模型误差5.1 用灵敏度扫描找出两套模型的差异来源概率最小路法和序贯蒙特卡洛法算出的SAIFI如果对不上先不要怀疑程序写错了而要想是哪层假设出了问题。最小路法在折算非最小路元件故障时隐含了故障事件相互独立的假设而时序模拟中同一时段多个元件同时故障虽然概率小但高渗透率DG接入后DG停运和线路故障的联合影响会被放大。一个常用的验证手段是扫描DG渗透率观察两种方法的指标差值曲线。把DG额定容量从总负荷的0%逐步加到50%分别运行两套评估程序绘制SAIFI随渗透率的变化曲线。penetration 0:0.1:0.5; for k 1:length(penetration) rated_DG penetration(k) * total_load; SAIFI_minpath(k) run_minpath(NODE, BRANCH, rated_DG); SAIFI_mc(k) run_mc(NODE, BRANCH, rated_DG, 2000); end plot(penetration, SAIFI_minpath, -o, penetration, SAIFI_mc, --s); xlabel(DG渗透率); ylabel(SAIFI次/用户·年); legend(概率最小路法, 序贯蒙特卡洛法);如果差值是稳定平移说明问题出在等值折算的精度上可以把isolation_factor细化到每个开关段如果差值随渗透率明显张开说明时序相关性起了主要作用这时应当以序贯蒙特卡洛结果为准最小路法只用于趋势分析。绘制曲线时可以顺带把SAIDI和ENS也画出来三者差值形态能进一步区分到底是故障率折算误差还是孤岛概率建模误差。另一个实用技巧是检查故障率单位。RBTS标准系统的lambda论文里有的给次/年有的给次/小时混用时指标会差8760倍。启动程序前先跑一段10年的小规模模拟把SAIFI手算一遍与程序输出对比验完单位再跑正式工况。本文还有配套的精品资源点击获取