结构可靠度

%% 谭斌，圆形截面钢拉杆可靠度分析
clc; clear; close all;
mu_D = 30;      % 直径均值 (mm)
sigma_D = 3;    % 直径标准差 (mm)
mu_r = 290;     % 强度均值 (N/mm^2)
sigma_r = 25;   % 强度标准差 (N/mm^2)
F = 1e5;        % 拉力 (N)

%% 中心点法计算
fprintf('中心点法计算\n');
mu_Z1 = (pi * mu_D^2 / 4) * mu_r - F;
dZ1_dr = pi * mu_D^2 / 4;      % ∂Z1/∂r
dZ1_dD = pi * mu_D * mu_r / 2;  % ∂Z1/∂D
sigma_Z1 = sqrt((dZ1_dr * sigma_r)^2 + (dZ1_dD * sigma_D)^2);
beta1_center = mu_Z1 / sigma_Z1;
R1_center = normcdf(beta1_center);

fprintf('功能函数1: Z = (πD²/4)*r - F\n');
fprintf('  μ_Z1 = %.2f N\n', mu_Z1);
fprintf('  σ_Z1 = %.2f N\n', sigma_Z1);
fprintf('  β1 = %.4f\n', beta1_center);
fprintf('  R1 = %.6f\n\n', R1_center);

mu_Z2 = mu_r - 4 * F / (pi * mu_D^2);
dZ2_dr = 1;  % ∂Z2/∂r
dZ2_dD = 8 * F / (pi * mu_D^3);  % ∂Z2/∂D
sigma_Z2 = sqrt((dZ2_dr * sigma_r)^2 + (dZ2_dD * sigma_D)^2);
beta2_center = mu_Z2 / sigma_Z2;
R2_center = normcdf(beta2_center);

fprintf('功能函数2: Z = r - 4F/(πD²)\n');
fprintf('  μ_Z2 = %.3f N/mm²\n', mu_Z2);
fprintf('  σ_Z2 = %.3f N/mm²\n', sigma_Z2);
fprintf('  β2 = %.4f\n', beta2_center);
fprintf('  R2 = %.6f\n\n', R2_center);

%% JC法计算
fprintf('JC法计算 \n');
tol = 1e-6;
max_iter = 50;

fprintf('使用功能函数1进行JC法迭代：\n');
beta_old = 0;
x = [mu_r, mu_D];  % [r, D]
mu = [mu_r, mu_D];
sigma = [sigma_r, sigma_D];
for iter = 1:max_iter
    r = x(1);
    D = x(2);
    Z = (pi * D^2 / 4) * r - F;
    % 梯度计算
    dZ_dr = pi * D^2 / 4;      % ∂Z/∂r
    dZ_dD = pi * D * r / 2;    % ∂Z/∂D
    grad = [dZ_dr, dZ_dD];
    % 灵敏度系数α
    sigma_grad = grad .* sigma;
    norm_sigma_grad = norm(sigma_grad);
    alpha = -sigma_grad / norm_sigma_grad;
    % 可靠指标β
    mu_Z = Z + grad * (mu - x)';
    sigma_Z = norm_sigma_grad;
    beta = mu_Z / sigma_Z;
    % 更新设计点
    x = mu + beta * (sigma .* alpha);
    % 检查收敛
    if iter > 1 && abs(beta - beta_old) < tol
        fprintf('  迭代 %d: β=%.6f, 收敛！\n', iter, beta);
        break;
    end
    beta_old = beta;
    if iter <= 5  % 显示前5次迭代
        fprintf('  迭代 %d: β=%.6f, r*=%.3f, D*=%.3f\n', iter, beta, x(1), x(2));
    end
    if iter == max_iter
        fprintf('  警告：达到最大迭代次数！\n');
    end
end

beta_JC = beta;
R_JC = normcdf(beta_JC);
design_point = x;

fprintf('\nJC法最终结果：\n');
fprintf('  可靠指标 β = %.6f\n', beta_JC);
fprintf('  可靠度 R = %.6f\n', R_JC);
fprintf('  设计点：r* = %.3f N/mm², D* = %.3f mm\n',design_point(1), design_point(2));

中心点法计算
功能函数1: Z = (πD²/4)*r - F
  μ_Z1 = 104988.92 N
  σ_Z1 = 44644.13 N
  β1 = 2.3517
  R1 = 0.990656

功能函数2: Z = r - 4F/(πD²)
  μ_Z2 = 148.529 N/mm²
  σ_Z2 = 37.757 N/mm²
  β2 = 3.9339
  R2 = 0.999958

JC法计算 
使用功能函数1进行JC法迭代：
  迭代 1: β=2.351685, r*=266.728, D*=23.521
  迭代 2: β=2.852905, r*=265.402, D*=21.966
  迭代 3: β=2.872293, r*=266.589, D*=21.854
  迭代 4: β=2.872213, r*=266.790, D*=21.846
  迭代 5: β=2.872212, r*=266.814, D*=21.845
  迭代 6: β=2.872212, 收敛！

JC法最终结果：
  可靠指标 β = 2.872212
  可靠度 R = 0.997962
  设计点：r* = 266.816 N/mm², D* = 21.845 mm

