subset simulation.rar_subset simulation_失效分析_子集模拟_小概率_数值模拟

function [minimum,solu,numTol,gRe,zRe] = SSO(objFun, n, boundaries, paras) % Algorithm parameters N = paras.NumSam; p0 = paras.CondPro; nt = round(N*p0); ns = ceil(1/p0-1); % Distributional parameters for truncated normal distribution means = mean(boundaries,2); stds = abs(boundaries(:,1)-boundaries(:,2))/6; %% Step 1 (Monte Carlo simulation) % Generate samples of truncated normal variables z = []; for i=1:n U(i) = normcdf(boundaries(i,2),means(i),stds(i)); L(i) = normcdf(boundaries(i,1),means(i),stds(i)); uu = unifrnd(0,1,N,1)*(U(i)-L(i))+L(i); tempz = norminv(uu,means(i),stds(i)); z = [z,tempz]; end % Calculate objective function for i=1:N; g(i) = feval(objFun,z(i,:)); end; % Sort samples and and objevtice function values [gSort,index]=sort(g); Z = z(index,:); % Record the calculated results numTol = N; gRe = gSort(1); zRe =Z(1,:); %% Step 2 (Markov Chain Monte Carlo) % preparation for MCMC sigma = std(Z,0,1); pfss = 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,:); gseed = g(i); for j = 1:ns % A Markov chain for k = 1:n % Component by component % Generate an candidate done = false; while ~done u(k) = seed(k)+(3*rand()-1)*sigma(k); if u(k)>= boundaries(k,1) && u(k)<= boundaries(k,2) done = true; end end % Calculate the acceptance probability pdf2 = exp(-0.5*((seed(k)-means(k))/stds(k)).^2); pdf1 = exp(-0.5*((u(k)-means(k))/stds(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(objFun,w(len,:)); % Accept or reject if gTemp <= gSort(nt+1); g(len) = gTemp; else g(len) = gseed; w(len,:) = seed; end seed = w(len,:);gseed = g(len); 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 objective function values [gSort,index]=sort(g); Z = w(index,:); sigma = std(Z,0,1); % Record the calculated results gRe = [gRe;gSort(1)]; zRe = [zRe;Z(1,:)]; numTol = numTol+N-nt; % Terminate simulation max_sigma = max(sigma); if max_sigma < paras.Tol || iter >= paras.MaxTry stopFlag = 1; break; else pfss = pfss*nt/N; end; end; minimum = gSort(1); solu = Z(1,:);