结构动力学仿真-主题071-结构动力学中的不确定性分析

kkchenkx

441人浏览 · 2026-04-09 21:48:39

kkchenkx · 2026-04-09 21:48:39 发布

主题071：结构动力学中的不确定性分析

1. 引言

在实际的工程结构分析中，不确定性无处不在。材料参数的离散性、几何尺寸的制造误差、荷载的随机性以及边界条件的模糊性，都使得确定性分析方法难以准确描述结构的真实行为。结构动力学中的不确定性分析，旨在量化这些不确定性对结构动力响应的影响，为结构的安全评估和可靠性设计提供科学依据。

不确定性通常分为两类：偶然不确定性（Aleatory Uncertainty）和认知不确定性（Epistemic Uncertainty）。偶然不确定性源于系统的固有随机性，如材料性能的统计变异；认知不确定性则源于知识的不完备，如模型简化带来的误差。本章将重点讨论偶然不确定性的分析方法，包括概率方法、随机过程理论以及蒙特卡洛模拟技术。
在这里插入图片描述

2. 不确定性的数学描述

2.1 随机变量及其分布

在结构动力学中，不确定性参数通常用随机变量来描述。设 $X$ 是一个随机变量，其概率密度函数（PDF）为 $f_X(x)$ ，累积分布函数（CDF）为 $FX(x)=P(X≤x)F_X(x) = P(X \leq x)$ 。

常见的随机变量分布包括：

