Truss_All_Methods.m

%% Reliability analysis of the 23-bar truss with five fundamental methods
%  Companion code of the tutorial paper. Requires FEA_truss.m and the
%  PlaneTruss*.m functions in the same folder, and the Statistics and
%  Machine Learning Toolbox for the AK-MCS section (fitrgp).
%
%  Random variables (10 in total)
%    x(1:6)  loads P1 to P6          Gumbel     mean 5e4 N      std 7.5e3 N
%    x(7)    area A1                 Lognormal  mean 2e-3 m^2   std 2e-4 m^2
%    x(8)    area A2                 Lognormal  mean 1e-3 m^2   std 1e-4 m^2
%    x(9)    modulus E1              Lognormal  mean 2.1e11 Pa  std 2.1e10 Pa
%    x(10)   modulus E2              Lognormal  mean 2.1e11 Pa  std 2.1e10 Pa
%  Limit state function
%    g(x) = 0.11 - abs(delta(x)),  delta = mid-node deflection from FEA_truss

clc; clear; close all;
rng(1)

%% Isoprobabilistic transformations of Table 1 (U-space to X-space)
mP = 5e4;  sP = 7.5e3;
bG = sP*sqrt(6)/pi;                       % Gumbel scale, 1/alpha
vG = mP - 0.5772156649*bG;                % Gumbel location
lnp = @(m,s) deal(log(m)-0.5*log(1+(s/m)^2), sqrt(log(1+(s/m)^2)));
[lA1, zA1] = lnp(2e-3, 2e-4);
[lA2, zA2] = lnp(1e-3, 1e-4);
[lE,  zE ] = lnp(2.1e11, 2.1e10);

U2X = @(U) [ vG - bG*log(-log(normcdf(U(:,1:6)))), ...
             exp(lA1 + zA1*U(:,7)), exp(lA2 + zA2*U(:,8)), ...
             exp(lE  + zE *U(:,9)), exp(lE  + zE *U(:,10)) ];

gX = @(X) 0.11 - abs(arrayfun(@(r) FEA_truss(X(r,:)), (1:size(X,1))'));
gU = @(U) gX(U2X(U));
n  = 10;

%% 1. Crude Monte Carlo simulation (reference solution)
N = 2e6;                                   % Reduce N for a quick test run
U = randn(N, n);
G = gU(U);
Pf_mcs  = mean(G <= 0);
beta_mcs = -norminv(Pf_mcs);
CoV_mcs  = sqrt((1 - Pf_mcs)/(N*Pf_mcs));
fprintf('MCS     Pf = %.3e  beta = %.3f  CoV = %.3f  calls = %d\n', ...
        Pf_mcs, beta_mcs, CoV_mcs, N)

%% 2. FORM with the HLRF algorithm
h = 1e-4;  u = zeros(1, n);  Ncall = 0;
for k = 1:100
    gk = gU(u);  Ncall = Ncall + 1;
    grad = zeros(1, n);
    for i = 1:n                            % Central finite differences
        e = zeros(1, n);  e(i) = h;
        grad(i) = (gU(u+e) - gU(u-e))/(2*h);
        Ncall = Ncall + 2;
    end
    alpha = -grad/norm(grad);
    beta  = alpha*u' + gk/norm(grad);
    unew  = alpha*beta;
    if norm(unew - u) < 1e-6, u = unew; break, end
    u = unew;
end
ustar = u;
Pf_form = normcdf(-beta);
fprintf('FORM    Pf = %.3e  beta = %.3f  calls = %d\n', Pf_form, beta, Ncall)

%% 3. Importance sampling centered at the design point
%  The analysis is repeated 10 times and the averaged result is reported
Nis = 5000;  Nruns = 10;
Pf_is_all = zeros(Nruns, 1);
for r = 1:Nruns
    rng(r)
    Uis = randn(Nis, n) + ustar;
    Gis = gU(Uis);
    w   = exp(-0.5*sum(Uis.^2, 2))./exp(-0.5*sum((Uis - ustar).^2, 2));
    Pf_is_all(r) = mean((Gis <= 0).*w);
end
Pf_is   = mean(Pf_is_all);                 % Average of the 10 runs
CoV_is  = std(Pf_is_all)/Pf_is;            % Variability across the runs
fprintf('IS      Pf = %.3e  beta = %.3f  CoVruns = %.3f  calls/run = %d\n', ...
        Pf_is, -norminv(Pf_is), CoV_is, Nis + Ncall)

%% 4. Subset simulation with the modified Metropolis algorithm
%  The analysis is repeated 10 times and the averaged result is reported
Nss = 1000;  p0 = 0.1;  Nruns = 10;
Pf_ss_all = zeros(Nruns, 1);
for r = 1:Nruns
rng(r)
U = randn(Nss, n);  G = gU(U);  Pf_ss = 1;  Nc2 = Nss;
for level = 1:20
    [Gs, idx] = sort(G);
    Nc = p0*Nss;  bthr = Gs(Nc);
    if bthr <= 0, Pf_ss = Pf_ss*mean(G <= 0); break, end
    Pf_ss = Pf_ss*p0;
    seeds = U(idx(1:Nc), :);  Gseed = Gs(1:Nc);
    Unew = zeros(Nss, n);  Gnew = zeros(Nss, 1);  t = 0;
    for s = 1:Nc
        cur = seeds(s, :);  gc = Gseed(s);
        for step = 1:Nss/Nc
            cand = cur + randn(1, n);
            r    = normpdf(cand)./normpdf(cur);
            acc  = rand(1, n) < min(1, r);
            prop = cur;  prop(acc) = cand(acc);
            gp   = gU(prop);  Nc2 = Nc2 + 1;
            if gp <= bthr, cur = prop;  gc = gp; end
            t = t + 1;  Unew(t, :) = cur;  Gnew(t) = gc;
        end
    end
    U = Unew;  G = Gnew;
end
Pf_ss_all(r) = Pf_ss;
end
Pf_ss   = mean(Pf_ss_all);                 % Average of the 10 runs
CoV_ss  = std(Pf_ss_all)/Pf_ss;            % Variability across the runs
fprintf('SS      Pf = %.3e  beta = %.3f  CoVruns = %.3f  calls/run = %d\n', ...
        Pf_ss, -norminv(Pf_ss), CoV_ss, Nc2)

%% 5. AK-MCS with the U learning function
Nmc = 1e5;
S   = randn(Nmc, n);
n0  = 15;
id  = randperm(Nmc, n0);
D   = S(id, :);  Y = gU(D);  Nak = n0;
for iter = 1:500
    mdl = fitrgp(D, Y, 'KernelFunction', 'ardsquaredexponential', ...
                 'Standardize', true);
    [mg, sg] = predict(mdl, S);
    sg = max(sg, 1e-12);
    Ufun = abs(mg)./sg;
    Pf_ak = mean(mg <= 0);
    if min(Ufun) >= 2, break, end
    [~, j] = min(Ufun);
    D = [D; S(j, :)];  Y = [Y; gU(S(j, :))];  Nak = Nak + 1;
end
fprintf('AK-MCS  Pf = %.3e  beta = %.3f  calls = %d\n', ...
        Pf_ak, -norminv(Pf_ak), Nak)