%% 谭斌，粒状土抗剪强度可靠度分析（中心点法）
clc; clear; close all;

%% 数据输入
tau = 55;           % 剪切应力 (kPa)
mu_w = 110;         % 法向应力均值 (kPa)
sigma_w = 25;       % 法向应力标准差 (kPa)
mu_phi_deg = 35;    % 摩擦角均值 (度)
sigma_phi_deg = 5;  % 摩擦角标准差 (度)
mu_phi = deg2rad(mu_phi_deg);       % 弧度
sigma_phi = deg2rad(sigma_phi_deg); % 弧度

fprintf('中心点法计算\n');
Z_mu = mu_w * tan(mu_phi) - tau;
dZ_dw = tan(mu_phi);
dZ_dphi = mu_w * (1/cos(mu_phi)^2);
mu_Z = Z_mu;
sigma_Z = sqrt((dZ_dw * sigma_w)^2 + (dZ_dphi * sigma_phi)^2);
% 计算可靠指标 β
beta = mu_Z / sigma_Z;
% 计算失效概率 p_f
p_f = normcdf(-beta);
R = normcdf(beta);

fprintf('失效概率与可靠度：\n');
fprintf('  可靠指标 β = %.4f\n', beta);
fprintf('  可靠度 R = Φ(β) = Φ(%.4f) = %.6f\n', beta, R);
fprintf('  失效概率 p_f = Φ(-β) = Φ(%.4f) = %.6f\n', -beta, p_f);
fprintf('  验证：p_f + R = %.6f + %.6f = %.6f ≈ 1\n', p_f, R, p_f+R);

中心点法计算
失效概率与可靠度：
  可靠指标 β = 0.9742
  可靠度 R = Φ(β) = Φ(0.9742) = 0.835009
  失效概率 p_f = Φ(-β) = Φ(-0.9742) = 0.164991
  验证：p_f + R = 0.164991 + 0.835009 = 1.000000 ≈ 1
>> 

%% 谭斌 钢梁抗弯可靠度分析（JC法）
clc; clear; close all;
M = 138;            % 弯矩 (kN·m)
mu_W = 900e-6;      % 截面模量均值 (m³)
delta_W = 0.05;     % W的变异系数
mu_f = 265e6;       % 强度均值 (Pa)
delta_f = 0.1;      % f的变异系数
sigma_W = mu_W * delta_W;
sigma_Y = sqrt(log(1 + delta_f^2));
mu_Y = log(mu_f) - sigma_Y^2/2;

%% JC法迭代计算
fprintf('JC法迭代计算\n');
tol = 1e-6;
max_iter = 50;
x = [mu_f, mu_W];      % 设计点 [f, W]
mu = [mu_f, mu_W];
beta_old = 0;
for iter = 1:max_iter
    % 当前设计点
    f_current = x(1);
    W_current = x(2);
    M_Nm = M * 1000;
    % 功能函数和梯度
    Z = f_current * W_current - M_Nm;
    grad = [W_current, f_current];  % [∂Z/∂f, ∂Z/∂W]
    % 当量正态化（对对数正态变量f）
    u_f = (log(f_current) - mu_Y) / sigma_Y;
    F_f = normcdf(u_f);
    f_f = (1/(f_current*sigma_Y*sqrt(2*pi))) * exp(-0.5*u_f^2);
    sigma_f_eq = normpdf(norminv(F_f)) / f_f;
    mu_f_eq = f_current - sigma_f_eq * norminv(F_f);
    % 当量正态化参数
    mu_eq = [mu_f_eq, mu_W];
    sigma_eq = [sigma_f_eq, sigma_W];
    % 灵敏度系数和可靠指标
    sigma_grad = grad .* sigma_eq;
    norm_sigma_grad = norm(sigma_grad);
    alpha = -sigma_grad / norm_sigma_grad;
    mu_Z = Z + grad * (mu_eq - x)';
    sigma_Z = norm_sigma_grad;
    beta = mu_Z / sigma_Z;
    % 更新设计点
    x = mu_eq + beta * (sigma_eq .* alpha);
    fprintf('迭代 %2d: β=%.6f, f*=%.1f MPa, W*=%.3e m³\n', iter, beta, x(1)/1e6, x(2));

    % 检查收敛
    if iter > 1 && abs(beta - beta_old) < tol
        fprintf('\n收敛于迭代 %d 次\n', iter);
        break;
    end

    beta_old = beta;

    if iter == max_iter
        fprintf('警告：达到最大迭代次数！\n');
    end
end

%% 最终结果
fprintf('\n最终结果\n');
fprintf('可靠指标: β = %.6f\n', beta);
fprintf('失效概率: P_f = Φ(-β) = %.6e\n', normcdf(-beta));
fprintf('可靠度: R = 1 - P_f = %.6f\n\n', 1 - normcdf(-beta));
fprintf('设计点（最可能失效点）：\n');
fprintf('  钢材强度 f* = %.1f MPa\n', x(1)/1e6);
fprintf('  截面模量 W* = %.3e m³\n', x(2));