正态分布（高斯分布）：
$fX(x)=1σ2πexp⁡(−(x−μ)22σ2)f_X(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$

其中， $μ\mu$ 为均值， $σ\sigma$ 为标准差。正态分布适用于描述许多自然现象的变异，如材料强度、弹性模量等。

对数正态分布：
$fX(x)=1xσ2πexp⁡(−(ln⁡x−μ)22σ2),x>0f_X(x) = \frac{1}{x\sigma\sqrt{2\pi}} \exp\left(-\frac{(\ln x - \mu)^2}{2\sigma^2}\right), \quad x > 0$

对数正态分布适用于描述只能取正值的参数，如材料的弹性模量、密度等。

均匀分布：
$fX(x)={1b−a,a≤x≤b0,其他f_X(x) = \begin{cases} \frac{1}{b-a}, & a \leq x \leq b \\ 0, & \text{其他} \end{cases}$

均匀分布用于描述在某一区间内等可能取值的参数。

威布尔分布：
$fX(x)=kλ(xλ)k−1exp⁡(−(xλ)k),x≥0f_X(x) = \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1} \exp\left(-\left(\frac{x}{\lambda}\right)^k\right), \quad x \geq 0$

威布尔分布常用于描述疲劳寿命和强度分布。

2.2 随机过程

当不确定性随时间变化时，需要用随机过程来描述。随机过程 ${X(t),t∈T}\{X(t), t \in T\}$ 是一族随机变量的集合，其中 $t$ 为时间参数。

平稳随机过程：
如果随机过程的统计特性不随时间平移而变化，即：
$\mu_X \quad \text{(常数)}$
$RX(t1,t2)=RX(τ),τ=t2−t1R_X(t_1, t_2) = R_X(\tau), \quad \tau = t_2 - t_1$

则称该过程为平稳随机过程。其中 $RX(τ)R_X(\tau)$ 为自相关函数。

高斯随机过程：
如果随机过程的任意有限维分布都是多维正态分布，则称为高斯随机过程。高斯过程完全由其均值函数和协方差函数确定。

白噪声过程：
白噪声是一种理想化的随机过程，其功率谱密度在所有频率上为常数：
$SXX(ω)=S0,−∞<ω<∞S_{XX}(\omega) = S_0, \quad -\infty < \omega < \infty$

自相关函数为狄拉克函数：
$RXX(τ)=2πS0δ(τ)R_{XX}(\tau) = 2\pi S_0 \delta(\tau)$

2.3 随机场

当不确定性在空间上变化时，需要用随机场来描述。随机场 ${X(x),x∈D}\{X(\mathbf{x}), \mathbf{x} \in D\}$ 是定义在空间域 $D$ 上的一族随机变量。

随机场的相关结构通常用相关函数来描述：
$R(x1,x2)=E[(X(x1)−μ(x1))(X(x2)−μ(x2))]R(\mathbf{x}_1, \mathbf{x}_2) = E[(X(\mathbf{x}_1) - \mu(\mathbf{x}_1))(X(\mathbf{x}_2) - \mu(\mathbf{x}_2))]$

常用的相关函数模型包括：

指数型相关函数：
$ρ(r)=exp⁡(−∣r∣b)\rho(r) = \exp\left(-\frac{|r|}{b}\right)$

高斯型相关函数：
$ρ(r)=exp⁡(−(rb)2)\rho(r) = \exp\left(-\left(\frac{r}{b}\right)^2\right)$

其中 $b$ 为相关长度，表征随机场的空间相关性尺度。

3. 结构参数不确定性分析

3.1 参数不确定性的来源

结构动力学中的参数不确定性主要来源于以下几个方面：

材料参数的不确定性：

弹性模量 $E$ ：受材料成分、加工工艺、温度等因素影响
密度 $ρ\rho$ ：受材料均匀性、孔隙率等因素影响
阻尼系数 $ξ\xi$ ：受材料内部摩擦、连接条件等因素影响

几何参数的不确定性：

截面尺寸：制造公差导致的尺寸偏差
结构长度：施工误差引起的长度变化
边界条件：实际约束与理想约束的差异

荷载参数的不确定性：

地震荷载：地震动强度、频谱特性的随机性
风荷载：风速、风向的随机波动
活荷载：人群、车辆荷载的随机分布

3.2 摄动方法

摄动方法是分析参数不确定性对结构响应影响的一种有效方法。基本思想是将不确定参数表示为确定性部分和随机扰动部分的叠加：

$X_0 + \epsilon X_1$

其中 $X_0$ 为均值， $ϵ\epsilon$ 为小参数， $X_1$ 为随机扰动。

一阶摄动法：

设结构的动力方程为：
$Mu¨+Cu˙+Ku=F(t)\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{F}(t)$

当质量矩阵、阻尼矩阵和刚度矩阵存在不确定性时，可以表示为：
$M=M0+ϵM1,C=C0+ϵC1,K=K0+ϵK1\mathbf{M} = \mathbf{M}_0 + \epsilon \mathbf{M}_1, \quad \mathbf{C} = \mathbf{C}_0 + \epsilon \mathbf{C}_1, \quad \mathbf{K} = \mathbf{K}_0 + \epsilon \mathbf{K}_1$

位移响应也进行摄动展开：
$u=u0+ϵu1+O(ϵ2)\mathbf{u} = \mathbf{u}_0 + \epsilon \mathbf{u}_1 + O(\epsilon^2)$

代入动力方程并比较 $ϵ\epsilon$ 的同次幂，得到：

零阶方程（确定性方程）：
$M0u¨0+C0u˙0+K0u0=F(t)\mathbf{M}_0\ddot{\mathbf{u}}_0 + \mathbf{C}_0\dot{\mathbf{u}}_0 + \mathbf{K}_0\mathbf{u}_0 = \mathbf{F}(t)$

一阶方程：
$M0u¨1+C0u˙1+K0u1=−M1u¨0−C1u˙0−K1u0\mathbf{M}_0\ddot{\mathbf{u}}_1 + \mathbf{C}_0\dot{\mathbf{u}}_1 + \mathbf{K}_0\mathbf{u}_1 = -\mathbf{M}_1\ddot{\mathbf{u}}_0 - \mathbf{C}_1\dot{\mathbf{u}}_0 - \mathbf{K}_1\mathbf{u}_0$

通过求解上述方程，可以得到响应的均值和方差。

响应统计量的计算：

响应的均值：
$E[u]≈u0E[\mathbf{u}] \approx \mathbf{u}_0$

响应的协方差：
$Cov(u)=E[(u−E[u])(u−E[u])T]≈ϵ2E[u1u1T]\text{Cov}(\mathbf{u}) = E[(\mathbf{u} - E[\mathbf{u}])(\mathbf{u} - E[\mathbf{u}])^T] \approx \epsilon^2 E[\mathbf{u}_1 \mathbf{u}_1^T]$

3.3 随机有限元方法

随机有限元方法（Stochastic Finite Element Method, SFEM）是有限元方法与概率理论的结合，用于分析具有随机参数的结构。

谱展开方法：

将随机场离散化为随机变量的函数。常用的方法包括：

Karhunen-Loève展开：
$X(x,ω)=∑i=1∞λiξi(ω)ϕi(x)X(\mathbf{x}, \omega) = \sum_{i=1}^{\infty} \sqrt{\lambda_i} \xi_i(\omega) \phi_i(\mathbf{x})$

其中 $λi\lambda_i$ 和 $ϕi(x)\phi_i(\mathbf{x})$ 是协方差核的特征值和特征函数， $ξi\xi_i$ 是互不相关的标准随机变量。

多项式混沌展开：
$Y(ω)=∑j=0PyjΨj(ξ)Y(\omega) = \sum_{j=0}^{P} y_j \Psi_j(\boldsymbol{\xi})$

其中 $Ψj\Psi_j$ 是多元正交多项式， $ξ\boldsymbol{\xi}$ 是标准随机变量向量。

4. 随机振动分析

4.1 线性系统的随机振动

对于线性时不变系统，在随机激励下的响应分析可以采用频域方法或时域方法。

频域分析方法：

设系统的频率响应函数为 $H(ω)H(\omega)$ ，输入随机过程的功率谱密度为 $SFF(ω)S_{FF}(\omega)$ ，则输出响应的功率谱密度为：

$SUU(ω)=∣H(ω)∣2SFF(ω)S_{UU}(\omega) = |H(\omega)|^2 S_{FF}(\omega)$

对于多自由度系统，频率响应函数矩阵为：
$H(ω)=(−ω2M+iωC+K)−1\mathbf{H}(\omega) = (-\omega^2 \mathbf{M} + i\omega \mathbf{C} + \mathbf{K})^{-1}$

响应功率谱密度矩阵为：
$SUU(ω)=H∗(ω)SFF(ω)HT(ω)\mathbf{S}_{UU}(\omega) = \mathbf{H}^*(\omega) \mathbf{S}_{FF}(\omega) \mathbf{H}^T(\omega)$

其中 $*$ 表示复共轭， $T$ 表示转置。

响应的统计量：

响应的均方值：
$E[U2]=∫−∞∞SUU(ω)dωE[U^2] = \int_{-\infty}^{\infty} S_{UU}(\omega) d\omega$

响应的标准差：
$σU=E[U2]−(E[U])2\sigma_U = \sqrt{E[U^2] - (E[U])^2}$

4.2 平稳随机响应分析

对于平稳随机激励，系统的响应也是平稳随机过程。

虚拟激励法：

虚拟激励法是一种高效的随机振动分析方法。基本思想是构造一个确定性虚拟激励：
$F~(t)=SFF(ω)eiωt\tilde{F}(t) = \sqrt{S_{FF}(\omega)} e^{i\omega t}$

通过求解确定性方程得到虚拟响应 $U~(t)\tilde{U}(t)$ ，则响应的功率谱密度为：
$SUU(ω)=∣U~(t)∣2S_{UU}(\omega) = |\tilde{U}(t)|^2$

虚拟激励法的优点：

将随机分析转化为确定性分析
适用于大规模结构
可以方便地处理多点相关激励

4.3 非平稳随机响应分析

当激励是非平稳随机过程时，需要采用时域方法或演化谱方法。

时域分析方法：

对于线性系统，响应可以表示为Duhamel积分：
$\int_{0}^{t} h(t-\tau) F(\tau) d\tau$

其中 $h (t)$ 是脉冲响应函数。

响应的均值：
$\int_{0}^{t} h(t-\tau) E[F(\tau)] d\tau$

响应的相关函数：
$RUU(t1,t2)=∫0t1∫0t2h(t1−τ1)h(t2−τ2)RFF(τ1,τ2)dτ1dτ2R_{UU}(t_1, t_2) = \int_{0}^{t_1} \int_{0}^{t_2} h(t_1-\tau_1) h(t_2-\tau_2) R_{FF}(\tau_1, \tau_2) d\tau_1 d\tau_2$

演化谱方法：

对于均匀调制非平稳过程，可以采用演化谱方法：
$SUU(ω,t)=∣A(ω,t)∣2SFF(ω)S_{UU}(\omega, t) = |A(\omega, t)|^2 S_{FF}(\omega)$

其中 $A(ω,t)A(\omega, t)$ 是调制函数。

5. 蒙特卡洛模拟方法

5.1 蒙特卡洛方法的基本原理

蒙特卡洛方法是一种基于随机抽样的统计模拟方法，通过大量随机试验来估计问题的解。

基本步骤：

定义输入随机变量：确定不确定性参数的分布类型和统计特性
生成随机样本：使用随机数生成器产生符合分布的样本
执行确定性分析：对每个样本进行确定性动力分析
统计分析：对结果样本进行统计分析，估计响应的统计特性

样本量的确定：

蒙特卡洛方法的精度与样本量 $N$ 有关。对于估计均值，标准误差为：
$σXˉ=σXN\sigma_{\bar{X}} = \frac{\sigma_X}{\sqrt{N}}$

要达到精度 $ϵ\epsilon$ ，需要的样本量为：
$\left(\frac{\sigma_X z_{\alpha/2}}{\epsilon}\right)^2$

其中 $zα/2z_{\alpha/2}$ 是标准正态分布的分位数。

5.2 随机数生成

伪随机数生成器：

常用的伪随机数生成器包括：

线性同余法：
$X_{n+1} = (a X_n + c) \mod m$

Mersenne Twister：
目前最常用的伪随机数生成器，周期长达 $2^{19937}-1$ 。

逆变换法：

通过均匀分布随机数生成其他分布的随机数：
$X = F^{-1}(U)$

其中 $\sim U(0,1)$ ， $F^{-1}$ 是目标分布的逆累积分布函数。

拒绝采样法：

对于复杂的分布，可以使用拒绝采样法：

从提议分布 $g (x)$ 生成样本 $x$
生成均匀随机数 $\sim U(0,1)$
如果 $\leq \frac{f(x)}{Mg(x)}$ ，接受样本；否则拒绝

5.3 方差缩减技术

为了提高蒙特卡洛方法的效率，可以采用方差缩减技术。

重要性采样：
$\int g(x) f(x) dx = \int g(x) \frac{f(x)}{h(x)} h(x) dx$

其中 $h (x)$ 是重要性采样密度函数，应使 $g (x) f (x) / h (x)$ 的方差最小。

分层采样：
将样本空间划分为若干层，在每层内独立采样：
$μ^=∑i=1LwiXˉi\hat{\mu} = \sum_{i=1}^{L} w_i \bar{X}_i$

其中 $w_i$ 是第 $i$ 层的权重， $Xˉi\bar{X}_i$ 是第 $i$ 层的样本均值。

拉丁超立方采样（LHS）：

LHS是一种高效的空间填充采样方法，可以保证样本在输入空间中均匀分布。

5.4 蒙特卡洛在结构动力学中的应用

应用步骤：

参数不确定性建模：
- 识别不确定性参数
- 确定参数的分布类型和统计特性
- 考虑参数之间的相关性
随机样本生成：
- 使用LHS或其他高效采样方法
- 生成 $N$ 组结构参数样本
确定性动力分析：
- 对每组样本进行动力分析
- 计算感兴趣的响应量
统计分析：
- 计算响应的均值、标准差
- 拟合响应的分布
- 进行可靠性分析

6. 可靠性分析

6.1 结构可靠性基本概念

结构可靠性是指结构在规定的时间内，在规定的条件下，完成预定功能的概率。

功能函数：
$Z = g(X_1, X_2, ..., X_n)$

其中 $X_i$ 是基本随机变量。当 $Z > 0$ 时，结构安全；当 $Z < 0$ 时，结构失效。

失效概率：
$Pf=P(Z<0)=∫g(x)<0fX(x)dxP_f = P(Z < 0) = \int_{g(x)<0} f_X(x) dx$

可靠指标：
$β=Φ−1(1−Pf)\beta = \Phi^{-1}(1 - P_f)$

其中 $Φ\Phi$ 是标准正态分布的累积分布函数。

6.2 一次二阶矩法（FORM）

FORM是一种近似计算失效概率的方法，通过将非线性功能函数在验算点处线性化。

标准正态空间变换：
将相关非正态随机变量变换为独立标准正态变量：
$Ui=Φ−1(FXi(Xi))U_i = \Phi^{-1}(F_{X_i}(X_i))$

功能函数的线性化：
在验算点 $u∗\mathbf{u}^*$ 处进行一阶泰勒展开：
$\approx g(\mathbf{u}^*) + \nabla g(\mathbf{u}^*)^T (\mathbf{u} - \mathbf{u}^*)$

可靠指标的计算：
$β=−g(u∗)+∇g(u∗)Tu∗∣∣∇g(u∗)∣∣\beta = \frac{-g(\mathbf{u}^*) + \nabla g(\mathbf{u}^*)^T \mathbf{u}^*}{||\nabla g(\mathbf{u}^*)||}$

验算点是标准正态空间中功能函数面上距离原点最近的点。

6.3 二次二阶矩法（SORM）

SORM在验算点处进行二次近似，可以提高计算精度。

功能函数的二次近似：
$\approx g(\mathbf{u}^*) + \nabla g(\mathbf{u}^*)^T (\mathbf{u} - \mathbf{u}^*) + \frac{1}{2}(\mathbf{u} - \mathbf{u}^*)^T \mathbf{H} (\mathbf{u} - \mathbf{u}^*)$

其中 $H\mathbf{H}$ 是Hessian矩阵。

失效概率的近似：
$Pf≈Φ(−β)∏i=1n−1(1−βκi)−1/2P_f \approx \Phi(-\beta) \prod_{i=1}^{n-1} (1 - \beta \kappa_i)^{-1/2}$

其中 $κi\kappa_i$ 是主曲率。

6.4 灵敏度分析

灵敏度分析用于评估各随机变量对失效概率的影响程度。

重要性测度：

方向余弦（重要性因子）：
$αi=−∂β∂ui=∂g/∂ui∣∣∇g∣∣\alpha_i = -\frac{\partial \beta}{\partial u_i} = \frac{\partial g/\partial u_i}{||\nabla g||}$

$αi2\alpha_i^2$ 表示第 $i$ 个变量对总不确定性的贡献比例。

弹性指标：
$∂Pf∂μi⋅μiPf\frac{\partial P_f}{\partial \mu_i} \cdot \frac{\mu_i}{P_f}$

表示参数均值变化1%时，失效概率变化的百分比。

基于蒙特卡洛的灵敏度分析：

通过比较不同参数分布下的失效概率变化，可以评估参数的重要性。

7. 案例分析

7.1 案例1：单自由度系统参数不确定性分析

考虑一个单自由度系统，其质量、刚度和阻尼系数存在不确定性。使用摄动方法和蒙特卡洛方法分析系统固有频率和动力响应的统计特性。

"""
案例1：单自由度系统参数不确定性分析
使用摄动方法和蒙特卡洛方法分析参数不确定性对系统动力响应的影响
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from scipy.special import erfcinv
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


def sdof_response(m, c, k, f, dt, u0=0, v0=0):
    """
    单自由度系统动力响应计算（Newmark-beta法）
    
    参数:
        m: 质量
        c: 阻尼系数
        k: 刚度
        f: 外力时程
        dt: 时间步长
        u0, v0: 初始位移和速度
    
    返回:
        u, v, a: 位移、速度、加速度时程
    """
    n_steps = len(f)
    u = np.zeros(n_steps)
    v = np.zeros(n_steps)
    a = np.zeros(n_steps)
    
    u[0] = u0
    v[0] = v0
    a[0] = (f[0] - c * v0 - k * u0) / m
    
    # Newmark-beta参数
    gamma = 0.5
    beta = 0.25
    
    # 等效刚度
    k_eff = k + gamma * c / (beta * dt) + m / (beta * dt**2)
    
    for i in range(n_steps - 1):
        # 等效荷载
        f_eff = (f[i+1] + 
                 m * (u[i] / (beta * dt**2) + v[i] / (beta * dt) + (1/(2*beta)-1) * a[i]) +
                 c * (gamma * u[i] / (beta * dt) + (gamma/beta - 1) * v[i] + 
                      dt * (gamma/(2*beta) - 1) * a[i]))
        
        u[i+1] = f_eff / k_eff
        a[i+1] = (u[i+1] - u[i]) / (beta * dt**2) - v[i] / (beta * dt) - (1/(2*beta)-1) * a[i]
        v[i+1] = v[i] + (1-gamma) * dt * a[i] + gamma * dt * a[i+1]
    
    return u, v, a


def perturbation_analysis(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt):
    """
    一阶摄动法分析参数不确定性
    
    参数:
        m0, c0, k0: 名义质量、阻尼、刚度
        sigma_m, sigma_c, sigma_k: 参数标准差
        f: 外力时程
        dt: 时间步长
    
    返回:
        u_mean: 位移均值
        u_std: 位移标准差
    """
    n_steps = len(f)
    
    # 名义响应
    u0, v0, a0 = sdof_response(m0, c0, k0, f, dt)
    
    # 计算灵敏度（有限差分法）
    delta = 0.01
    
    # 对m的灵敏度
    u_m_plus, _, _ = sdof_response(m0 * (1 + delta), c0, k0, f, dt)
    u_m_minus, _, _ = sdof_response(m0 * (1 - delta), c0, k0, f, dt)
    du_dm = (u_m_plus - u_m_minus) / (2 * delta * m0)
    
    # 对c的灵敏度
    u_c_plus, _, _ = sdof_response(m0, c0 * (1 + delta), k0, f, dt)
    u_c_minus, _, _ = sdof_response(m0, c0 * (1 - delta), k0, f, dt)
    du_dc = (u_c_plus - u_c_minus) / (2 * delta * c0)
    
    # 对k的灵敏度
    u_k_plus, _, _ = sdof_response(m0, c0, k0 * (1 + delta), f, dt)
    u_k_minus, _, _ = sdof_response(m0, c0, k0 * (1 - delta), f, dt)
    du_dk = (u_k_plus - u_k_minus) / (2 * delta * k0)
    
    # 响应方差（假设参数独立）
    u_var = (du_dm * sigma_m)**2 + (du_dc * sigma_c)**2 + (du_dk * sigma_k)**2
    u_std = np.sqrt(u_var)
    
    return u0, u_std, du_dm, du_dc, du_dk


def monte_carlo_analysis(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples=1000):
    """
    蒙特卡洛模拟分析参数不确定性
    
    参数:
        m0, c0, k0: 名义质量、阻尼、刚度
        sigma_m, sigma_c, sigma_k: 参数标准差
        f: 外力时程
        dt: 时间步长
        n_samples: 样本数量
    
    返回:
        u_mean: 位移均值
        u_std: 位移标准差
        u_samples: 所有样本的位移响应
    """
    n_steps = len(f)
    u_samples = np.zeros((n_samples, n_steps))
    
    # 生成随机样本（对数正态分布）
    np.random.seed(42)
    
    # 对数正态分布参数
    zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
    lambda_m = np.log(m0) - 0.5 * zeta_m**2
    
    zeta_c = np.sqrt(np.log(1 + (sigma_c/c0)**2))
    lambda_c = np.log(c0) - 0.5 * zeta_c**2
    
    zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
    lambda_k = np.log(k0) - 0.5 * zeta_k**2
    
    for i in range(n_samples):
        m_sample = np.random.lognormal(lambda_m, zeta_m)
        c_sample = np.random.lognormal(lambda_c, zeta_c)
        k_sample = np.random.lognormal(lambda_k, zeta_k)
        
        u, _, _ = sdof_response(m_sample, c_sample, k_sample, f, dt)
        u_samples[i, :] = u
    
    u_mean = np.mean(u_samples, axis=0)
    u_std = np.std(u_samples, axis=0)
    
    return u_mean, u_std, u_samples


def latin_hypercube_sampling(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples=1000):
    """
    拉丁超立方采样（LHS）分析参数不确定性
    
    参数:
        m0, c0, k0: 名义质量、阻尼、刚度
        sigma_m, sigma_c, sigma_k: 参数标准差
        f: 外力时程
        dt: 时间步长
        n_samples: 样本数量
    
    返回:
        u_mean: 位移均值
        u_std: 位移标准差
        u_samples: 所有样本的位移响应
    """
    n_steps = len(f)
    u_samples = np.zeros((n_samples, n_steps))
    
    # LHS采样
    np.random.seed(42)
    
    # 生成拉丁超立方样本（标准正态空间）
    def lhs_normal(n_samples, n_vars):
        """生成LHS标准正态样本"""
        samples = np.zeros((n_samples, n_vars))
        for i in range(n_vars):
            # 分层采样
            strata = np.random.permutation(n_samples)
            for j in range(n_samples):
                # 在每一层内均匀采样，然后转换为正态分布
                u = (strata[j] + np.random.rand()) / n_samples
                samples[j, i] = np.sqrt(2) * erfcinv(2 * (1 - u))
        return samples
    
    lhs_samples = lhs_normal(n_samples, 3)
    
    # 转换为对数正态分布
    zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
    lambda_m = np.log(m0) - 0.5 * zeta_m**2
    
    zeta_c = np.sqrt(np.log(1 + (sigma_c/c0)**2))
    lambda_c = np.log(c0) - 0.5 * zeta_c**2
    
    zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
    lambda_k = np.log(k0) - 0.5 * zeta_k**2
    
    for i in range(n_samples):
        m_sample = np.exp(lambda_m + zeta_m * lhs_samples[i, 0])
        c_sample = np.exp(lambda_c + zeta_c * lhs_samples[i, 1])
        k_sample = np.exp(lambda_k + zeta_k * lhs_samples[i, 2])
        
        u, _, _ = sdof_response(m_sample, c_sample, k_sample, f, dt)
        u_samples[i, :] = u
    
    u_mean = np.mean(u_samples, axis=0)
    u_std = np.std(u_samples, axis=0)
    
    return u_mean, u_std, u_samples


# ==================== 主程序 ====================
print("=" * 60)
print("案例1：单自由度系统参数不确定性分析")
print("=" * 60)

# 系统参数（名义值）
m0 = 1000.0      # 质量 (kg)
c0 = 2000.0      # 阻尼系数 (N·s/m)
k0 = 100000.0    # 刚度 (N/m)

# 不确定性水平（变异系数）
cov_m = 0.1      # 质量变异系数
cov_c = 0.15     # 阻尼变异系数
cov_k = 0.1      # 刚度变异系数

sigma_m = m0 * cov_m
sigma_c = c0 * cov_c
sigma_k = k0 * cov_k

print(f"\n系统名义参数:")
print(f"  质量 m0 = {m0:.2f} kg")
print(f"  阻尼 c0 = {c0:.2f} N·s/m")
print(f"  刚度 k0 = {k0:.2f} N/m")

print(f"\n不确定性水平（变异系数）:")
print(f"  质量 COV = {cov_m*100:.1f}%")
print(f"  阻尼 COV = {cov_c*100:.1f}%")
print(f"  刚度 COV = {cov_k*100:.1f}%")

# 计算名义系统的固有特性
omega_n = np.sqrt(k0 / m0)  # 固有频率
xi = c0 / (2 * np.sqrt(k0 * m0))  # 阻尼比
T_n = 2 * np.pi / omega_n  # 固有周期

print(f"\n名义系统动力特性:")
print(f"  固有频率 fn = {omega_n/(2*np.pi):.3f} Hz")
print(f"  阻尼比 xi = {xi:.4f}")
print(f"  固有周期 Tn = {T_n:.4f} s")

# 时间参数
dt = 0.01
t_max = 10.0
t = np.arange(0, t_max, dt)
n_steps = len(t)

# 定义外力（脉冲荷载）
f = np.zeros(n_steps)
pulse_start = int(0.5 / dt)
pulse_end = int(1.0 / dt)
f[pulse_start:pulse_end] = 10000.0  # 脉冲幅值

print(f"\n荷载条件:")
print(f"  脉冲荷载幅值 = 10000 N")
print(f"  脉冲持续时间 = 0.5 s")

# ==================== 摄动法分析 ====================
print("\n" + "=" * 60)
print("一阶摄动法分析")
print("=" * 60)

u_pert, u_pert_std, du_dm, du_dc, du_dk = perturbation_analysis(
    m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt
)

print(f"\n摄动法结果:")
print(f"  最大位移均值 = {np.max(np.abs(u_pert)):.6f} m")
print(f"  最大位移标准差 = {np.max(u_pert_std):.6f} m")
print(f"  位移变异系数 = {np.max(u_pert_std)/np.max(np.abs(u_pert))*100:.2f}%")

# ==================== 蒙特卡洛分析 ====================
print("\n" + "=" * 60)
print("蒙特卡洛模拟分析")
print("=" * 60)

n_samples_mc = 2000
print(f"\n样本数量: {n_samples_mc}")

u_mc_mean, u_mc_std, u_mc_samples = monte_carlo_analysis(
    m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples_mc
)

print(f"\n蒙特卡洛结果:")
print(f"  最大位移均值 = {np.max(np.abs(u_mc_mean)):.6f} m")
print(f"  最大位移标准差 = {np.max(u_mc_std):.6f} m")
print(f"  位移变异系数 = {np.max(u_mc_std)/np.max(np.abs(u_mc_mean))*100:.2f}%")

# ==================== LHS分析 ====================
print("\n" + "=" * 60)
print("拉丁超立方采样（LHS）分析")
print("=" * 60)

u_lhs_mean, u_lhs_std, u_lhs_samples = latin_hypercube_sampling(
    m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples_mc
)

print(f"\nLHS结果:")
print(f"  最大位移均值 = {np.max(np.abs(u_lhs_mean)):.6f} m")
print(f"  最大位移标准差 = {np.max(u_lhs_std):.6f} m")
print(f"  位移变异系数 = {np.max(u_lhs_std)/np.max(np.abs(u_lhs_mean))*100:.2f}%")

# ==================== 方法对比 ====================
print("\n" + "=" * 60)
print("方法对比")
print("=" * 60)

print(f"\n最大位移统计量对比:")
print(f"{'方法':<20} {'均值(m)':<15} {'标准差(m)':<15} {'变异系数(%)':<15}")
print("-" * 65)
print(f"{'摄动法':<20} {np.max(np.abs(u_pert)):<15.6f} {np.max(u_pert_std):<15.6f} {np.max(u_pert_std)/np.max(np.abs(u_pert))*100:<15.2f}")
print(f"{'蒙特卡洛':<20} {np.max(np.abs(u_mc_mean)):<15.6f} {np.max(u_mc_std):<15.6f} {np.max(u_mc_std)/np.max(np.abs(u_mc_mean))*100:<15.2f}")
print(f"{'LHS':<20} {np.max(np.abs(u_lhs_mean)):<15.6f} {np.max(u_lhs_std):<15.6f} {np.max(u_lhs_std)/np.max(np.abs(u_lhs_mean))*100:<15.2f}")

# ==================== 绘制结果 ====================
print("\n正在生成可视化结果...")

# 图1：位移响应均值和标准差
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：位移均值对比
ax = axes[0, 0]
ax.plot(t, u_pert, 'b-', linewidth=2, label='Perturbation Method')
ax.plot(t, u_mc_mean, 'r--', linewidth=2, label='Monte Carlo')
ax.plot(t, u_lhs_mean, 'g:', linewidth=2, label='LHS')
ax.set_xlabel('Time (s)', fontsize=11)
ax.set_ylabel('Displacement (m)', fontsize=11)
ax.set_title('Mean Displacement Response', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图2：位移标准差对比
ax = axes[0, 1]
ax.plot(t, u_pert_std, 'b-', linewidth=2, label='Perturbation Method')
ax.plot(t, u_mc_std, 'r--', linewidth=2, label='Monte Carlo')
ax.plot(t, u_lhs_std, 'g:', linewidth=2, label='LHS')
ax.set_xlabel('Time (s)', fontsize=11)
ax.set_ylabel('Standard Deviation (m)', fontsize=11)
ax.set_title('Displacement Standard Deviation', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图3：蒙特卡洛样本云图
ax = axes[1, 0]
for i in range(min(100, n_samples_mc)):
    ax.plot(t, u_mc_samples[i, :], 'gray', alpha=0.1, linewidth=0.5)
ax.plot(t, u_mc_mean, 'r-', linewidth=2, label='Mean')
ax.plot(t, u_mc_mean + 2*u_mc_std, 'r--', linewidth=1.5, label='Mean ± 2σ')
ax.plot(t, u_mc_mean - 2*u_mc_std, 'r--', linewidth=1.5)
ax.set_xlabel('Time (s)', fontsize=11)
ax.set_ylabel('Displacement (m)', fontsize=11)
ax.set_title('Monte Carlo Response Samples (100 samples)', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图4：响应概率密度（特定时刻）
ax = axes[1, 1]
time_idx = np.argmax(np.abs(u_mc_mean))
response_at_peak = u_mc_samples[:, time_idx]
ax.hist(response_at_peak, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
ax.axvline(u_mc_mean[time_idx], color='r', linestyle='--', linewidth=2, label=f'Mean = {u_mc_mean[time_idx]:.4f} m')
ax.axvline(u_mc_mean[time_idx] + 2*u_mc_std[time_idx], color='g', linestyle=':', linewidth=2, label=f'Mean + 2σ = {u_mc_mean[time_idx] + 2*u_mc_std[time_idx]:.4f} m')
ax.axvline(u_mc_mean[time_idx] - 2*u_mc_std[time_idx], color='g', linestyle=':', linewidth=2)
ax.set_xlabel('Displacement (m)', fontsize=11)
ax.set_ylabel('Probability Density', fontsize=11)
ax.set_title(f'Response PDF at t = {t[time_idx]:.2f} s (Peak Response)', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case1_parameter_uncertainty_analysis.png', dpi=150, bbox_inches='tight')
plt.close()

print("  已保存: case1_parameter_uncertainty_analysis.png")

# 图2：参数灵敏度分析
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：位移对质量的灵敏度
ax = axes[0, 0]
ax.plot(t, du_dm, 'b-', linewidth=2)
ax.set_xlabel('Time (s)', fontsize=11)
ax.set_ylabel('du/dm (m/kg)', fontsize=11)
ax.set_title('Sensitivity of Displacement to Mass', fontsize=12)
ax.grid(True, alpha=0.3)

# 子图2：位移对阻尼的灵敏度
ax = axes[0, 1]
ax.plot(t, du_dc, 'r-', linewidth=2)
ax.set_xlabel('Time (s)', fontsize=11)
ax.set_ylabel('du/dc (m/(N·s/m))', fontsize=11)
ax.set_title('Sensitivity of Displacement to Damping', fontsize=12)
ax.grid(True, alpha=0.3)

# 子图3：位移对刚度的灵敏度
ax = axes[1, 0]
ax.plot(t, du_dk, 'g-', linewidth=2)
ax.set_xlabel('Time (s)', fontsize=11)
ax.set_ylabel('du/dk (m/(N/m))', fontsize=11)
ax.set_title('Sensitivity of Displacement to Stiffness', fontsize=12)
ax.grid(True, alpha=0.3)

# 子图4：参数不确定性贡献
ax = axes[1, 1]
# 计算各参数对响应方差的贡献
var_m = (du_dm * sigma_m)**2
var_c = (du_dc * sigma_c)**2
var_k = (du_dk * sigma_k)**2
var_total = var_m + var_c + var_k

# 绘制贡献比例随时间的变化
ax.plot(t, var_m/var_total * 100, 'b-', linewidth=2, label='Mass')
ax.plot(t, var_c/var_total * 100, 'r-', linewidth=2, label='Damping')
ax.plot(t, var_k/var_total * 100, 'g-', linewidth=2, label='Stiffness')
ax.set_xlabel('Time (s)', fontsize=11)
ax.set_ylabel('Contribution to Variance (%)', fontsize=11)
ax.set_title('Parameter Uncertainty Contribution', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case1_sensitivity_analysis.png', dpi=150, bbox_inches='tight')
plt.close()

print("  已保存: case1_sensitivity_analysis.png")

# 图3：固有频率的不确定性
print("\n正在分析固有频率的不确定性...")

# 生成大量样本计算固有频率
n_freq_samples = 5000
np.random.seed(123)

zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
lambda_m = np.log(m0) - 0.5 * zeta_m**2
zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
lambda_k = np.log(k0) - 0.5 * zeta_k**2

m_samples = np.random.lognormal(lambda_m, zeta_m, n_freq_samples)
k_samples = np.random.lognormal(lambda_k, zeta_k, n_freq_samples)
freq_samples = np.sqrt(k_samples / m_samples) / (2 * np.pi)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 子图1：固有频率分布
ax = axes[0]
ax.hist(freq_samples, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
ax.axvline(omega_n/(2*np.pi), color='r', linestyle='--', linewidth=2, label=f'Nominal = {omega_n/(2*np.pi):.3f} Hz')
ax.axvline(np.mean(freq_samples), color='g', linestyle=':', linewidth=2, label=f'Mean = {np.mean(freq_samples):.3f} Hz')
ax.set_xlabel('Natural Frequency (Hz)', fontsize=11)
ax.set_ylabel('Probability Density', fontsize=11)
ax.set_title('Uncertainty in Natural Frequency', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图2：固有频率统计量
ax = axes[1]
stats_text = f"Natural Frequency Statistics:\n"
stats_text += f"  Nominal: {omega_n/(2*np.pi):.4f} Hz\n"
stats_text += f"  Mean: {np.mean(freq_samples):.4f} Hz\n"
stats_text += f"  Std Dev: {np.std(freq_samples):.4f} Hz\n"
stats_text += f"  COV: {np.std(freq_samples)/np.mean(freq_samples)*100:.2f}%\n\n"
stats_text += f"  5th percentile: {np.percentile(freq_samples, 5):.4f} Hz\n"
stats_text += f"  95th percentile: {np.percentile(freq_samples, 95):.4f} Hz\n\n"
stats_text += f"Damping Ratio Statistics:\n"
c_samples = c0 * np.random.lognormal(np.log(1)-0.5*0.15**2, 0.15, n_freq_samples)
xi_samples = c_samples / (2 * np.sqrt(k_samples * m_samples))
stats_text += f"  Nominal: {xi:.4f}\n"
stats_text += f"  Mean: {np.mean(xi_samples):.4f}\n"
stats_text += f"  Std Dev: {np.std(xi_samples):.4f}\n"
stats_text += f"  COV: {np.std(xi_samples)/np.mean(xi_samples)*100:.2f}%"

ax.text(0.1, 0.5, stats_text, fontsize=11, verticalalignment='center',
        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
ax.axis('off')

plt.tight_layout()
plt.savefig('case1_natural_frequency_uncertainty.png', dpi=150, bbox_inches='tight')
plt.close()

print("  已保存: case1_natural_frequency_uncertainty.png")

# 创建动画：响应样本的时间演化
print("\n正在生成响应样本动画...")

fig, ax = plt.subplots(figsize=(10, 6))
ax.set_xlim(0, t_max)
ax.set_ylim(np.min(u_mc_samples) * 1.1, np.max(u_mc_samples) * 1.1)
ax.set_xlabel('Time (s)', fontsize=12)
ax.set_ylabel('Displacement (m)', fontsize=12)
ax.set_title('SDOF Response with Parameter Uncertainty', fontsize=14)
ax.grid(True, alpha=0.3)

# 绘制均值和置信区间
line_mean, = ax.plot([], [], 'r-', linewidth=2.5, label='Mean')
line_upper, = ax.plot([], [], 'r--', linewidth=1.5, alpha=0.7, label='Mean ± 2σ')
line_lower, = ax.plot([], [], 'r--', linewidth=1.5, alpha=0.7)

# 绘制样本
sample_lines = []
for i in range(20):
    line, = ax.plot([], [], 'gray', alpha=0.3, linewidth=0.8)
    sample_lines.append(line)

ax.legend(loc='upper right')

# 外力曲线
ax_force = ax.twinx()
line_force, = ax_force.plot([], [], 'b:', linewidth=1.5, alpha=0.5, label='Force')
ax_force.set_ylabel('Force (N)', fontsize=12, color='blue')
ax_force.tick_params(axis='y', labelcolor='blue')

def init():
    line_mean.set_data([], [])
    line_upper.set_data([], [])
    line_lower.set_data([], [])
    line_force.set_data([], [])
    for line in sample_lines:
        line.set_data([], [])
    return [line_mean, line_upper, line_lower, line_force] + sample_lines

def update(frame):
    end_idx = frame + 1
    line_mean.set_data(t[:end_idx], u_mc_mean[:end_idx])
    line_upper.set_data(t[:end_idx], (u_mc_mean + 2*u_mc_std)[:end_idx])
    line_lower.set_data(t[:end_idx], (u_mc_mean - 2*u_mc_std)[:end_idx])
    line_force.set_data(t[:end_idx], f[:end_idx] / 10000 * np.max(u_mc_samples) * 0.5)
    
    for i, line in enumerate(sample_lines):
        line.set_data(t[:end_idx], u_mc_samples[i, :end_idx])
    
    return [line_mean, line_upper, line_lower, line_force] + sample_lines

# 每5帧保存一次，减少动画大小
frame_indices = np.arange(0, n_steps, 5)
anim = FuncAnimation(fig, update, frames=frame_indices, init_func=init, blit=True)

writer = PillowWriter(fps=20)
anim.save('case1_response_animation.gif', writer=writer)
plt.close()

print("  已保存: case1_response_animation.gif")

print("\n" + "=" * 60)
print("案例1分析完成!")
print("=" * 60)

7.2 案例2：多自由度系统随机振动响应分析

考虑一个多自由度结构在随机地震激励下的响应。使用虚拟激励法和时域分析方法计算响应的功率谱密度和统计量。

"""
案例2：多自由度系统随机振动响应分析
使用虚拟激励法和时域分析方法计算随机地震激励下的响应统计量
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from scipy import linalg
from scipy.integrate import trapezoid
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


def build_mdof_system(n_dof, m_elem, k_elem, damping_ratio=0.05):
    """
    建立多自由度系统（剪切型结构）
    
    参数:
        n_dof: 自由度数
        m_elem: 单元质量
        k_elem: 单元刚度
        damping_ratio: 阻尼比
    
    返回:
        M, C, K: 质量、阻尼、刚度矩阵
        Phi: 模态矩阵
        omega: 固有频率
    """
    # 质量矩阵（集中质量）
    M = np.diag([m_elem] * n_dof)
    
    # 刚度矩阵（三对角）
    K = np.zeros((n_dof, n_dof))
    for i in range(n_dof):
        K[i, i] = 2 * k_elem if i < n_dof - 1 else k_elem
        if i > 0:
            K[i, i] += k_elem
            K[i, i-1] = -k_elem
            K[i-1, i] = -k_elem
    
    # 模态分析
    eigenvalues, eigenvectors = linalg.eigh(K, M)
    omega = np.sqrt(eigenvalues)
    Phi = eigenvectors
    
    # 瑞利阻尼
    omega1, omega2 = omega[0], omega[1]
    a0 = 2 * damping_ratio * omega1 * omega2 / (omega1 + omega2)
    a1 = 2 * damping_ratio / (omega1 + omega2)
    C = a0 * M + a1 * K
    
    return M, C, K, Phi, omega


def kanai_tajimi_psd(omega, omega_g, xi_g, S0):
    """
    Kanai-Tajimi地震动功率谱密度模型
    
    参数:
        omega: 频率
        omega_g: 场地特征频率
        xi_g: 场地阻尼比
        S0: 谱强度因子
    
    返回:
        S: 功率谱密度
    """
    S = S0 * (1 + (2 * xi_g * omega / omega_g)**2) / \
        ((1 - (omega / omega_g)**2)**2 + (2 * xi_g * omega / omega_g)**2)
    return S


def clough_penzien_psd(omega, omega_g, xi_g, omega_f, xi_f, S0):
    """
    Clough-Penzien地震动功率谱密度模型
    
    参数:
        omega: 频率
        omega_g: 场地特征频率
        xi_g: 场地阻尼比
        omega_f: 低频滤波频率
        xi_f: 低频滤波阻尼比
        S0: 谱强度因子
    
    返回:
        S: 功率谱密度
    """
    S_kt = kanai_tajimi_psd(omega, omega_g, xi_g, S0)
    H_f = (omega / omega_f)**4 / ((1 - (omega / omega_f)**2)**2 + (2 * xi_f * omega / omega_f)**2)
    S = S_kt * H_f
    return S


def pseudo_excitation_method(M, C, K, Phi, omega_n, S_ff, omega_list):
    """
    虚拟激励法计算随机振动响应
    
    参数:
        M, C, K: 质量、阻尼、刚度矩阵
        Phi: 模态矩阵
        omega_n: 固有频率
        S_ff: 激励功率谱密度函数
        omega_list: 频率列表
    
    返回:
        S_uu: 位移功率谱密度
        S_vv: 速度功率谱密度
        S_aa: 加速度功率谱密度
    """
    n_dof = M.shape[0]
    n_freq = len(omega_list)
    
    S_uu = np.zeros((n_dof, n_freq))
    S_vv = np.zeros((n_dof, n_freq))
    S_aa = np.zeros((n_dof, n_freq))
    
    # 影响向量（地震激励）
    r = np.ones(n_dof)
    
    for i, omega in enumerate(omega_list):
        # 虚拟激励幅值
        sqrt_Sff = np.sqrt(S_ff(omega))
        
        # 频率响应函数
        H = np.linalg.inv(-omega**2 * M + 1j * omega * C + K)
        
        # 虚拟响应
        u_tilde = H @ r * sqrt_Sff
        v_tilde = 1j * omega * u_tilde
        a_tilde = -omega**2 * u_tilde
        
        # 功率谱密度
        S_uu[:, i] = np.abs(u_tilde)**2
        S_vv[:, i] = np.abs(v_tilde)**2
        S_aa[:, i] = np.abs(a_tilde)**2
    
    return S_uu, S_vv, S_aa


def modal_superposition_psd(Phi, omega_n, damping_ratio, S_ff, omega_list):
    """
    模态叠加法计算功率谱密度
    
    参数:
        Phi: 模态矩阵
        omega_n: 固有频率
        damping_ratio: 模态阻尼比
        S_ff: 激励功率谱密度函数
        omega_list: 频率列表
    
    返回:
        S_uu_modal: 各模态的位移功率谱密度贡献
    """
    n_dof = Phi.shape[0]
    n_modes = len(omega_n)
    n_freq = len(omega_list)
    
    S_uu_modal = np.zeros((n_modes, n_freq))
    
    for i in range(n_modes):
        omega_i = omega_n[i]
        xi_i = damping_ratio
        
        for j, omega in enumerate(omega_list):
            # 模态频率响应函数
            H_i = 1 / (omega_i**2 - omega**2 + 2j * xi_i * omega_i * omega)
            
            # 模态功率谱密度
            S_uu_modal[i, j] = np.abs(H_i)**2 * S_ff(omega)
    
    return S_uu_modal


def generate_earthquake_acceleration(t, dt, omega_g, xi_g, S0, seed=None):
    """
    生成人工地震动加速度时程（基于Kanai-Tajimi模型）
    
    参数:
        t: 时间向量
        dt: 时间步长
        omega_g: 场地特征频率
        xi_g: 场地阻尼比
        S0: 谱强度因子
        seed: 随机种子
    
    返回:
        acc: 加速度时程
    """
    if seed is not None:
        np.random.seed(seed)
    
    n_steps = len(t)
    
    # 白噪声
    w = np.random.randn(n_steps)
    
    # Kanai-Tajimi滤波器（时域实现）
    # 滤波器参数
    omega_d = omega_g * np.sqrt(1 - xi_g**2)
    
    # 离散化滤波器
    acc = np.zeros(n_steps)
    u = np.zeros(n_steps)  # 滤波器状态
    v = np.zeros(n_steps)  # 滤波器速度
    
    # 滤波器系数
    a1 = 2 * np.exp(-xi_g * omega_g * dt) * np.cos(omega_d * dt)
    a2 = -np.exp(-2 * xi_g * omega_g * dt)
    b0 = 2 * xi_g * omega_g * dt
    b1 = omega_g**2 * dt
    
    # 简化实现：使用状态空间方法
    for i in range(1, n_steps):
        # 滤波器方程
        du = v[i-1]
        dv = -omega_g**2 * u[i-1] - 2 * xi_g * omega_g * v[i-1] + w[i]
        
        u[i] = u[i-1] + du * dt
        v[i] = v[i-1] + dv * dt
        
        # 输出加速度
        acc[i] = -omega_g**2 * u[i] - 2 * xi_g * omega_g * v[i] + w[i]
    
    # 归一化
    acc = acc * np.sqrt(S0 * np.pi / (2 * xi_g * omega_g**3))
    
    return acc


def newmark_beta_mdof(M, C, K, f, dt, u0=None, v0=None):
    """
    Newmark-beta法求解多自由度系统动力响应
    
    参数:
        M, C, K: 质量、阻尼、刚度矩阵
        f: 外力矩阵 (n_steps x n_dof)
        dt: 时间步长
        u0, v0: 初始位移和速度
    
    返回:
        u, v, a: 位移、速度、加速度时程
    """
    n_steps, n_dof = f.shape
    
    u = np.zeros((n_steps, n_dof))
    v = np.zeros((n_steps, n_dof))
    a = np.zeros((n_steps, n_dof))
    
    if u0 is not None:
        u[0] = u0
    if v0 is not None:
        v[0] = v0
    
    # 初始加速度
    a[0] = np.linalg.solve(M, f[0] - C @ v[0] - K @ u[0])
    
    # Newmark参数
    gamma = 0.5
    beta = 0.25
    
    # 等效刚度
    K_eff = K + gamma * C / (beta * dt) + M / (beta * dt**2)
    
    for i in range(n_steps - 1):
        # 等效荷载
        f_eff = (f[i+1] + 
                 M @ (u[i] / (beta * dt**2) + v[i] / (beta * dt) + (1/(2*beta)-1) * a[i]) +
                 C @ (gamma * u[i] / (beta * dt) + (gamma/beta - 1) * v[i] + 
                      dt * (gamma/(2*beta) - 1) * a[i]))
        
        u[i+1] = np.linalg.solve(K_eff, f_eff)
        a[i+1] = (u[i+1] - u[i]) / (beta * dt**2) - v[i] / (beta * dt) - (1/(2*beta)-1) * a[i]
        v[i+1] = v[i] + (1-gamma) * dt * a[i] + gamma * dt * a[i+1]
    
    return u, v, a


# ==================== 主程序 ====================
print("=" * 60)
print("案例2：多自由度系统随机振动响应分析")
print("=" * 60)

# 结构参数
n_dof = 5
m_elem = 1000.0  # 每层质量 (kg)
k_elem = 50000.0  # 每层刚度 (N/m)
damping_ratio = 0.05  # 阻尼比

print(f"\n结构参数:")
print(f"  层数: {n_dof}")
print(f"  每层质量: {m_elem} kg")
print(f"  每层刚度: {k_elem} N/m")
print(f"  阻尼比: {damping_ratio}")

# 建立结构模型
M, C, K, Phi, omega = build_mdof_system(n_dof, m_elem, k_elem, damping_ratio)

print(f"\n结构动力特性:")
for i in range(min(3, n_dof)):
    print(f"  第{i+1}阶频率: {omega[i]/(2*np.pi):.3f} Hz, 周期: {2*np.pi/omega[i]:.3f} s")

# 地震动参数
omega_g = 2 * np.pi * 2.0  # 场地特征频率 (2 Hz)
xi_g = 0.6  # 场地阻尼比
S0 = 0.01  # 谱强度因子

print(f"\n地震动参数 (Kanai-Tajimi模型):")
print(f"  场地特征频率: {omega_g/(2*np.pi):.1f} Hz")
print(f"  场地阻尼比: {xi_g}")
print(f"  谱强度因子: {S0}")

# 频率范围
omega_min = 0.1
omega_max = 50.0
n_freq = 500
omega_list = np.linspace(omega_min, omega_max, n_freq)

# 定义功率谱密度函数
def S_ff_func(omega):
    return kanai_tajimi_psd(omega, omega_g, xi_g, S0)

# ==================== 虚拟激励法分析 ====================
print("\n" + "=" * 60)
print("虚拟激励法分析")
print("=" * 60)

S_uu, S_vv, S_aa = pseudo_excitation_method(M, C, K, Phi, omega, S_ff_func, omega_list)

# 计算响应的均方值
sigma_u = np.sqrt(trapezoid(S_uu, omega_list, axis=1))
sigma_v = np.sqrt(trapezoid(S_vv, omega_list, axis=1))
sigma_a = np.sqrt(trapezoid(S_aa, omega_list, axis=1))

print(f"\n各层位移响应标准差:")
for i in range(n_dof):
    print(f"  第{i+1}层: {sigma_u[i]*1000:.4f} mm")

print(f"\n各层加速度响应标准差:")
for i in range(n_dof):
    print(f"  第{i+1}层: {sigma_a[i]:.4f} m/s² ({sigma_a[i]/9.81*100:.2f}%g)")

# ==================== 模态叠加法分析 ====================
print("\n" + "=" * 60)
print("模态叠加法分析")
print("=" * 60)

S_uu_modal = modal_superposition_psd(Phi, omega, damping_ratio, S_ff_func, omega_list)

print(f"\n各模态的位移响应贡献:")
for i in range(min(3, n_dof)):
    sigma_modal = np.sqrt(trapezoid(S_uu_modal[i, :], omega_list))
    print(f"  第{i+1}阶模态: {sigma_modal*1000:.4f} mm")

# ==================== 时域蒙特卡洛验证 ====================
print("\n" + "=" * 60)
print("时域蒙特卡洛验证")
print("=" * 60)

# 时间参数
dt = 0.01
t_max = 30.0
t = np.arange(0, t_max, dt)
n_steps = len(t)

# 生成多个地震动样本进行蒙特卡洛模拟
n_samples = 100
u_samples = np.zeros((n_samples, n_steps, n_dof))

print(f"\n生成 {n_samples} 条地震动样本进行时域分析...")

for i in range(n_samples):
    # 生成地震动加速度
    acc = generate_earthquake_acceleration(t, dt, omega_g, xi_g, S0, seed=i*100)
    
    # 地震力（假设质量集中在各层）
    f = np.outer(acc, M @ np.ones(n_dof))
    
    # 计算响应
    u, v, a = newmark_beta_mdof(M, C, K, f, dt)
    u_samples[i, :, :] = u

# 统计响应
u_mean_mc = np.mean(u_samples, axis=0)
u_std_mc = np.std(u_samples, axis=0)

print(f"\n时域蒙特卡洛结果 (各层位移标准差):")
for i in range(n_dof):
    max_std = np.max(u_std_mc[:, i])
    print(f"  第{i+1}层: {max_std*1000:.4f} mm")

# ==================== 绘制结果 ====================
print("\n正在生成可视化结果...")

# 图1：功率谱密度分析
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：地震动功率谱密度
ax = axes[0, 0]
S_ff_values = [S_ff_func(om) for om in omega_list]
ax.semilogy(omega_list / (2*np.pi), S_ff_values, 'b-', linewidth=2)
ax.axvline(omega_g / (2*np.pi), color='r', linestyle='--', label=f'Site freq = {omega_g/(2*np.pi):.1f} Hz')
ax.set_xlabel('Frequency (Hz)', fontsize=11)
ax.set_ylabel('PSD ((m/s²)²/Hz)', fontsize=11)
ax.set_title('Kanai-Tajimi PSD Model', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 10)

# 子图2：各层位移功率谱密度
ax = axes[0, 1]
colors = plt.cm.viridis(np.linspace(0, 1, n_dof))
for i in range(n_dof):
    ax.semilogy(omega_list / (2*np.pi), S_uu[i, :], color=colors[i], linewidth=1.5, label=f'Floor {i+1}')
ax.set_xlabel('Frequency (Hz)', fontsize=11)
ax.set_ylabel('Displacement PSD (m²/Hz)', fontsize=11)
ax.set_title('Displacement PSD (Pseudo-Excitation)', fontsize=12)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 10)

