基于序贯蒙塔卡罗方法进行元件状态抽样附matlab代码

clc; clear all; tic; t0 = clock; rand('state',0) [ NB, NL, NG, databus, dataline, datagen ] = RBTS; databus(:,5) = databus(:,5) .* 1.0; PeakLoad = sum(databus(:,5)); FaultRate = datagen(:,5)'; RepairRate = datagen(:,6)'; Tstate = zeros(1,NG); NumState = 1; DurationComponents = - log( rand(1,NG) ) ./ FaultRate; SystemLoad = PeakLoad; [MinD IndexD] = min(DurationComponents); DurationState (NumState ) = MinD; Bs = 1; NumYear = 0; NumStateSum = 0; EESum = 0; ESum = 0; LOLPSum = 0; EENSSum = 0;IsLossLoad = 0 ;LOLFSum = 0; tic; Iteration_Error = 0.05; while ( Bs>Iteration_Error ) NumYear = NumYear + 1; EnergySheddingSum = 0; DurationSheddingSum = 0; LOLF_Level = 0 ; DurationState = zeros(1,10000); while ( sum(DurationState)<=1 ) GenPower = ~Tstate * datagen(:,3); LoadShedding = SystemLoad - GenPower; [MinD IndexD] = min(DurationComponents); DurationState(NumState) = MinD -sum(DurationState); if ( LoadShedding > 0 ) DurationSheddingSum = DurationSheddingSum + DurationState(NumState); EnergySheddingSum = EnergySheddingSum + LoadShedding * DurationState(NumState) * 8760; if(IsLossLoad == 0 ) LOLF_Level=LOLF_Level+1; IsLossLoad = 1; end else IsLossLoad = 0; end Tstate(IndexD) = ~Tstate(IndexD); if (Tstate(IndexD)==1) DurationComponents(IndexD) = DurationComponents(IndexD) + (-log(rand)/RepairRate(IndexD)); else DurationComponents(IndexD) = DurationComponents(IndexD) + (-log(rand)/FaultRate(IndexD)); end NumState = NumState + 1; end NumStateSum = NumStateSum + (NumState-1); LOLPSum = LOLPSum + DurationSheddingSum; LOLP = LOLPSum / NumYear; EENSSum = EENSSum + EnergySheddingSum; EENS = EENSSum / NumYear; LOLFSum = LOLFSum +LOLF_Level; LOLF = LOLFSum /NumYear ; LOLD = LOLP / LOLF ; EESum = EESum + EnergySheddingSum^2; ESum = ESum + EnergySheddingSum; VAR = EESum/NumYear - (ESum/NumYear)^2; if ( NumYear>10 ) Bs = sqrt( VAR/(NumYear-1) ) / EENS if ( ~mod(NumYear,100) ) NumYear LOLP EENS Bs end end DurationState = zeros(1,10000); DurationComponents = DurationComponents - 1; NumState = 1; end toc; t1 = etime(clock,t0); result_du = [LOLP LOLF EENS Bs t1 NumYear ] fprintf('the Monte Carlo Assessment is over!');