JC法迭代计算
迭代  1: β=3.731892, f*=175.5 MPa, W*=8.247e-04 m³
迭代  2: β=4.792922, f*=173.3 MPa, W*=7.965e-04 m³
迭代  3: β=4.796909, f*=173.9 MPa, W*=7.936e-04 m³
迭代  4: β=4.796847, f*=174.0 MPa, W*=7.933e-04 m³
迭代  5: β=4.796846, f*=174.0 MPa, W*=7.933e-04 m³

收敛于迭代 5 次

最终结果
可靠指标: β = 4.796846
失效概率: P_f = Φ(-β) = 8.059160e-07
可靠度: R = 1 - P_f = 0.999999

设计点（最可能失效点）：
  钢材强度 f* = 174.0 MPa
  截面模量 W* = 7.933e-04 m³
>> 

%% 谭斌 钢筋混凝土受压短柱可靠度设计
clc; clear; close all;
fprintf('钢筋混凝土受压短柱可靠度设计（JC法\n');
beta_target = 3.7;
mu_G = 636;      % kN
sigma_G = 0.07;  % kN
mu_Q = 840;      % kN
delta_Q = 0.29;
sigma_Q = mu_Q * delta_Q;  % kN
delta_R = 0.18;

%% 改进的JC法迭代求解μ_R
fprintf('JC法迭代求解μ_R\n');
alpha_Q = 1.28255 / sigma_Q;
u_Q = mu_Q - 0.57722 / alpha_Q;
tol = 1e-4;
max_iter = 50;
% 确定μ_R的搜索范围
mu_R_low = mu_G + mu_Q;  % 下限：至少大于恒载+活载
mu_R_high = 2 * (mu_G + mu_Q) * (1 + 5*beta_target*delta_R);  % 上限
fprintf('搜索范围: μ_R ∈ [%.0f, %.0f] kN\n\n', mu_R_low, mu_R_high);

% 记录迭代过程
mu_R_history = [];
beta_history = [];

for iter = 1:max_iter
    % 当前μ_R（使用二分法）
    mu_R_current = (mu_R_low + mu_R_high) / 2;

    % 对于给定的μ_R，使用JC法计算实际β
    sigma_R = mu_R_current * delta_R;

    % 对数正态分布参数
    sigma_Y_R = sqrt(log(1 + delta_R^2));
    mu_Y_R = log(mu_R_current) - sigma_Y_R^2/2;

    % 初始化设计点
    x = [mu_R_current, mu_G, mu_Q];

    % JC法内层迭代（设计点迭代）
    for inner = 1:10
        R = x(1);
        G = x(2);
        Q = x(3);

        % 功能函数
        Z = R - G - Q;
        grad = [1, -1, -1];

        % 当量正态化：R（对数正态）
        u_R = (log(R) - mu_Y_R) / sigma_Y_R;
        F_R = normcdf(u_R);
        f_R = (1/(R*sigma_Y_R*sqrt(2*pi))) * exp(-0.5*u_R^2);
        sigma_R_eq = normpdf(norminv(F_R)) / f_R;
        mu_R_eq = R - sigma_R_eq * norminv(F_R);

        % 当量正态化：Q（极值I型）
        F_Q = exp(-exp(-alpha_Q*(Q - u_Q)));
        f_Q = alpha_Q * exp(-alpha_Q*(Q - u_Q)) * F_Q;
        sigma_Q_eq = normpdf(norminv(F_Q)) / f_Q;
        mu_Q_eq = Q - sigma_Q_eq * norminv(F_Q);

        % G是正态分布，参数不变
        mu_G_eq = mu_G;
        sigma_G_eq = sigma_G;

        % 当量正态化参数
        mu_eq = [mu_R_eq, mu_G_eq, mu_Q_eq];
        sigma_eq = [sigma_R_eq, sigma_G_eq, sigma_Q_eq];

        % 计算β
        sigma_grad = grad .* sigma_eq;
        norm_sigma_grad = norm(sigma_grad);
        alpha = -sigma_grad / norm_sigma_grad;

        mu_Z = Z + grad * (mu_eq - x)';
        sigma_Z = norm_sigma_grad;
        beta_current = mu_Z / sigma_Z;

        % 更新设计点
        x_new = mu_eq + beta_current * (sigma_eq .* alpha);

        % 检查内层收敛
        if norm(x_new - x) < 1e-6
            x = x_new;
            break;
        end
        x = x_new;
    end

    % 记录结果
    mu_R_history(end+1) = mu_R_current;
    beta_history(end+1) = beta_current;

    fprintf('迭代 %2d: μ_R = %.1f kN, β = %.6f\n', iter, mu_R_current, beta_current);

    % 检查是否达到目标β
    if abs(beta_current - beta_target) < tol
        fprintf('\n收敛于迭代 %d 次\n', iter);
        mu_R_final = mu_R_current;
        break;
    end

    % 根据β的差异调整搜索区间
    if beta_current < beta_target
        mu_R_low = mu_R_current;  % β太小，需要增加μ_R
    else
        mu_R_high = mu_R_current;  % β太大，需要减小μ_R
    end

    if iter == max_iter
        fprintf('\n达到最大迭代次数，取最接近的值\n');
        % 找到最接近目标β的μ_R
        [~, idx] = min(abs(beta_history - beta_target));
        mu_R_final = mu_R_history(idx);
    end