# 子图3：各层加速度功率谱密度
ax = axes[1, 0]
for i in range(n_dof):
    ax.semilogy(omega_list / (2*np.pi), S_aa[i, :], color=colors[i], linewidth=1.5, label=f'Floor {i+1}')
ax.set_xlabel('Frequency (Hz)', fontsize=11)
ax.set_ylabel('Acceleration PSD ((m/s²)²/Hz)', fontsize=11)
ax.set_title('Acceleration PSD (Pseudo-Excitation)', fontsize=12)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 10)

# 子图4：模态贡献
ax = axes[1, 1]
for i in range(min(3, n_dof)):
    ax.semilogy(omega_list / (2*np.pi), S_uu_modal[i, :], linewidth=1.5, label=f'Mode {i+1}')
ax.set_xlabel('Frequency (Hz)', fontsize=11)
ax.set_ylabel('Modal PSD (m²/Hz)', fontsize=11)
ax.set_title('Modal Contribution to Displacement PSD', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 10)

plt.tight_layout()
plt.savefig('case2_psd_analysis.png', dpi=150, bbox_inches='tight')
plt.close()

print("  已保存: case2_psd_analysis.png")

# 图2：响应统计量对比
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：各层位移标准差对比
ax = axes[0, 0]
floors = np.arange(1, n_dof + 1)
ax.bar(floors - 0.2, sigma_u * 1000, 0.4, label='Pseudo-Excitation', color='blue', alpha=0.7)
ax.bar(floors + 0.2, np.max(u_std_mc, axis=0) * 1000, 0.4, label='Time Domain MC', color='red', alpha=0.7)
ax.set_xlabel('Floor', fontsize=11)
ax.set_ylabel('Displacement Std Dev (mm)', fontsize=11)
ax.set_title('Displacement Response Standard Deviation', fontsize=12)
ax.set_xticks(floors)
ax.legend()
ax.grid(True, alpha=0.3, axis='y')

