Matlab code of Subset Simulation with Monte-Carlo Simulation for reliability analysis

Matlab code of Subset Simulation with Monte-Carlo Simulation for reliability analysis

function [pfss, gRe, cdf, zRe] = SS(Fun, n, paras)

% Algorithm parameters
N = paras.NumSam;
p0 = paras.CondPro;
nt = round(N*p0);
ns = ceil(1/p0-1);
%% Step 1 (Monte Carlo simulation)
z = normrnd(0,1,N,n); % Generate standard normal samples
% Calculate LSF
for i=1:N; 
    g(i) = feval(Fun,z(i,:));
end;
% Sort samples and LSF values
[gSort,index]=sort(g);
Z = z(index,:);
% Record the calculated results
pfss = nt/N;
gRe = gSort;
cdf = (N-[N:-1:1]+1)/N;
zRe = Z;
%% Step 2 (Markov Chain Monte Carlo)
% preparation for MCMC
sigma = std(Z(1:nt,:),0,1);
stopFlag = 0;
iter = 1;
while (stopFlag == 0)
    w = Z(1:nt,:); 
    g = gSort(1:nt);
    iter = iter+1;
    len = nt+1;
    for i=1:nt
        seed = w(i,:); 
        seed_g = g(i);
        for j = 1:ns % A Markov chain
            for k = 1:n % Component by component
                % Generate an candidate            
                u(k) = seed(k)+(2*rand()-1)*sigma(k);
                % Calculate the acceptance probability
                pdf2 = exp(-0.5*seed(k).^2);
                pdf1 = exp(-0.5*u(k).^2);
                alpha = pdf1/pdf2;
                % Accept u(k) with a probability of alpha
                if  alpha > rand(); 
                    w(len,k)= u(k);
                else 
                    w(len,k)= seed(k);
                end;
            end 
            gTemp = feval(Fun,w(len,:));
            % Accept or reject 
            if gTemp <= gSort(nt+1);
                g(len) = gTemp;
                seed = w(len,:); 
                seed_g = g(len);
            else
                g(len) = seed_g;
                w(len,:) = seed;              
            end
            len = len+1;
            % terminate the simulation of a Markov chain
            if len > N break; end
        end
        if len > N break; end
    end
    % Sort samples and LSF values
    [gSort,index]=sort(g);
    Z = w(index,:); 
    % Record the calculated results
    gRe = [gSort,gRe(nt+1:end)];
    cdf = [(N-[N:-1:1]+1)/N*(nt/N)^(iter-1),cdf(nt+1:end)];
    zRe = [Z;zRe(nt+1:end,:)];
    % Terminate simulation   
    if iter >= paras.MaxTry || gSort(nt+1) <= 0
        pfss = pfss*(size(find(g<0),2)/N);
        stopFlag = 1;    
        break;
    else
        pfss = pfss*(nt/N);
    end;
end;

إرسال تعليق

أحدث أقدم