end

%% 输出最终结果
fprintf('\n最终结果\n');
fprintf('当β_target = %.1f时，所需的抗力平均值为：\n', beta_target);
fprintf('μ_R = %.1f kN\n\n', mu_R_final);

% 计算最终设计点
sigma_R_final = mu_R_final * delta_R;
sigma_Y_R = sqrt(log(1 + delta_R^2));
mu_Y_R = log(mu_R_final) - sigma_Y_R^2/2;

% 使用JC法计算设计点
x_final = [mu_R_final, mu_G, mu_Q];
for inner = 1:10
    R = x_final(1); G = x_final(2); Q = x_final(3);

    % 当量正态化
    u_R = (log(R) - mu_Y_R) / sigma_Y_R;
    F_R = normcdf(u_R);
    f_R = (1/(R*sigma_Y_R*sqrt(2*pi))) * exp(-0.5*u_R^2);
    sigma_R_eq = normpdf(norminv(F_R)) / f_R;
    mu_R_eq = R - sigma_R_eq * norminv(F_R);

    F_Q = exp(-exp(-alpha_Q*(Q - u_Q)));
    f_Q = alpha_Q * exp(-alpha_Q*(Q - u_Q)) * F_Q;
    sigma_Q_eq = normpdf(norminv(F_Q)) / f_Q;
    mu_Q_eq = Q - sigma_Q_eq * norminv(F_Q);

    mu_eq = [mu_R_eq, mu_G, mu_Q_eq];
    sigma_eq = [sigma_R_eq, sigma_G, sigma_Q_eq];

    Z = R - G - Q;
    grad = [1, -1, -1];
    sigma_grad = grad .* sigma_eq;
    norm_sigma_grad = norm(sigma_grad);
    alpha = -sigma_grad / norm_sigma_grad;

    mu_Z = Z + grad * (mu_eq - x_final)';
    sigma_Z = norm_sigma_grad;
    beta_final = mu_Z / sigma_Z;

    x_final = mu_eq + beta_final * (sigma_eq .* alpha);
end

fprintf('设计点（最可能失效点）：\n');
fprintf('  R* = %.1f kN\n', x_final(1));
fprintf('  G* = %.2f kN\n', x_final(2));
fprintf('  Q* = %.1f kN\n', x_final(3));
fprintf('  验证: R* - G* - Q* = %.2e ≈ 0\n', x_final(1)-x_final(2)-x_final(3));
fprintf('  实际可靠指标: β = %.6f\n\n', beta_final);

%% 计算受压钢筋面积
fprintf('计算受压钢筋面积\n');
b = 400;    % mm
h = 400;    % mm
K_R = 1.33;
f_c = 14.3;  % C30混凝土 (N/mm²)
f_y = 360;   % HRB400钢筋 (N/mm²)

% 计算设计抗力和钢筋面积
R_d = mu_R_final / K_R;
A_c = b * h;
N_c = 0.9 * f_c * A_c / 1000;
N_s = R_d - N_c;
A_s = max(0, N_s * 1000 / (0.9 * f_y));

fprintf('截面尺寸: %d×%d mm\n', b, h);
fprintf('设计抗力: R_d = μ_R/K_R = %.1f/%.2f = %.1f kN\n', mu_R_final, K_R, R_d);
fprintf('混凝土承载力: N_c = 0.9×%.1f×%d/1000 = %.1f kN\n', f_c, A_c, N_c);
fprintf('需要钢筋承担的荷载: N_s = %.1f - %.1f = %.1f kN\n', R_d, N_c, N_s);

if A_s > 0
    fprintf('受压钢筋面积: A_s = %.1f×1000/(0.9×%.0f) = %.0f mm²\n', N_s, f_y, A_s);
    % 钢筋配置
    d = 20;  % 假设直径20mm
    A_s_one = pi * d^2 / 4;
    n = ceil(A_s / A_s_one);
    fprintf('可选择 %dΦ%d (A_s = %d×%.1f = %.0f mm²)\n', n, d, n, A_s_one, n*A_s_one);
    fprintf('配筋率: ρ = %.2f%%\n', (n*A_s_one)/A_c*100);
else
    fprintf('按构造配筋: A_s_min = 0.006×%d = %.0f mm²\n', A_c, 0.006*A_c);
end

钢筋混凝土受压短柱可靠度设计（JC法
JC法迭代求解μ_R
搜索范围: μ_R ∈ [1476, 12782] kN

迭代  1: μ_R = 7129.1 kN, β = 5.697893
迭代  2: μ_R = 4302.5 kN, β = 4.010849
迭代  3: μ_R = 2889.3 kN, β = 2.651828
迭代  4: μ_R = 3595.9 kN, β = 3.407055
迭代  5: μ_R = 3949.2 kN, β = 3.723581
迭代  6: μ_R = 3772.6 kN, β = 3.569393
迭代  7: μ_R = 3860.9 kN, β = 3.647448
迭代  8: μ_R = 3905.1 kN, β = 3.685748
迭代  9: μ_R = 3927.1 kN, β = 3.704722
迭代 10: μ_R = 3916.1 kN, β = 3.695249
迭代 11: μ_R = 3921.6 kN, β = 3.699989