# 子图2：时域响应样本（顶层）
ax = axes[0, 1]
for i in range(min(10, n_samples)):
    ax.plot(t, u_samples[i, :, -1] * 1000, 'gray', alpha=0.3, linewidth=0.5)
ax.plot(t, u_mean_mc[:, -1] * 1000, 'r-', linewidth=2, label='Mean')
ax.plot(t, (u_mean_mc[:, -1] + 2*u_std_mc[:, -1]) * 1000, 'r--', linewidth=1.5, label='Mean ± 2σ')
ax.plot(t, (u_mean_mc[:, -1] - 2*u_std_mc[:, -1]) * 1000, 'r--', linewidth=1.5)
ax.set_xlabel('Time (s)', fontsize=11)
ax.set_ylabel('Displacement (mm)', fontsize=11)
ax.set_title('Top Floor Response Samples', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, t_max)

# 子图3：响应概率密度（顶层，t=10s）
ax = axes[1, 0]
time_idx = int(10 / dt)
response_at_t = u_samples[:, time_idx, -1] * 1000
ax.hist(response_at_t, bins=30, density=True, alpha=0.7, color='blue', edgecolor='black')
ax.axvline(np.mean(response_at_t), color='r', linestyle='--', linewidth=2, label=f'Mean = {np.mean(response_at_t):.2f} mm')
ax.axvline(np.mean(response_at_t) + 2*np.std(response_at_t), color='g', linestyle=':', linewidth=2)
ax.axvline(np.mean(response_at_t) - 2*np.std(response_at_t), color='g', linestyle=':', linewidth=2, label=f'±2σ = {2*np.std(response_at_t):.2f} mm')
ax.set_xlabel('Displacement (mm)', fontsize=11)
ax.set_ylabel('Probability Density', fontsize=11)
ax.set_title(f'Response PDF at t = 10s (Top Floor)', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图4：频域和时域方法对比
ax = axes[1, 1]
methods = ['Pseudo-Excitation', 'Time Domain MC']
displacement_top = [sigma_u[-1] * 1000, np.max(u_std_mc[:, -1]) * 1000]
ax.bar(methods, displacement_top, color=['blue', 'red'], alpha=0.7)
ax.set_ylabel('Top Floor Displacement Std Dev (mm)', fontsize=11)
ax.set_title('Method Comparison', fontsize=12)
ax.grid(True, alpha=0.3, axis='y')
for i, v in enumerate(displacement_top):
    ax.text(i, v + 0.5, f'{v:.2f}', ha='center', fontsize=11)

plt.tight_layout()
plt.savefig('case2_response_statistics.png', dpi=150, bbox_inches='tight')
plt.close()

print("  已保存: case2_response_statistics.png")

# 图3：人工地震动样本
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 生成3条地震动样本
for idx in range(4):
    ax = axes[idx // 2, idx % 2]
    
    if idx < 3:
        acc = generate_earthquake_acceleration(t, dt, omega_g, xi_g, S0, seed=idx*100)
        ax.plot(t, acc / 9.81, 'b-', linewidth=0.8)
        ax.set_ylabel('Acceleration (g)', fontsize=11)
        ax.set_title(f'Artificial Earthquake Sample {idx+1}', fontsize=12)
        
        # 添加统计信息
        rms = np.sqrt(np.mean(acc**2))
        pga = np.max(np.abs(acc))
        textstr = f'RMS: {rms/9.81:.3f}g\nPGA: {pga/9.81:.3f}g'
        ax.text(0.02, 0.95, textstr, transform=ax.transAxes, fontsize=9,
                verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    else:
        # 第4个子图显示平均反应谱
        n_periods = 50
        periods = np.linspace(0.1, 5.0, n_periods)
        spectra = np.zeros((n_samples, n_periods))
        
        for i in range(n_samples):
            acc = generate_earthquake_acceleration(t, dt, omega_g, xi_g, S0, seed=i*100)
            for j, T in enumerate(periods):
                omega_T = 2 * np.pi / T
                # 单自由度响应
                m_sdof = 1.0
                c_sdof = 2 * damping_ratio * omega_T * m_sdof
                k_sdof = omega_T**2 * m_sdof
                f_sdof = -acc * m_sdof
                u_sdof, _, a_sdof = newmark_beta_mdof(
                    np.array([[m_sdof]]), 
                    np.array([[c_sdof]]), 
                    np.array([[k_sdof]]), 
                    f_sdof.reshape(-1, 1), 
                    dt
                )
                spectra[i, j] = np.max(np.abs(a_sdof + acc.reshape(-1, 1)))
        
        mean_spectrum = np.mean(spectra, axis=0)
        ax.plot(periods, mean_spectrum / 9.81, 'b-', linewidth=2, label='Mean Spectrum')
        ax.plot(periods, np.percentile(spectra, 84, axis=0) / 9.81, 'r--', linewidth=1.5, label='84th Percentile')
        ax.set_xlabel('Period (s)', fontsize=11)
        ax.set_ylabel('Pseudo-Acceleration (g)', fontsize=11)
        ax.set_title('Mean Response Spectrum (5% Damping)', fontsize=12)
        ax.legend()
    
    ax.set_xlabel('Time (s)' if idx < 3 else 'Period (s)', fontsize=11)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, t_max if idx < 3 else 5.0)

plt.tight_layout()
plt.savefig('case2_earthquake_samples.png', dpi=150, bbox_inches='tight')
plt.close()

print("  已保存: case2_earthquake_samples.png")

# 创建动画：结构响应时程
print("\n正在生成结构响应动画...")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 子图1：变形动画
ax1 = axes[0]
ax1.set_xlim(-0.5, 0.5)
ax1.set_ylim(0, n_dof + 1)
ax1.set_xlabel('Displacement (m)', fontsize=12)
ax1.set_ylabel('Floor', fontsize=12)
ax1.set_title('Structure Deformation', fontsize=14)
ax1.grid(True, alpha=0.3)
ax1.axvline(0, color='k', linestyle='-', linewidth=0.5)

# 绘制楼层线
floor_lines = []
for i in range(n_dof):
    line, = ax1.plot([0, 0], [i+0.5, i+1.5], 'b-', linewidth=3)
    floor_lines.append(line)

# 子图2：顶层位移时程
ax2 = axes[1]
ax2.set_xlim(0, t_max)
ax2.set_ylim(np.min(u_mean_mc[:, -1]) * 1.5 * 1000, np.max(u_mean_mc[:, -1]) * 1.5 * 1000)
ax2.set_xlabel('Time (s)', fontsize=12)
ax2.set_ylabel('Displacement (mm)', fontsize=12)
ax2.set_title('Top Floor Displacement', fontsize=14)
ax2.grid(True, alpha=0.3)

line_disp, = ax2.plot([], [], 'b-', linewidth=2)
line_current, = ax2.plot([], [], 'ro', markersize=8)

def init():
    for line in floor_lines:
        line.set_data([0, 0], [0, 0])
    line_disp.set_data([], [])
    line_current.set_data([], [])
    return floor_lines + [line_disp, line_current]

def update(frame):
    # 更新变形
    scale = 50  # 放大系数
    for i, line in enumerate(floor_lines):
        disp = u_mean_mc[frame, i] * scale
        line.set_data([0 + disp, 0 + disp], [i+0.5, i+1.5])
    
    # 更新时程曲线
    line_disp.set_data(t[:frame+1], u_mean_mc[:frame+1, -1] * 1000)
    line_current.set_data([t[frame]], [u_mean_mc[frame, -1] * 1000])
    
    return floor_lines + [line_disp, line_current]

# 每3帧保存一次
frame_indices = np.arange(0, n_steps, 3)
anim = FuncAnimation(fig, update, frames=frame_indices, init_func=init, blit=True)

writer = PillowWriter(fps=30)
anim.save('case2_structure_response.gif', writer=writer)
plt.close()

print("  已保存: case2_structure_response.gif")

print("\n" + "=" * 60)
print("案例2分析完成!")
print("=" * 60)

7.3 案例3：蒙特卡洛模拟在结构可靠性分析中的应用

使用蒙特卡洛方法分析一个结构的失效概率。比较不同样本量下的计算精度和效率，并应用方差缩减技术。

"""
案例3：蒙特卡洛模拟在结构可靠性分析中的应用
比较不同样本量下的计算精度和效率，并应用方差缩减技术
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from scipy import stats
from scipy.integrate import trapezoid
import time
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


def sdof_response(m, c, k, f, dt):
    """单自由度系统动力响应计算"""
    n_steps = len(f)
    u = np.zeros(n_steps)
    v = np.zeros(n_steps)
    a = np.zeros(n_steps)
    
    a[0] = f[0] / m
    gamma, beta = 0.5, 0.25
    k_eff = k + gamma * c / (beta * dt) + m / (beta * dt**2)
    
    for i in range(n_steps - 1):
        f_eff = (f[i+1] + 
                 m * (u[i] / (beta * dt**2) + v[i] / (beta * dt) + (1/(2*beta)-1) * a[i]) +
                 c * (gamma * u[i] / (beta * dt) + (gamma/beta - 1) * v[i] + 
                      dt * (gamma/(2*beta) - 1) * a[i]))
        u[i+1] = f_eff / k_eff
        a[i+1] = (u[i+1] - u[i]) / (beta * dt**2) - v[i] / (beta * dt) - (1/(2*beta)-1) * a[i]
        v[i+1] = v[i] + (1-gamma) * dt * a[i] + gamma * dt * a[i+1]
    
    return u, v, a


def crude_monte_carlo(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples):
    """ crude蒙特卡洛模拟 """
    max_displacements = np.zeros(n_samples)
    
    zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
    lambda_m = np.log(m0) - 0.5 * zeta_m**2
    zeta_c = np.sqrt(np.log(1 + (sigma_c/c0)**2))
    lambda_c = np.log(c0) - 0.5 * zeta_c**2
    zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
    lambda_k = np.log(k0) - 0.5 * zeta_k**2
    
    for i in range(n_samples):
        m = np.random.lognormal(lambda_m, zeta_m)
        c = np.random.lognormal(lambda_c, zeta_c)
        k = np.random.lognormal(lambda_k, zeta_k)
        u, _, _ = sdof_response(m, c, k, f, dt)
        max_displacements[i] = np.max(np.abs(u))
    
    return max_displacements


def importance_sampling(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples, threshold):
    """重要性采样 """
    max_displacements = np.zeros(n_samples)
    weights = np.ones(n_samples)
    
    # 调整采样分布，增加大位移的概率
    shift_factor = 1.5  # 向刚度降低方向偏移
    
    zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
    lambda_m = np.log(m0) - 0.5 * zeta_m**2
    zeta_c = np.sqrt(np.log(1 + (sigma_c/c0)**2))
    lambda_c = np.log(c0) - 0.5 * zeta_c**2
    zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
    lambda_k = np.log(k0) - 0.5 * zeta_k**2
    
    # 重要性采样分布参数（降低刚度均值）
    lambda_k_is = lambda_k - 0.5
    
    for i in range(n_samples):
        m = np.random.lognormal(lambda_m, zeta_m)
        c = np.random.lognormal(lambda_c, zeta_c)
        k = np.random.lognormal(lambda_k_is, zeta_k)
        
        u, _, _ = sdof_response(m, c, k, f, dt)
        max_displacements[i] = np.max(np.abs(u))
        
        # 计算权重（原始PDF / 采样PDF）
        pdf_original = stats.lognorm.pdf(k, s=zeta_k, scale=np.exp(lambda_k))
        pdf_sampling = stats.lognorm.pdf(k, s=zeta_k, scale=np.exp(lambda_k_is))
        weights[i] = pdf_original / pdf_sampling
    
    return max_displacements, weights


def latin_hypercube_sampling_mc(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples):
    """拉丁超立方采样蒙特卡洛 """
    max_displacements = np.zeros(n_samples)
    
    # LHS采样
    def lhs_normal(n_samples, n_vars):
        from scipy.special import erfcinv
        samples = np.zeros((n_samples, n_vars))
        for i in range(n_vars):
            strata = np.random.permutation(n_samples)
            for j in range(n_samples):
                u = (strata[j] + np.random.rand()) / n_samples
                samples[j, i] = np.sqrt(2) * erfcinv(2 * (1 - u))
        return samples
    
    lhs_samples = lhs_normal(n_samples, 3)
    
    zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
    lambda_m = np.log(m0) - 0.5 * zeta_m**2
    zeta_c = np.sqrt(np.log(1 + (sigma_c/c0)**2))
    lambda_c = np.log(c0) - 0.5 * zeta_c**2
    zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
    lambda_k = np.log(k0) - 0.5 * zeta_k**2
    
    for i in range(n_samples):
        m = np.exp(lambda_m + zeta_m * lhs_samples[i, 0])
        c = np.exp(lambda_c + zeta_c * lhs_samples[i, 1])
        k = np.exp(lambda_k + zeta_k * lhs_samples[i, 2])
        u, _, _ = sdof_response(m, c, k, f, dt)
        max_displacements[i] = np.max(np.abs(u))
    
    return max_displacements


# ==================== 主程序 ====================
print("=" * 60)
print("案例3：蒙特卡洛模拟在结构可靠性分析中的应用")
print("=" * 60)

# 系统参数
m0 = 1000.0
c0 = 2000.0
k0 = 100000.0
cov_m, cov_c, cov_k = 0.1, 0.15, 0.1
sigma_m = m0 * cov_m
sigma_c = c0 * cov_c
sigma_k = k0 * cov_k

# 荷载条件
dt = 0.01
t_max = 10.0
t = np.arange(0, t_max, dt)
f = np.zeros(len(t))
pulse_start, pulse_end = int(0.5/dt), int(1.0/dt)
f[pulse_start:pulse_end] = 10000.0

# 失效阈值
threshold = 0.25  # 最大允许位移 (m)

print(f"\n系统参数:")
print(f"  质量: {m0:.0f} kg (COV={cov_m*100:.0f}%)")
print(f"  阻尼: {c0:.0f} N·s/m (COV={cov_c*100:.0f}%)")
print(f"  刚度: {k0:.0f} N/m (COV={cov_k*100:.0f}%)")
print(f"\n失效阈值: 最大位移 > {threshold*1000:.0f} mm")

# ==================== 样本量收敛性分析 ====================
print("\n" + "=" * 60)
print("样本量收敛性分析")
print("=" * 60)

sample_sizes = [100, 500, 1000, 2000, 5000, 10000]
results_convergence = []

np.random.seed(42)

for n in sample_sizes:
    start_time = time.time()
    max_disp = crude_monte_carlo(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n)
    elapsed = time.time() - start_time
    
    pf = np.sum(max_disp > threshold) / n
    cov_pf = np.sqrt((1 - pf) / (pf * n)) if pf > 0 else 0
    
    results_convergence.append({
        'n_samples': n,
        'pf': pf,
        'cov': cov_pf,
        'time': elapsed,
        'mean': np.mean(max_disp),
        'std': np.std(max_disp)
    })
    
    print(f"n={n:5d}: Pf={pf:.4f}, COV={cov_pf:.3f}, Time={elapsed:.2f}s")

# ==================== 方差缩减技术对比 ====================
print("\n" + "=" * 60)
print("方差缩减技术对比")
print("=" * 60)

n_samples_comparison = 5000

# Crude MC
np.random.seed(123)
start = time.time()
max_disp_crude = crude_monte_carlo(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples_comparison)
time_crude = time.time() - start
pf_crude = np.sum(max_disp_crude > threshold) / n_samples_comparison
var_crude = np.var(max_disp_crude)

# LHS
np.random.seed(123)
start = time.time()
max_disp_lhs = latin_hypercube_sampling_mc(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples_comparison)
time_lhs = time.time() - start
pf_lhs = np.sum(max_disp_lhs > threshold) / n_samples_comparison
var_lhs = np.var(max_disp_lhs)

# Importance Sampling
np.random.seed(123)
start = time.time()
max_disp_is, weights_is = importance_sampling(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, n_samples_comparison, threshold)
time_is = time.time() - start

# 加权失效概率
failures_is = max_disp_is > threshold
pf_is = np.sum(weights_is[failures_is]) / np.sum(weights_is)
var_is = np.var(max_disp_is * weights_is) / np.mean(weights_is)**2

print(f"\n方法对比 (n={n_samples_comparison}):")
print(f"{'方法':<20} {'失效概率':<12} {'方差':<12} {'计算时间(s)':<12} {'效率提升'}")
print("-" * 70)
print(f"{'Crude MC':<20} {pf_crude:<12.4f} {var_crude:<12.6f} {time_crude:<12.2f} {'1.0x'}")
print(f"{'LHS':<20} {pf_lhs:<12.4f} {var_lhs:<12.6f} {time_lhs:<12.2f} {var_crude/var_lhs:.2f}x")
print(f"{'Importance Sampling':<20} {pf_is:<12.4f} {var_is:<12.6f} {time_is:<12.2f} {var_crude/var_is:.2f}x")

# ==================== 绘制结果 ====================
print("\n正在生成可视化结果...")

# 图1：收敛性分析
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

n_vals = [r['n_samples'] for r in results_convergence]
pf_vals = [r['pf'] for r in results_convergence]
cov_vals = [r['cov'] for r in results_convergence]
time_vals = [r['time'] for r in results_convergence]

ax = axes[0, 0]
ax.semilogx(n_vals, pf_vals, 'bo-', linewidth=2, markersize=8)
ax.axhline(y=pf_crude, color='r', linestyle='--', label=f'Reference Pf={pf_crude:.4f}')
ax.set_xlabel('Number of Samples', fontsize=11)
ax.set_ylabel('Failure Probability', fontsize=11)
ax.set_title('Convergence of Failure Probability', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

ax = axes[0, 1]
ax.loglog(n_vals, cov_vals, 'rs-', linewidth=2, markersize=8)
ax.set_xlabel('Number of Samples', fontsize=11)
ax.set_ylabel('Coefficient of Variation', fontsize=11)
ax.set_title('COV of Pf Estimate vs Sample Size', fontsize=12)
ax.grid(True, alpha=0.3)

ax = axes[1, 0]
ax.plot(n_vals, time_vals, 'g^-', linewidth=2, markersize=8)
ax.set_xlabel('Number of Samples', fontsize=11)
ax.set_ylabel('Computation Time (s)', fontsize=11)
ax.set_title('Computation Time vs Sample Size', fontsize=12)
ax.grid(True, alpha=0.3)

ax = axes[1, 1]
ax.loglog(n_vals, [c*t for c, t in zip(cov_vals, time_vals)], 'mv-', linewidth=2, markersize=8)
ax.set_xlabel('Number of Samples', fontsize=11)
ax.set_ylabel('COV × Time (effort)', fontsize=11)
ax.set_title('Computational Effort vs Sample Size', fontsize=12)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case3_convergence_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: case3_convergence_analysis.png")

# 图2：方法对比
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 响应分布对比
ax = axes[0, 0]
ax.hist(max_disp_crude, bins=50, density=True, alpha=0.5, label='Crude MC', color='blue')
ax.hist(max_disp_lhs, bins=50, density=True, alpha=0.5, label='LHS', color='green')
ax.axvline(threshold, color='r', linestyle='--', linewidth=2, label=f'Threshold={threshold}m')
ax.set_xlabel('Maximum Displacement (m)', fontsize=11)
ax.set_ylabel('Probability Density', fontsize=11)
ax.set_title('Response Distribution Comparison', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 重要性采样权重分布
ax = axes[0, 1]
ax.hist(weights_is, bins=50, density=True, alpha=0.7, color='orange', edgecolor='black')
ax.set_xlabel('Importance Weight', fontsize=11)
ax.set_ylabel('Probability Density', fontsize=11)
ax.set_title('Importance Sampling Weights Distribution', fontsize=12)
ax.grid(True, alpha=0.3)

# 方法效率对比
ax = axes[1, 0]
methods = ['Crude MC', 'LHS', 'Importance\nSampling']
efficiency = [1.0, var_crude/var_lhs, var_crude/var_is]
colors = ['blue', 'green', 'orange']
bars = ax.bar(methods, efficiency, color=colors, alpha=0.7)
ax.set_ylabel('Variance Reduction Factor', fontsize=11)
ax.set_title('Variance Reduction Efficiency', fontsize=12)
ax.grid(True, alpha=0.3, axis='y')
for bar, eff in zip(bars, efficiency):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, 
            f'{eff:.2f}x', ha='center', fontsize=11)

# 失效概率置信区间
ax = axes[1, 1]
pf_values = [pf_crude, pf_lhs, pf_is]
errors = [1.96*np.sqrt(pf_crude*(1-pf_crude)/n_samples_comparison),
          1.96*np.sqrt(pf_lhs*(1-pf_lhs)/n_samples_comparison),
          1.96*np.sqrt(pf_is*(1-pf_is)/n_samples_comparison)]
x_pos = np.arange(len(methods))
ax.bar(x_pos, pf_values, yerr=errors, capsize=5, color=colors, alpha=0.7)
ax.set_xticks(x_pos)
ax.set_xticklabels(methods)
ax.set_ylabel('Failure Probability', fontsize=11)
ax.set_title('Failure Probability with 95% Confidence Interval', fontsize=12)
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('case3_method_comparison.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: case3_method_comparison.png")

# 图3：响应统计量
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 累积分布函数
ax = axes[0, 0]
sorted_crude = np.sort(max_disp_crude)
sorted_lhs = np.sort(max_disp_lhs)
y_vals = np.arange(1, len(sorted_crude)+1) / len(sorted_crude)
ax.plot(sorted_crude, y_vals, 'b-', linewidth=2, label='Crude MC')
ax.plot(sorted_lhs, y_vals, 'g--', linewidth=2, label='LHS')
ax.axvline(threshold, color='r', linestyle=':', linewidth=2, label='Threshold')
ax.set_xlabel('Maximum Displacement (m)', fontsize=11)
ax.set_ylabel('Cumulative Probability', fontsize=11)
ax.set_title('Cumulative Distribution Function', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 概率纸图
ax = axes[0, 1]
# 转换为标准正态分位数
z_crude = stats.norm.ppf(y_vals)
z_lhs = stats.norm.ppf(y_vals)
ax.plot(sorted_crude, z_crude, 'b.', markersize=3, alpha=0.5, label='Crude MC')
ax.plot(sorted_lhs, z_lhs, 'g.', markersize=3, alpha=0.5, label='LHS')
ax.axvline(threshold, color='r', linestyle='--', linewidth=2)
ax.set_xlabel('Maximum Displacement (m)', fontsize=11)
ax.set_ylabel('Standard Normal Quantile', fontsize=11)
ax.set_title('Normal Probability Plot', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 分位数对比
ax = axes[1, 0]
quantiles = [0.5, 0.75, 0.9, 0.95, 0.99]
q_crude = [np.percentile(max_disp_crude, q*100) for q in quantiles]
q_lhs = [np.percentile(max_disp_lhs, q*100) for q in quantiles]
x = np.arange(len(quantiles))
width = 0.35
ax.bar(x - width/2, q_crude, width, label='Crude MC', color='blue', alpha=0.7)
ax.bar(x + width/2, q_lhs, width, label='LHS', color='green', alpha=0.7)
ax.set_xticks(x)
ax.set_xticklabels([f'{int(q*100)}%' for q in quantiles])
ax.set_ylabel('Displacement (m)', fontsize=11)
ax.set_title('Quantile Comparison', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3, axis='y')

# 统计量汇总
ax = axes[1, 1]
ax.axis('off')
stats_text = "Statistical Summary (Crude MC):\n\n"
stats_text += f"Mean: {np.mean(max_disp_crude):.4f} m\n"
stats_text += f"Std Dev: {np.std(max_disp_crude):.4f} m\n"
stats_text += f"COV: {np.std(max_disp_crude)/np.mean(max_disp_crude)*100:.2f}%\n"
stats_text += f"Min: {np.min(max_disp_crude):.4f} m\n"
stats_text += f"Max: {np.max(max_disp_crude):.4f} m\n\n"
stats_text += f"Failure Probability: {pf_crude:.4f}\n"
stats_text += f"95% CI: [{pf_crude-1.96*np.sqrt(pf_crude*(1-pf_crude)/n_samples_comparison):.4f}, "
stats_text += f"{pf_crude+1.96*np.sqrt(pf_crude*(1-pf_crude)/n_samples_comparison):.4f}]\n\n"
stats_text += f"Reliability Index: {-stats.norm.ppf(pf_crude):.3f}"

ax.text(0.1, 0.5, stats_text, fontsize=11, verticalalignment='center',
        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.tight_layout()
plt.savefig('case3_response_statistics.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: case3_response_statistics.png")

# 创建动画：样本收敛过程
print("\n正在生成收敛动画...")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

ax1 = axes[0]
ax1.set_xlim(0, np.max(max_disp_crude) * 1.1)
ax1.set_ylim(0, 15)
ax1.set_xlabel('Maximum Displacement (m)', fontsize=12)
ax1.set_ylabel('Probability Density', fontsize=12)
ax1.set_title('Response Distribution Evolution', fontsize=14)
ax1.grid(True, alpha=0.3)
ax1.axvline(threshold, color='r', linestyle='--', linewidth=2, label='Threshold')

ax2 = axes[1]
ax2.set_xlim(0, n_samples_comparison)
ax2.set_ylim(0, max(pf_crude * 2, 0.1))
ax2.set_xlabel('Number of Samples', fontsize=12)
ax2.set_ylabel('Failure Probability', fontsize=12)
ax2.set_title('Failure Probability Convergence', fontsize=14)
ax2.grid(True, alpha=0.3)
ax2.axhline(y=pf_crude, color='r', linestyle='--', linewidth=2, label='Final Pf')

line_pdf, = ax1.plot([], [], 'b-', linewidth=2)
line_pf, = ax2.plot([], [], 'g-', linewidth=2)

# 预计算
sample_steps = np.logspace(2, np.log10(n_samples_comparison), 50).astype(int)
sample_steps = np.unique(sample_steps)

pf_history = []
for n in sample_steps:
    pf_n = np.sum(max_disp_crude[:n] > threshold) / n
    pf_history.append(pf_n)

def init():
    line_pdf.set_data([], [])
    line_pf.set_data([], [])
    return [line_pdf, line_pf]

def update(frame):
    n = sample_steps[frame]
    
    # 更新PDF
    data = max_disp_crude[:n]
    hist, bins = np.histogram(data, bins=30, density=True)
    bin_centers = (bins[:-1] + bins[1:]) / 2
    line_pdf.set_data(bin_centers, hist)
    
    # 更新Pf
    ax2.clear()
    ax2.set_xlim(0, n_samples_comparison)
    ax2.set_ylim(0, max(max(pf_history), 0.1))
    ax2.set_xlabel('Number of Samples', fontsize=12)
    ax2.set_ylabel('Failure Probability', fontsize=12)
    ax2.set_title(f'Failure Probability Convergence (n={n})', fontsize=14)
    ax2.grid(True, alpha=0.3)
    ax2.axhline(y=pf_crude, color='r', linestyle='--', linewidth=2, label='Final Pf')
    ax2.plot(sample_steps[:frame+1], pf_history[:frame+1], 'g-', linewidth=2)
    ax2.legend()
    
    return [line_pdf]

anim = FuncAnimation(fig, update, frames=len(sample_steps), init_func=init)
writer = PillowWriter(fps=5)
anim.save('case3_convergence_animation.gif', writer=writer)
plt.close()
print("  已保存: case3_convergence_animation.gif")

print("\n" + "=" * 60)
print("案例3分析完成!")
print("=" * 60)

7.4 案例4：基于FORM的可靠性指标计算

使用一次二阶矩法计算结构的可靠指标。进行灵敏度分析，识别对结构可靠性影响最大的参数。

"""
案例4：基于FORM的可靠性指标计算与灵敏度分析
使用一次二阶矩法计算结构的可靠指标，识别对结构可靠性影响最大的参数
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from scipy import stats, optimize
from scipy.integrate import trapezoid
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


def sdof_response_max(m, c, k, f, dt):
    """计算单自由度系统最大位移响应"""
    n_steps = len(f)
    u = np.zeros(n_steps)
    v = np.zeros(n_steps)
    a = np.zeros(n_steps)
    
    a[0] = f[0] / m
    gamma, beta = 0.5, 0.25
    k_eff = k + gamma * c / (beta * dt) + m / (beta * dt**2)
    
    for i in range(n_steps - 1):
        f_eff = (f[i+1] + 
                 m * (u[i] / (beta * dt**2) + v[i] / (beta * dt) + (1/(2*beta)-1) * a[i]) +
                 c * (gamma * u[i] / (beta * dt) + (gamma/beta - 1) * v[i] + 
                      dt * (gamma/(2*beta) - 1) * a[i]))
        u[i+1] = f_eff / k_eff
        a[i+1] = (u[i+1] - u[i]) / (beta * dt**2) - v[i] / (beta * dt) - (1/(2*beta)-1) * a[i]
        v[i+1] = v[i] + (1-gamma) * dt * a[i] + gamma * dt * a[i+1]
    
    return np.max(np.abs(u))


def limit_state_function(x, f, dt, threshold):
    """
    功能函数
    x = [ln(m), ln(c), ln(k)] - 对数空间的标准正态变量
    """
    # 转换回物理参数
    m = np.exp(x[0])
    c = np.exp(x[1])
    k = np.exp(x[2])
    
    # 计算最大位移
    max_disp = sdof_response_max(m, c, k, f, dt)
    
    # 功能函数：阈值 - 最大位移
    g = threshold - max_disp
    
    return g


def form_analysis(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, tol=1e-6, max_iter=100):
    """
    一次二阶矩法（FORM）计算可靠指标
    
    参数:
        m0, c0, k0: 名义参数
        sigma_m, sigma_c, sigma_k: 参数标准差
        f, dt: 荷载参数
        threshold: 失效阈值
    
    返回:
        beta: 可靠指标
        u_star: 验算点（标准正态空间）
        alpha: 方向余弦（重要性因子）
    """
    # 对数正态分布参数
    zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
    lambda_m = np.log(m0) - 0.5 * zeta_m**2
    
    zeta_c = np.sqrt(np.log(1 + (sigma_c/c0)**2))
    lambda_c = np.log(c0) - 0.5 * zeta_c**2
    
    zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
    lambda_k = np.log(k0) - 0.5 * zeta_k**2
    
    # 初始验算点（均值点）
    u = np.array([0.0, 0.0, 0.0])
    
    # 迭代求解
    for iteration in range(max_iter):
        # 转换到物理空间
        x = np.array([lambda_m + zeta_m * u[0],
                      lambda_c + zeta_c * u[1],
                      lambda_k + zeta_k * u[2]])
        
        # 计算功能函数值
        g_val = limit_state_function(x, f, dt, threshold)
        
        # 数值计算梯度
        eps = 0.01
        grad = np.zeros(3)
        for i in range(3):
            x_plus = x.copy()
            x_plus[i] += eps
            x_minus = x.copy()
            x_minus[i] -= eps
            grad[i] = (limit_state_function(x_plus, f, dt, threshold) - 
                       limit_state_function(x_minus, f, dt, threshold)) / (2 * eps)
        
        # 链式法则：dg/du = dg/dx * dx/du
        grad_u = grad * np.array([zeta_m, zeta_c, zeta_k])
        
        # 方向余弦
        grad_norm = np.linalg.norm(grad_u)
        if grad_norm < 1e-10:
            print("Warning: Gradient too small")
            break
        
        alpha = -grad_u / grad_norm
        
        # 更新验算点
        beta_new = (g_val - np.dot(grad_u, u)) / grad_norm
        u_new = beta_new * alpha
        
        # 检查收敛
        if np.linalg.norm(u_new - u) < tol:
            u = u_new
            beta = beta_new
            break
        
        u = u_new
        beta = beta_new
    
    # 计算最终的方向余弦
    x = np.array([lambda_m + zeta_m * u[0],
                  lambda_c + zeta_c * u[1],
                  lambda_k + zeta_k * u[2]])
    
    eps = 0.01
    grad = np.zeros(3)
    for i in range(3):
        x_plus = x.copy()
        x_plus[i] += eps
        x_minus = x.copy()
        x_minus[i] -= eps
        grad[i] = (limit_state_function(x_plus, f, dt, threshold) - 
                   limit_state_function(x_minus, f, dt, threshold)) / (2 * eps)
    
    grad_u = grad * np.array([zeta_m, zeta_c, zeta_k])
    grad_norm = np.linalg.norm(grad_u)
    alpha = -grad_u / grad_norm if grad_norm > 1e-10 else np.array([0, 0, 0])
    
    return beta, u, alpha, iteration + 1


def monte_carlo_reliability(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=10000):
    """蒙特卡洛模拟计算失效概率"""
    failures = 0
    max_disps = np.zeros(n_samples)
    
    zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
    lambda_m = np.log(m0) - 0.5 * zeta_m**2
    zeta_c = np.sqrt(np.log(1 + (sigma_c/c0)**2))
    lambda_c = np.log(c0) - 0.5 * zeta_c**2
    zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
    lambda_k = np.log(k0) - 0.5 * zeta_k**2
    
    for i in range(n_samples):
        m = np.random.lognormal(lambda_m, zeta_m)
        c = np.random.lognormal(lambda_c, zeta_c)
        k = np.random.lognormal(lambda_k, zeta_k)
        
        max_disp = sdof_response_max(m, c, k, f, dt)
        max_disps[i] = max_disp
        
        if max_disp > threshold:
            failures += 1
    
    pf = failures / n_samples
    beta_mc = -stats.norm.ppf(pf) if pf > 0 else 10.0
    
    return pf, beta_mc, max_disps


def sensitivity_analysis(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold):
    """
    灵敏度分析
    """
    # 基准失效概率
    pf_base, _, _ = monte_carlo_reliability(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=5000)
    
    # 参数扰动
    delta = 0.1
    sensitivities = {}
    
    # 质量灵敏度
    pf_m_plus, _, _ = monte_carlo_reliability(m0*(1+delta), c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=3000)
    pf_m_minus, _, _ = monte_carlo_reliability(m0*(1-delta), c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=3000)
    sensitivities['mass'] = (pf_m_plus - pf_m_minus) / (2 * delta * m0) * m0 / pf_base if pf_base > 0 else 0
    
    # 阻尼灵敏度
    pf_c_plus, _, _ = monte_carlo_reliability(m0, c0*(1+delta), k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=3000)
    pf_c_minus, _, _ = monte_carlo_reliability(m0, c0*(1-delta), k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=3000)
    sensitivities['damping'] = (pf_c_plus - pf_c_minus) / (2 * delta * c0) * c0 / pf_base if pf_base > 0 else 0
    
    # 刚度灵敏度
    pf_k_plus, _, _ = monte_carlo_reliability(m0, c0, k0*(1+delta), sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=3000)
    pf_k_minus, _, _ = monte_carlo_reliability(m0, c0, k0*(1-delta), sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=3000)
    sensitivities['stiffness'] = (pf_k_plus - pf_k_minus) / (2 * delta * k0) * k0 / pf_base if pf_base > 0 else 0
    
    return sensitivities


# ==================== 主程序 ====================
print("=" * 60)
print("案例4：基于FORM的可靠性指标计算与灵敏度分析")
print("=" * 60)

# 系统参数
m0 = 1000.0
c0 = 2000.0
k0 = 100000.0
cov_m, cov_c, cov_k = 0.1, 0.15, 0.1
sigma_m = m0 * cov_m
sigma_c = c0 * cov_c
sigma_k = k0 * cov_k

# 荷载条件
dt = 0.01
t_max = 10.0
t = np.arange(0, t_max, dt)
f = np.zeros(len(t))
pulse_start, pulse_end = int(0.5/dt), int(1.0/dt)
f[pulse_start:pulse_end] = 10000.0

# 失效阈值
threshold = 0.25

print(f"\n系统参数:")
print(f"  质量: {m0:.0f} kg (COV={cov_m*100:.0f}%)")
print(f"  阻尼: {c0:.0f} N·s/m (COV={cov_c*100:.0f}%)")
print(f"  刚度: {k0:.0f} N/m (COV={cov_k*100:.0f}%)")
print(f"\n失效阈值: 最大位移 > {threshold*1000:.0f} mm")

# ==================== FORM分析 ====================
print("\n" + "=" * 60)
print("FORM分析")
print("=" * 60)

beta_form, u_star, alpha, n_iter = form_analysis(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold)
pf_form = stats.norm.cdf(-beta_form)

print(f"\nFORM结果:")
print(f"  迭代次数: {n_iter}")
print(f"  可靠指标 β: {beta_form:.4f}")
print(f"  失效概率 Pf: {pf_form:.6f}")

print(f"\n验算点（标准正态空间）:")
print(f"  u_m* = {u_star[0]:.4f}")
print(f"  u_c* = {u_star[1]:.4f}")
print(f"  u_k* = {u_star[2]:.4f}")

print(f"\n方向余弦（重要性因子）:")
print(f"  α_m = {alpha[0]:.4f} (α² = {alpha[0]**2*100:.2f}%)")
print(f"  α_c = {alpha[1]:.4f} (α² = {alpha[1]**2*100:.2f}%)")
print(f"  α_k = {alpha[2]:.4f} (α² = {alpha[2]**2*100:.2f}%)")

# 转换到物理空间
zeta_m = np.sqrt(np.log(1 + (sigma_m/m0)**2))
lambda_m = np.log(m0) - 0.5 * zeta_m**2
zeta_c = np.sqrt(np.log(1 + (sigma_c/c0)**2))
lambda_c = np.log(c0) - 0.5 * zeta_c**2
zeta_k = np.sqrt(np.log(1 + (sigma_k/k0)**2))
lambda_k = np.log(k0) - 0.5 * zeta_k**2

m_star = np.exp(lambda_m + zeta_m * u_star[0])
c_star = np.exp(lambda_c + zeta_c * u_star[1])
k_star = np.exp(lambda_k + zeta_k * u_star[2])

print(f"\n验算点（物理空间）:")
print(f"  m* = {m_star:.2f} kg")
print(f"  c* = {c_star:.2f} N·s/m")
print(f"  k* = {k_star:.2f} N/m")

# ==================== 蒙特卡洛验证 ====================
print("\n" + "=" * 60)
print("蒙特卡洛验证")
print("=" * 60)

pf_mc, beta_mc, max_disps = monte_carlo_reliability(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=10000)

print(f"\n蒙特卡洛结果:")
print(f"  失效概率 Pf: {pf_mc:.6f}")
print(f"  可靠指标 β: {beta_mc:.4f}")

print(f"\n方法对比:")
print(f"  FORM: β = {beta_form:.4f}, Pf = {pf_form:.6f}")
print(f"  MC:   β = {beta_mc:.4f}, Pf = {pf_mc:.6f}")
print(f"  相对误差: {abs(beta_form - beta_mc) / beta_mc * 100:.2f}%")

# ==================== 绘制结果 ====================
print("\n正在生成可视化结果...")

# 图1：FORM结果可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：标准正态空间中的验算点
ax = axes[0, 0]
# 绘制单位球
phi = np.linspace(0, 2*np.pi, 100)
ax.plot(np.cos(phi), np.sin(phi), 'k--', linewidth=1, alpha=0.5, label='Unit Sphere')
# 绘制验算点
ax.plot(u_star[0], u_star[1], 'ro', markersize=12, label='Design Point')
ax.arrow(0, 0, u_star[0], u_star[1], head_width=0.1, head_length=0.1, fc='red', ec='red', alpha=0.7)
ax.set_xlabel('u_m (Mass)', fontsize=11)
ax.set_ylabel('u_c (Damping)', fontsize=11)
ax.set_title('Design Point in Standard Normal Space (u_m - u_c)', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)
ax.axis('equal')

ax = axes[0, 1]
ax.plot(np.cos(phi), np.sin(phi), 'k--', linewidth=1, alpha=0.5, label='Unit Sphere')
ax.plot(u_star[0], u_star[2], 'ro', markersize=12, label='Design Point')
ax.arrow(0, 0, u_star[0], u_star[2], head_width=0.1, head_length=0.1, fc='red', ec='red', alpha=0.7)
ax.set_xlabel('u_m (Mass)', fontsize=11)
ax.set_ylabel('u_k (Stiffness)', fontsize=11)
ax.set_title('Design Point in Standard Normal Space (u_m - u_k)', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)
ax.axis('equal')

# 子图2：重要性因子
ax = axes[1, 0]
params = ['Mass', 'Damping', 'Stiffness']
importance = [alpha[0]**2, alpha[1]**2, alpha[2]**2]
colors = ['blue', 'green', 'orange']
bars = ax.bar(params, importance, color=colors, alpha=0.7)
ax.set_ylabel('Importance Factor α²', fontsize=11)
ax.set_title('Parameter Importance Factors', fontsize=12)
ax.grid(True, alpha=0.3, axis='y')
for bar, imp in zip(bars, importance):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01, 
            f'{imp*100:.1f}%', ha='center', fontsize=11)

# 子图3：可靠性指标对比
ax = axes[1, 1]
methods = ['FORM', 'Monte Carlo']
betas = [beta_form, beta_mc]
colors = ['blue', 'red']
bars = ax.bar(methods, betas, color=colors, alpha=0.7)
ax.set_ylabel('Reliability Index β', fontsize=11)
ax.set_title('Reliability Index Comparison', fontsize=12)
ax.grid(True, alpha=0.3, axis='y')
for bar, beta in zip(bars, betas):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.05, 
            f'{beta:.3f}', ha='center', fontsize=11)

