之前的文章简单介绍了配电网可靠性评估解析法和模拟法的基本原理，分别以最小路法和非序贯蒙特卡洛模拟法为例进行介绍。
一、序贯和非序贯蒙特卡罗模拟法的区别
发出来之后很多朋友私信我有没有序贯蒙特卡洛模拟法的程序，以及非序贯蒙特卡洛模拟法和序贯蒙特卡洛模拟法的区别，闲下来就写了这篇博客，跟大家一起交流学习。其实非序贯和序贯蒙特卡洛模拟法原理和步骤比较相似，主要的区别有下面这些：
(1)顾名思义，序贯蒙特卡洛模拟法是基于时序进行模拟的，而非序贯模拟法是基于事件进行模拟的，不考虑时序关系。近年来配电网的研究基本都离不开分布式电源、储能这些时序特性非常明显的设备，因此现在基本上文章都是采用序贯蒙特卡洛模拟法。
(2)在故障场景抽样时，非序贯蒙特卡洛模拟法的每次模拟对应的时间都是一整年，序贯蒙特卡洛模拟法每次模拟的时间不等，是基于设备的TTR和TTF变化的。另外，统计负荷可靠性指标时，两种算法的计算公式也存在明显区别。
具体来说，非序贯每次仿真时，要对每个设备产生一个随机数，如果这个随机数小于设备的故障率，就假定本次模拟中该设备发生故障，也就是在一年内设备会发生故障，这是重点，非序贯蒙特卡罗模拟法每次模拟对应的是时间是一年！！！最后统计负荷可靠性指标时，用累计的停电次数除以总仿真次数就可以得到年平均停电次数（因为一次仿真其实就对应一年），停电时间、失负荷量之类的指标也是同理。
而采用序贯抽样时，每次仿真的时间是不确定的，取决于设备的TTF。对故障场景进行抽样时，故障设备的确定是基于自己设定的规则。一般的文章中都假设抽样时TTF最小的设备为故障元件，这是只考虑了一阶故障。也可以认为TTF最小的前n个设备为故障元件，这样考虑的就是n阶故障。最后统计负荷可靠性指标时，负荷的故障率=仿真中总停电次数/负荷的总正常工作时间，负荷的平均故障持续时间=仿真中总停电时间/仿真中总停电次数。
二、序贯蒙特卡洛模拟法的基本步骤
用序贯蒙特卡洛模拟法进行配电网可靠性评估的步骤如下：
(1)初始化负荷可靠性指标和仿真参数；
(2)计算每个设备的正常工作时间TTF，利用一定规则确定故障的设备；
(3)确定故障设备的修复时间或者上级开关设备的切换时间；
(4)确定因设备故障而停运的负荷，以及因开关设备切换而受影响的负荷，累加负荷的正常工作时间，停运次数、停运时间等。
(5)仿真时间加上TTF和故障持续时间，判断是否达到仿真阈值，若是则进行下一步，否则返回步骤(2)；
(6)统计负荷可靠性指标和系统可靠性指标。
matlab代码：
clc clear close all %% 读取数据 [line_num,line,load_num,load] = IEEE_RBTS_BUS6_F4; %生成图 g=digraph(line(:,1),line(:,2),line(:,3)); %% 利用节点网络分析法确定影响负荷的元件 node_analyse; %% 序贯蒙特卡洛法 load_users=load(:,5); load_fault=zeros(1,load_num);%故障次数 load_fault_time=zeros(1,load_num);%故障时间 load_power=zeros(1,load_num);%失负荷量 Lambda=(line(:,3).*line(:,4))'; Gamma1=line(:,5); Gamma2=line(:,8); power=load(:,3); N=500;%仿真的年数 time=0;%初始化仿真时间 load_normal=zeros(1,load_num);%正常工作时间 while time<N*8760 进行故障抽样 %累加仿真时间 time=time+min_TTF+max([TTR1,TTR2]); end %负荷停电次数 Lambda_load=load_fault./load_normal; %负荷停电时间 Gamma_load=load_fault_time./load_fault; U_load=Lambda_load.*Gamma_load; disp('**************系统可靠性指标**************') SAIFI=Lambda_load*load_users/(sum(load_users)); disp(['SAIFI=',num2str(SAIFI),'次/（年·户）']) SAIDI=U_load*load_users/(sum(load_users)); CAIDI=SAIDI/SAIFI; disp(['SAIDI=',num2str(SAIDI),'小时/（年·户）']) EENS=load_power*load_users/(sum(load_users))/N; disp(['EENS=',num2str(EENS),'MW/（年·户）']) ASAI=1-U_load*load_users/(sum(load_users*8760)); disp(['ASAI=',num2str(ASAI*100),'%']) figure(1) bar(Lambda_load,'stack') title('负荷可靠性评估结果') xlabel('负荷节点');ylabel('年平均停电次数/次'); figure(2) bar(U_load,'stack') title('负荷可靠性评估结果') xlabel('负荷节点');ylabel('年平均停电时间/h'); Z=[SAIFI SAIDI CAIDI ASAI]; figure(3) bar(Z,0.2) title('系统可靠性评估结果') set(gca,'XTickLabel',{'SAIFI','SAIDI','CAIDI','ASAI'}); Lambda_load0=[1.6725 1.6725 1.6725 1.6725 1.6725 1.7115 1.7213 2.5370 2.5890 2.5370 2.5370 2.5370 1.6725 1.7115 1.6725 2.5110 2.5598 2.5110 2.5110 2.5110 2.2250 2.2250 2.2250 ]; U_load0=[8.4015 8.4015 8.4015 8.4015 8.4015 8.5965 8.6453 12.7240 12.9840 12.7240 12.7240 12.7240 11.2875 11.4825 11.2875 15.4800 15.7238 15.4800 15.4800 15.4800 14.0500 14.0500 14.0500 ]; SAIFI0=1.97781318681319; SAIDI0=11.074659340659341; CAIDI0=5.599446608253094; ASAI0=0.998735769481660; figure(4) bar(Lambda_load-Lambda_load0,'stack') title('模拟法误差分析-负荷年平均停电次数') xlabel('负荷节点');ylabel('年平均停电次数的误差'); figure(5) bar(U_load-U_load0,'stack') title('模拟法误差分析-负荷年平均停电时间') xlabel('负荷节点');ylabel('年平均停电时间的误差'); figure(6) Z0=[SAIFI0 SAIDI0 CAIDI0 ASAI0]; bar(Z-Z0,0.4) title('系统可靠性评估误差') set(gca,'XTickLabel',{'SAIFI误差','SAIDI误差','CAIDI误差','ASAI误差'});
三、运行结果