收敛于迭代 11 次

最终结果
当β_target = 3.7时，所需的抗力平均值为：
μ_R = 3921.6 kN

设计点（最可能失效点）：
  R* = 2595.1 kN
  G* = 636.00 kN
  Q* = 1959.1 kN
  验证: R* - G* - Q* = -6.82e-13 ≈ 0
  实际可靠指标: β = 3.699989

计算受压钢筋面积
截面尺寸: 400×400 mm
设计抗力: R_d = μ_R/K_R = 3921.6/1.33 = 2948.6 kN
混凝土承载力: N_c = 0.9×14.3×160000/1000 = 2059.2 kN
需要钢筋承担的荷载: N_s = 2948.6 - 2059.2 = 889.4 kN
受压钢筋面积: A_s = 889.4×1000/(0.9×360) = 2745 mm²
可选择 9Φ20 (A_s = 9×314.2 = 2827 mm²)
配筋率: ρ = 1.77%
>> 

clear; clc;

%% 给定参数
fprintf('===== 蒙特卡洛法计算结构构件可靠度 =====\n');

% 各随机变量的均值和标准差
mu_x1 = 25;
sigma_x1 = 5.75;

mu_x2 = 0.0113;
sigma_x2 = 0.00339;

mu_x3 = 0.0006;
sigma_x3 = 0.00018;

mu_x4 = 0;
sigma_x4 = 0.1;

% 极限状态方程：Z = 1 - x1*x2 - x1^2*x3 - x4

%% 蒙特卡洛模拟参数设置
N = 1000000;  % 模拟次数
failure_count = 0;  % 失效次数计数器

fprintf('开始蒙特卡洛模拟，样本量：%d\n', N);

% 为了监控进度，设置进度显示间隔
progress_interval = floor(N/10);

%% 生成随机样本并进行模拟
for i = 1:N
    % 1. 生成x1（对数正态分布）
    % 计算对数正态分布的对数空间参数
    mu_lnx1 = log(mu_x1^2 / sqrt(sigma_x1^2 + mu_x1^2));
    sigma_lnx1 = sqrt(log(1 + (sigma_x1^2 / mu_x1^2)));
    x1 = lognrnd(mu_lnx1, sigma_lnx1);
    % 2. 生成x2（正态分布）
    x2 = normrnd(mu_x2, sigma_x2);
    % 3. 生成x3（正态分布）
    x3 = normrnd(mu_x3, sigma_x3);
    % 4. 生成x4（正态分布）
    x4 = normrnd(mu_x4, sigma_x4);

    % 5. 计算极限状态函数值
    Z = 1 - x1*x2 - x1^2*x3 - x4;
    % 6. 判断是否失效
    if Z <= 0
        failure_count = failure_count + 1;
    end
    % 显示进度
    if mod(i, progress_interval) == 0
        fprintf('  进度：%d/%d (%.0f%%)\n', i, N, i/N*100);
    end
end

%% 计算失效概率和可靠指标
Pf = failure_count / N;  % 失效概率
beta = -norminv(Pf);     % 可靠指标

%% 结果输出
fprintf('\n===== 模拟结果 =====\n');
fprintf('总样本数：%d\n', N);
fprintf('失效次数：%d\n', failure_count);
fprintf('失效概率 Pf = %.6f\n', Pf);
fprintf('可靠指标 β = %.4f\n', beta);
fprintf('可靠度 = %.6f\n', 1-Pf);

===== 蒙特卡洛法计算结构构件可靠度 =====
开始蒙特卡洛模拟，样本量：1000000
  进度：100000/1000000 (10%)
  进度：200000/1000000 (20%)
  进度：300000/1000000 (30%)
  进度：400000/1000000 (40%)
  进度：500000/1000000 (50%)
  进度：600000/1000000 (60%)
  进度：700000/1000000 (70%)
  进度：800000/1000000 (80%)
  进度：900000/1000000 (90%)
  进度：1000000/1000000 (100%)

===== 模拟结果 =====
总样本数：1000000
失效次数：136681
失效概率 Pf = 0.136681
可靠指标 β = 1.0954
可靠度 = 0.863319

%% 钢筋混凝土框架底层柱可靠度分析（蒙特卡洛法）
clear; clc;

%% 已知参数

% R: 抗力，对数正态分布
mu_R = 8.5;
delta_R = 0.15;
sigma_R = mu_R * delta_R;

% 对数正态分布的对数空间参数
sigma_lnR = sqrt(log(1 + delta_R^2));
mu_lnR = log(mu_R) - 0.5 * sigma_lnR^2;

% S_G: 永久荷载效应，正态分布
mu_SG = 2.3;
delta_SG = 0.05;
sigma_SG = mu_SG * delta_SG;