plt.tight_layout()
plt.savefig('case4_form_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: case4_form_analysis.png")

# 图2：失效概率分析
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：响应分布与失效区域
ax = axes[0, 0]
ax.hist(max_disps, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
ax.axvline(threshold, color='r', linestyle='--', linewidth=2, label=f'Threshold = {threshold}m')
ax.axvline(np.mean(max_disps), color='g', linestyle=':', linewidth=2, label=f'Mean = {np.mean(max_disps):.3f}m')
ax.fill_betweenx([0, ax.get_ylim()[1]], threshold, ax.get_xlim()[1], alpha=0.3, color='red', label='Failure Region')
ax.set_xlabel('Maximum Displacement (m)', fontsize=11)
ax.set_ylabel('Probability Density', fontsize=11)
ax.set_title('Response Distribution and Failure Region', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图2：累积分布函数
ax = axes[0, 1]
sorted_disps = np.sort(max_disps)
y_vals = np.arange(1, len(sorted_disps)+1) / len(sorted_disps)
ax.plot(sorted_disps, y_vals, 'b-', linewidth=2)
ax.axvline(threshold, color='r', linestyle='--', linewidth=2)
ax.axhline(pf_mc, color='r', linestyle=':', linewidth=2, label=f'Pf = {pf_mc:.4f}')
ax.set_xlabel('Maximum Displacement (m)', fontsize=11)
ax.set_ylabel('Cumulative Probability', fontsize=11)
ax.set_title('Cumulative Distribution Function', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图3：不同阈值下的可靠指标
ax = axes[1, 0]
thresholds = np.linspace(0.15, 0.35, 20)
betas_form_range = []
betas_mc_range = []

for thresh in thresholds:
    beta_f, _, _, _ = form_analysis(m0, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, thresh)
    betas_form_range.append(beta_f)
    
    # 简化的MC估计
    failures = np.sum(max_disps > thresh)
    pf = failures / len(max_disps)
    beta_m = -stats.norm.ppf(pf) if pf > 0 else 5.0
    betas_mc_range.append(beta_m)

ax.plot(thresholds*1000, betas_form_range, 'b-', linewidth=2, label='FORM')
ax.plot(thresholds*1000, betas_mc_range, 'r--', linewidth=2, label='Monte Carlo')
ax.set_xlabel('Threshold (mm)', fontsize=11)
ax.set_ylabel('Reliability Index β', fontsize=11)
ax.set_title('Reliability Index vs Threshold', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图4：失效概率vs阈值
ax = axes[1, 1]
pf_form_range = [stats.norm.cdf(-beta) for beta in betas_form_range]
pf_mc_range = [stats.norm.cdf(-beta) for beta in betas_mc_range]
ax.semilogy(thresholds*1000, pf_form_range, 'b-', linewidth=2, label='FORM')
ax.semilogy(thresholds*1000, pf_mc_range, 'r--', linewidth=2, label='Monte Carlo')
ax.set_xlabel('Threshold (mm)', fontsize=11)
ax.set_ylabel('Failure Probability Pf', fontsize=11)
ax.set_title('Failure Probability vs Threshold', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case4_failure_probability.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: case4_failure_probability.png")

# 图3：参数敏感性分析
print("\n正在生成参数敏感性分析...")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：不同质量下的失效概率
ax = axes[0, 0]
m_range = np.linspace(m0*0.7, m0*1.3, 20)
pf_m_range = []
for m in m_range:
    pf, _, _ = monte_carlo_reliability(m, c0, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=2000)
    pf_m_range.append(pf)

ax.plot(m_range, pf_m_range, 'b-', linewidth=2)
ax.axvline(m0, color='r', linestyle='--', linewidth=2, label=f'Nominal m = {m0}')
ax.axvline(m_star, color='g', linestyle=':', linewidth=2, label=f'Design m* = {m_star:.1f}')
ax.set_xlabel('Mass (kg)', fontsize=11)
ax.set_ylabel('Failure Probability', fontsize=11)
ax.set_title('Failure Probability vs Mass', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图2：不同阻尼下的失效概率
ax = axes[0, 1]
c_range = np.linspace(c0*0.5, c0*2.0, 20)
pf_c_range = []
for c in c_range:
    pf, _, _ = monte_carlo_reliability(m0, c, k0, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=2000)
    pf_c_range.append(pf)

ax.plot(c_range, pf_c_range, 'g-', linewidth=2)
ax.axvline(c0, color='r', linestyle='--', linewidth=2, label=f'Nominal c = {c0}')
ax.axvline(c_star, color='orange', linestyle=':', linewidth=2, label=f'Design c* = {c_star:.1f}')
ax.set_xlabel('Damping (N·s/m)', fontsize=11)
ax.set_ylabel('Failure Probability', fontsize=11)
ax.set_title('Failure Probability vs Damping', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图3：不同刚度下的失效概率
ax = axes[1, 0]
k_range = np.linspace(k0*0.5, k0*2.0, 20)
pf_k_range = []
for k in k_range:
    pf, _, _ = monte_carlo_reliability(m0, c0, k, sigma_m, sigma_c, sigma_k, f, dt, threshold, n_samples=2000)
    pf_k_range.append(pf)

ax.plot(k_range, pf_k_range, 'orange', linewidth=2)
ax.axvline(k0, color='r', linestyle='--', linewidth=2, label=f'Nominal k = {k0}')
ax.axvline(k_star, color='purple', linestyle=':', linewidth=2, label=f'Design k* = {k_star:.1f}')
ax.set_xlabel('Stiffness (N/m)', fontsize=11)
ax.set_ylabel('Failure Probability', fontsize=11)
ax.set_title('Failure Probability vs Stiffness', fontsize=12)
ax.legend()
ax.grid(True, alpha=0.3)

# 子图4：三维参数空间中的失效边界
ax = axes[1, 1]
ax.axis('off')
summary_text = "FORM Analysis Summary:\n\n"
summary_text += f"Reliability Index: β = {beta_form:.4f}\n"
summary_text += f"Failure Probability: Pf = {pf_form:.6f}\n"
summary_text += f"Monte Carlo Pf: {pf_mc:.6f}\n\n"
summary_text += "Design Point (Physical Space):\n"
summary_text += f"  m* = {m_star:.2f} kg\n"
summary_text += f"  c* = {c_star:.2f} N·s/m\n"
summary_text += f"  k* = {k_star:.2f} N/m\n\n"
summary_text += "Importance Factors:\n"
summary_text += f"  Mass: {alpha[0]**2*100:.1f}%\n"
summary_text += f"  Damping: {alpha[1]**2*100:.1f}%\n"
summary_text += f"  Stiffness: {alpha[2]**2*100:.1f}%\n\n"
summary_text += "Most Critical Parameter:\n"
max_idx = np.argmax(np.abs(alpha))
param_names = ['Mass', 'Damping', 'Stiffness']
summary_text += f"  {param_names[max_idx]} (α = {alpha[max_idx]:.4f})"

ax.text(0.1, 0.5, summary_text, fontsize=11, verticalalignment='center',
        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.tight_layout()
plt.savefig('case4_sensitivity_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: case4_sensitivity_analysis.png")

# 创建动画：可靠指标的迭代收敛
print("\n正在生成FORM迭代动画...")

fig, ax = plt.subplots(figsize=(10, 8))
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_xlabel('u_m (Mass)', fontsize=12)
ax.set_ylabel('u_c (Damping)', fontsize=12)
ax.set_title('FORM Iteration Process', fontsize=14)
ax.grid(True, alpha=0.3)

# 绘制单位圆
phi = np.linspace(0, 2*np.pi, 100)
ax.plot(np.cos(phi), np.sin(phi), 'k--', linewidth=1, alpha=0.5, label='Unit Circle')

# 模拟迭代过程
n_iter_anim = 10
beta_history = np.linspace(0, beta_form, n_iter_anim)
u_history = np.array([[beta * alpha[0], beta * alpha[1]] for beta in beta_history])

line_path, = ax.plot([], [], 'b-', linewidth=1.5, alpha=0.5)
point_current, = ax.plot([], [], 'ro', markersize=10)
point_final, = ax.plot([], [], 'g*', markersize=15, label='Design Point')

beta_text = ax.text(0.02, 0.98, '', transform=ax.transAxes, fontsize=11,
                    verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

def init():
    line_path.set_data([], [])
    point_current.set_data([], [])
    point_final.set_data([u_star[0]], [u_star[1]])
    beta_text.set_text('')
    return [line_path, point_current, point_final, beta_text]

def update(frame):
    line_path.set_data(u_history[:frame+1, 0], u_history[:frame+1, 1])
    point_current.set_data([u_history[frame, 0]], [u_history[frame, 1]])
    beta_text.set_text(f'Iteration: {frame+1}\nβ = {beta_history[frame]:.4f}')
    return [line_path, point_current, point_final, beta_text]

anim = FuncAnimation(fig, update, frames=n_iter_anim, init_func=init, blit=True)
writer = PillowWriter(fps=2)
anim.save('case4_form_iteration.gif', writer=writer)
plt.close()
print("  已保存: case4_form_iteration.gif")

print("\n" + "=" * 60)
print("案例4分析完成!")
print("=" * 60)