一、序贯和非序贯蒙特卡罗模拟法的区别

发出来之后很多朋友私信我有没有序贯蒙特卡洛模拟法的程序，以及非序贯蒙特卡洛模拟法和序贯蒙特卡洛模拟法的区别，闲下来就写了这篇博客，跟大家一起交流学习。其实非序贯和序贯蒙特卡洛模拟法原理和步骤比较相似，主要的区别有下面这些：

(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误差'});