% S_Q1, S_Q2, S_Q3: 可变荷载效应，极值I型分布
% S_Q1
mu_SQ1 = 0.8;
delta_SQ1 = 0.415;
sigma_SQ1 = mu_SQ1 * delta_SQ1;

% S_Q2
mu_SQ2 = 0.6;
delta_SQ2 = 0.415;
sigma_SQ2 = mu_SQ2 * delta_SQ2;

% S_Q3
mu_SQ3 = 0.5;
delta_SQ3 = 0.29;
sigma_SQ3 = mu_SQ3 * delta_SQ3;

% 极值I型分布参数计算
% 对于极值I型分布：CDF = exp(-exp(-α(x-u)))，其中：
% α = π/(√6σ) ≈ 1.28255/σ
% u = μ - 0.57722/α

alpha_SQ1 = 1.28255 / sigma_SQ1;
u_SQ1 = mu_SQ1 - 0.57722 / alpha_SQ1;

alpha_SQ2 = 1.28255 / sigma_SQ2;
u_SQ2 = mu_SQ2 - 0.57722 / alpha_SQ2;

alpha_SQ3 = 1.28255 / sigma_SQ3;
u_SQ3 = mu_SQ3 - 0.57722 / alpha_SQ3;

%% 2. 蒙特卡洛模拟参数
N = 1000000;  % 模拟次数
fprintf('\n开始蒙特卡洛模拟，样本量：N = %d\n\n', N);

% 预分配数组
Z_values = zeros(N, 1);
failure_count = 0;

% 进度显示
progress_interval = floor(N/10);
fprintf('模拟进度：\n');

%% 3. 蒙特卡洛模拟
start_time = tic;

for i = 1:N
    % 生成R的样本（对数正态分布）
    R = lognrnd(mu_lnR, sigma_lnR);
    % 生成S_G的样本（正态分布）
    SG = mu_SG + sigma_SG * randn;
    % 生成S_Q1, S_Q2, S_Q3的样本（极值I型分布）
    % 极值I型分布：使用逆变换法生成
    U1 = rand;
    SQ1 = u_SQ1 - (1/alpha_SQ1) * log(-log(U1));
    U2 = rand;
    SQ2 = u_SQ2 - (1/alpha_SQ2) * log(-log(U2));
    U3 = rand;
    SQ3 = u_SQ3 - (1/alpha_SQ3) * log(-log(U3));
    % 计算功能函数值
    Z = R - SG - SQ1 - SQ2 - SQ3;
    Z_values(i) = Z;
    % 判断是否失效
    if Z <= 0
        failure_count = failure_count + 1;
    end
    % 显示进度
    if mod(i, progress_interval) == 0
        fprintf('  进度：%d/%d (%.0f%%)\n', i, N, i/N*100);
    end
end

elapsed_time = toc(start_time);
fprintf('模拟完成！耗时：%.2f 秒\n\n', elapsed_time);

%% 4. 计算失效概率和可靠指标
Pf = failure_count / N;           % 失效概率
R_val = 1 - Pf;                   % 可靠度
beta = -norminv(Pf);              % 可靠指标

%% 5. 结果输出
fprintf('===== 模拟结果 =====\n');
fprintf('总样本数：N = %d\n', N);
fprintf('失效次数：N_f = %d\n', failure_count);
fprintf('失效概率：Pf = %d/%d = %.6e\n', failure_count, N, Pf);
fprintf('可靠度：R = 1 - Pf = %.6f\n', R_val);
fprintf('可靠指标：β = -Φ⁻¹(Pf) = %.4f\n\n', beta);

% 功能函数Z的统计特性
Z_mean = mean(Z_values);
Z_std = std(Z_values);
Z_min = min(Z_values);
Z_max = max(Z_values);

fprintf('功能函数Z的统计特性：\n');
fprintf('  均值 μ_Z = %.4f\n', Z_mean);
fprintf('  标准差 σ_Z = %.4f\n', Z_std);
fprintf('  最小值 Z_min = %.4f\n', Z_min);
fprintf('  最大值 Z_max = %.4f\n', Z_max);
fprintf('  安全裕度 μ_Z/σ_Z = %.4f\n', Z_mean/Z_std);

开始蒙特卡洛模拟，样本量：N = 1000000

模拟进度：
  进度：100000/1000000 (10%)
  进度：200000/1000000 (20%)
  进度：300000/1000000 (30%)
  进度：400000/1000000 (40%)
  进度：500000/1000000 (50%)
  进度：600000/1000000 (60%)
  进度：700000/1000000 (70%)
  进度：800000/1000000 (80%)
  进度：900000/1000000 (90%)
  进度：1000000/1000000 (100%)
模拟完成！耗时：6.07 秒

===== 模拟结果 =====
总样本数：N = 1000000
失效次数：N_f = 121
失效概率：Pf = 121/1000000 = 1.210000e-04
可靠度：R = 1 - Pf = 0.999879
可靠指标：β = -Φ⁻¹(Pf) = 3.6706

功能函数Z的统计特性：
  均值 μ_Z = 4.2995
  标准差 σ_Z = 1.3526
  最小值 Z_min = -1.2769
  最大值 Z_max = 12.5964
  安全裕度 μ_Z/σ_Z = 3.1788

%% 钢筋混凝土轴心受压短柱可靠度分析（修正版）
clear; clc; close all;

%% 1. 已知参数
fprintf('===== 钢筋混凝土轴心受压短柱可靠度分析 =====\n');

% R: 抗力，对数正态分布
mu_R = 4560;      % kN
sigma_R = 729.6;  % kN
delta_R = sigma_R / mu_R;

% 对数正态分布的对数空间参数
sigma_lnR = sqrt(log(1 + delta_R^2));
mu_lnR = log(mu_R) - 0.5 * sigma_lnR^2;

% N_G: 恒载，正态分布
mu_NG = 1159.1;   % kN
sigma_NG = 81.1;  % kN

% N_L: 活载，极值I型分布
mu_NL = 765.5;    % kN
sigma_NL = 222;   % kN

% 修正的极值I型分布参数（更精确的计算）
alpha_NL = pi/(sigma_NL*sqrt(6));  % 修正：使用精确公式
gamma = 0.5772156649;  % 欧拉常数
u_NL = mu_NL - gamma/alpha_NL;

fprintf('随机变量参数：\n');
fprintf('R (抗力，对数正态): μ=%.1f kN, σ=%.1f kN, δ=%.4f\n', mu_R, sigma_R, delta_R);
fprintf('   对数空间参数: μ_ln=%.6f, σ_ln=%.6f\n', mu_lnR, sigma_lnR);

fprintf('N_G (恒载，正态): μ=%.1f kN, σ=%.1f kN\n', mu_NG, sigma_NG);
fprintf('N_L (活载，极值I型): μ=%.1f kN, σ=%.1f kN\n', mu_NL, sigma_NL);
fprintf('   极值I型参数: α=%.6f, u=%.4f\n\n', alpha_NL, u_NL);

%% 2. JC法计算
fprintf('-------- JC法计算 --------\n');

tol = 1e-6;
max_iter = 50;
x = [mu_R, mu_NG, mu_NL];  % [R, N_G, N_L]
beta_old = 0;

fprintf('JC法迭代过程：\n');
for iter = 1:max_iter
    R = x(1);
    NG = x(2);
    NL = x(3);
    % 功能函数
    Z = R - NG - NL;
    grad = [1, -1, -1];
    % 当量正态化：R（对数正态）
    u_R = (log(R) - mu_lnR) / sigma_lnR;
    F_R = normcdf(u_R);
    f_R = (1/(R*sigma_lnR*sqrt(2*pi))) * exp(-0.5*u_R^2);
    sigma_R_eq = normpdf(norminv(F_R)) / f_R;
    mu_R_eq = R - sigma_R_eq * norminv(F_R);
    % 当量正态化：N_L（极值I型）
    F_NL = exp(-exp(-alpha_NL*(NL - u_NL)));
    if F_NL > 0 && F_NL < 1
        f_NL = alpha_NL * exp(-alpha_NL*(NL - u_NL)) * F_NL;
        sigma_NL_eq = normpdf(norminv(F_NL)) / f_NL;
        mu_NL_eq = NL - sigma_NL_eq * norminv(F_NL);
    else
        % 边界情况处理
        sigma_NL_eq = sigma_NL;  % 近似为正态分布
        mu_NL_eq = mu_NL;
    end
    % N_G是正态分布，参数不变
    mu_NG_eq = mu_NG;
    sigma_NG_eq = sigma_NG;
    % 当量正态化参数
    mu_eq = [mu_R_eq, mu_NG_eq, mu_NL_eq];
    sigma_eq = [sigma_R_eq, sigma_NG_eq, sigma_NL_eq];
    % 可靠指标计算
    sigma_grad = grad .* sigma_eq;
    norm_sigma_grad = norm(sigma_grad);
    alpha = -sigma_grad / norm_sigma_grad;
    mu_Z = Z + grad * (mu_eq - x)';
    sigma_Z = norm_sigma_grad;
    beta = mu_Z / sigma_Z;
    % 更新设计点
    x_new = mu_eq + beta * (sigma_eq .* alpha);
    fprintf('迭代 %2d: β=%.6f, R*=%.1f, N_G*=%.1f, N_L*=%.1f\n',iter, beta, x(1), x(2), x(3));
    % 检查收敛
    if iter > 1 && abs(beta - beta_old) < tol
        fprintf('\nJC法收敛于迭代 %d 次\n', iter);
        break;
    end
    beta_old = beta;
    x = x_new;
    if iter == max_iter
        fprintf('警告：JC法达到最大迭代次数！\n');
    end
end

beta_JC = beta;
Pf_JC = normcdf(-beta_JC);
design_point_JC = x;

fprintf('\nJC法最终结果：\n');
fprintf('  可靠指标: β_JC = %.6f\n', beta_JC);
fprintf('  失效概率: P_f = Φ(-β) = %.6e\n', Pf_JC);
fprintf('  可靠度: R = 1 - P_f = %.6f\n', 1-Pf_JC);
fprintf('  设计验算点:\n');
fprintf('    R* = %.1f kN\n', design_point_JC(1));
fprintf('    N_G* = %.1f kN\n', design_point_JC(2));
fprintf('    N_L* = %.1f kN\n', design_point_JC(3));
fprintf('  验证: Z(R*,N_G*,N_L*) = %.2e ≈ 0\n\n', design_point_JC(1)-design_point_JC(2)-design_point_JC(3));

%% 3. 蒙特卡罗法计算（修正版）
fprintf('-------- 蒙特卡罗法计算 --------\n');

N_MC = 1000000;
fprintf('蒙特卡罗模拟，样本量：N = %d\n\n', N_MC);

% 更高效的向量化实现
rng(42);  % 设置随机种子保证可重复性
fprintf('生成随机样本...\n');

% 对数正态分布抽样
U_R = randn(N_MC, 1);
R_samples = exp(mu_lnR + sigma_lnR * U_R);

% 正态分布抽样
NG_samples = mu_NG + sigma_NG * randn(N_MC, 1);

% 极值I型分布抽样（向量化）
U_NL = rand(N_MC, 1);
NL_samples = u_NL - (1/alpha_NL) * log(-log(U_NL));

% 计算功能函数值（向量化）
Z_samples = R_samples - NG_samples - NL_samples;

% 统计失效
failure_mask = Z_samples <= 0;
failure_count = sum(failure_mask);

% 计算结果
Pf_MC = failure_count / N_MC;
R_MC = 1 - Pf_MC;
beta_MC = -norminv(Pf_MC);

% 精度分析
if failure_count > 0
    cov_Pf = sqrt((1-Pf_MC)/(N_MC*Pf_MC));
    z_95 = 1.96;
    CI_lower = Pf_MC * (1 - z_95 * cov_Pf);
    CI_upper = Pf_MC * (1 + z_95 * cov_Pf);
else
    cov_Pf = 0;
    CI_lower = 0;
    CI_upper = 0;
end

fprintf('蒙特卡罗法结果：\n');
fprintf('  抽样次数: N = %d\n', N_MC);
fprintf('  失效次数: N_f = %d\n', failure_count);
fprintf('  失效概率: P_f = %.6e\n', Pf_MC);
fprintf('  可靠度: R = %.6f\n', R_MC);
fprintf('  可靠指标: β_MC = %.6f\n', beta_MC);
fprintf('  精度指标: COV(P_f) = %.4f\n', cov_Pf);
fprintf('  95%%置信区间: [%.6e, %.6e]\n', CI_lower, CI_upper);

% 统计特性
Z_mean = mean(Z_samples);
Z_std = std(Z_samples);
fprintf('  功能函数统计: μ_Z = %.2f, σ_Z = %.2f, μ_Z/σ_Z = %.4f',Z_mean, Z_std, Z_mean/Z_std);

===== 钢筋混凝土轴心受压短柱可靠度分析 =====
随机变量参数：
R (抗力，对数正态): μ=4560.0 kN, σ=729.6 kN, δ=0.1600
   对数空间参数: μ_ln=8.412439, σ_ln=0.158990
N_G (恒载，正态): μ=1159.1 kN, σ=81.1 kN
N_L (活载，极值I型): μ=765.5 kN, σ=222.0 kN
   极值I型参数: α=0.005777, u=665.5882

-------- JC法计算 --------
JC法迭代过程：
迭代  1: β=3.442367, R*=4560.0, N_G*=1159.1, N_L*=765.5
迭代  2: β=4.217026, R*=2120.9, N_G*=1188.9, N_L*=932.0
迭代  3: β=3.991285, R*=2632.8, N_G*=1221.9, N_L*=1410.9
迭代  4: β=3.960488, R*=2913.4, N_G*=1201.6, N_L*=1711.8
迭代  5: β=3.959016, R*=2996.0, N_G*=1195.4, N_L*=1800.6
迭代  6: β=3.958978, R*=3009.8, N_G*=1194.1, N_L*=1815.7
迭代  7: β=3.958977, R*=3011.8, N_G*=1193.9, N_L*=1818.0

JC法收敛于迭代 7 次

JC法最终结果：
  可靠指标: β_JC = 3.958977
  失效概率: P_f = Φ(-β) = 3.763578e-05
  可靠度: R = 1 - P_f = 0.999962
  设计验算点:
    R* = 3011.8 kN
    N_G* = 1193.9 kN
    N_L* = 1818.0 kN
  验证: Z(R*,N_G*,N_L*) = -4.55e-13 ≈ 0

-------- 蒙特卡罗法计算 --------
蒙特卡罗模拟，样本量：N = 1000000

生成随机样本...
蒙特卡罗法结果：
  抽样次数: N = 1000000
  失效次数: N_f = 40
  失效概率: P_f = 4.000000e-05
  可靠度: R = 0.999960
  可靠指标: β_MC = 3.944400
  精度指标: COV(P_f) = 0.1581
  95%置信区间: [2.760412e-05, 5.239588e-05]
  功能函数统计: μ_Z = 2635.88, σ_Z = 766.47, μ_Z/σ_Z = 3.4390>>