第012篇：结构动力学中的不确定性分析

摘要

结构动力学中的不确定性分析是现代工程设计与安全评估的核心组成部分。在实际工程中，结构系统的材料属性、几何尺寸、边界条件以及外部荷载往往存在不可避免的随机性和认知不确定性。本教程系统地介绍结构动力学中不确定性分析的理论基础、数学方法和数值实现技术。首先阐述不确定性的分类与来源，包括偶然不确定性（Aleatory Uncertainty）和认知不确定性（Epistemic Uncertainty）。然后详细介绍概率分析方法，包括蒙特卡洛模拟、拉丁超立方抽样、响应面法等随机分析方法。接着探讨基于可靠性的设计优化方法，以及敏感性分析在结构动力学中的应用。通过Python实现多个典型案例，包括单自由度系统的随机响应分析、多自由度系统的参数敏感性研究、以及基于可靠度的结构优化设计，帮助读者掌握不确定性分析在实际工程问题中的应用技巧。

关键词

不确定性分析，蒙特卡洛模拟，拉丁超立方抽样，可靠性分析，敏感性分析，响应面法，结构可靠性，随机振动，概率设计

---

1. 不确定性理论基础

1.1 不确定性的分类与来源

在结构动力学分析中，不确定性是指我们对系统行为、参数或模型认知的不完全性。根据不确定性的本质，可以将其分为两大类：

偶然不确定性（Aleatory Uncertainty），也称为固有不确定性或随机不确定性，是指自然界中固有的、不可消除的随机性。这种不确定性来源于物理系统的内在变异性，即使我们拥有完美的知识和无限的数据，也无法消除这种不确定性。在结构动力学中，偶然不确定性的典型来源包括：

• 材料属性的随机性：混凝土强度、钢材屈服强度、弹性模量等力学性能参数存在批次差异和位置变异
• 几何尺寸的制造误差：构件截面尺寸、长度、厚度等存在加工公差
• 外部荷载的随机性：地震动、风荷载、波浪荷载等自然灾害具有显著的随机特征
• 环境条件的波动：温度、湿度、腐蚀环境等随时间和空间变化

认知不确定性（Epistemic Uncertainty），也称为知识不确定性或系统不确定性，是指由于知识不足、数据缺乏或模型简化导致的不确定性。这种不确定性原则上可以通过收集更多数据、改进测量技术或完善理论模型来减少。认知不确定性的主要来源包括：

• 模型形式的不确定性：数学模型对真实物理系统的近似程度
• 参数识别的不确定性：由于测试数据有限导致的参数估计误差
• 边界条件的不确定性：实际约束条件与理想化假设的差异
• 数值方法的不确定性：离散化误差、截断误差、收敛准则等

从数学角度，不确定性可以用概率论、模糊理论、区间分析或证据理论等不同框架来描述。概率论是最常用和最成熟的工具，适用于具有足够数据支持随机变量分布假设的情况。

1.2 概率论基础

概率论为不确定性量化提供了严格的数学框架。在结构动力学不确定性分析中，我们需要掌握以下核心概念：

随机变量是描述不确定量的基本工具。根据取值特点，随机变量可分为：

• 连续型随机变量：如材料强度、构件尺寸，用概率密度函数（PDF）描述
• 离散型随机变量：如缺陷数量、失效模式，用概率质量函数（PMF）描述

概率分布描述了随机变量的统计规律。工程中最常用的分布包括：

1. 正态分布（Normal/Gaussian Distribution）

   正态分布是最常见的连续概率分布，其概率密度函数为：

   $$f_X(x)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$

   其中，$\mu$ 为均值，$\sigma$ 为标准差。正态分布适用于描述许多自然现象和工程参数，如材料强度、测量误差等。

2. 对数正态分布（Lognormal Distribution）

   如果随机变量 $Y=\ln (X)$ 服从正态分布，则 $X$ 服从对数正态分布。其概率密度函数为：

   $$f_X(x)=\frac{1}{x{\sigma }_Y\sqrt{2\pi }}\exp \left(-\frac{(\ln x-{\mu }_Y{)}^2}{2{\sigma }_Y^2}\right),x>0$$

   对数正态分布适用于描述非负随机变量，如材料强度、疲劳寿命等。其特点是右偏分布，即出现大值的概率比正态分布更高。

3. 威布尔分布（Weibull Distribution）

   威布尔分布是可靠性工程中最常用的分布之一，其累积分布函数为：

   $$F_X(x)=1-\exp \left[-{\left(\frac{x-\gamma }{\eta }\right)}^{\beta }\right],x\ge \gamma$$

   其中，$\beta$ 为形状参数，$\eta$ 为尺度参数，$\gamma$ 为位置参数。威布尔分布的灵活性使其能够描述各种类型的失效数据。

4. 均匀分布（Uniform Distribution）

   当只知道随机变量的取值范围而缺乏其他信息时，常采用均匀分布：

   $$f_X(x)=\{\begin{array}{ll}\frac{1}{b-a},& a\le x\le b\\ 0,& \text{其他}\end{array}$$

统计矩是描述随机变量分布特征的重要指标：

• 均值（一阶原点矩）：${\mu }_X=E[X]={\int }_{-\infty }^{\infty }x{f}_X(x)dx$
• 方差（二阶中心矩）：${\sigma }_X^2=E[(X-{\mu }_X{)}^2]={\int }_{-\infty }^{\infty }(x-{\mu }_X{)}^2{f}_X(x)dx$
• 变异系数：${\delta }_X=\frac{{\sigma }_X}{{\mu }_X}$，表示相对离散程度
• 偏度（三阶标准化矩）：${\gamma }_1=\frac{E[(X-{\mu }_X{)}^3]}{{\sigma }_X^3}$，描述分布的不对称性
• 峰度（四阶标准化矩）：${\gamma }_2=\frac{E[(X-{\mu }_X{)}^4]}{{\sigma }_X^4}$，描述分布的尖锐程度

1.3 结构可靠性基本概念

结构可靠性是指结构在规定的时间内、规定的条件下完成预定功能的概率。与可靠性相对的概念是失效概率，两者满足：

$$P_s+P_f=1$$

其中，$P_s$ 为可靠度，$P_f$ 为失效概率。

在结构可靠性分析中，通常定义极限状态函数（Limit State Function）来区分安全域和失效域：

$$g(X)=R-S$$

其中，$R$ 为结构抗力（Resistance），$S$ 为荷载效应（Load Effect），$X$ 为随机变量向量。

• 当 $g(X)>0$ 时，结构处于安全状态
• 当 $g(X)=0$ 时，结构处于极限状态
• 当 $g(X)<0$ 时，结构处于失效状态

失效概率可以表示为：

$$P_f=P(g(X)<0)={\int }_{g(X)<0}{f}_X(x)dx$$

由于直接计算上述多维积分通常非常困难，工程实践中发展了多种近似和数值方法。

可靠指标（Reliability Index）$\beta$ 是衡量结构安全性的重要指标，定义为标准正态空间中原点到极限状态曲面的最短距离。对于线性极限状态函数和正态随机变量，可靠指标与失效概率的关系为：

$$P_f=\Phi (-\beta )$$

$$\beta =-{\Phi }^{-1}(P_f)$$

其中，$\Phi (\cdot )$ 为标准正态累积分布函数。

1.4 不确定性传播

在结构动力学中，输入参数的随机性会通过系统传递，导致响应也具有随机性，这就是不确定性传播问题。对于单自由度系统：

$$m\ddot{X}+c\dot{X}+kX=F(t)$$

如果质量 $m$、阻尼 $c$、刚度 $k$ 或外力 $F(t)$ 是随机变量，则位移响应 $X(t)$ 也是随机过程。

不确定性传播的数学描述包括：

矩传播法：通过泰勒展开近似计算响应的统计矩

$$E[Y]\approx g({\mu }_X)+\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n\frac{{\partial }^{2}g}{\partial X_i\partial X_j}{|}_{{\mu }_X}\text{Cov}(X_i,X_j)$$

$$\text{Var}[Y]\approx \sum _{i=1}^n\sum _{j=1}^n\frac{\partial g}{\partial X_i}{|}_{{\mu }_X}\frac{\partial g}{\partial X_j}{|}_{{\mu }_X}\text{Cov}(X_i,X_j)$$

概率密度演化法：通过求解Fokker-Planck-Kolmogorov方程获得响应的概率密度函数

$$\frac{\partial p_{X\dot{X}}}{\partial t}=-\dot{x}\frac{\partial p_{X\dot{X}}}{\partial x}+\frac{\partial }{\partial \dot{x}}\left[\left(\frac{c}{m}\dot{x}+\frac{k}{m}x-\frac{f(t)}{m}\right)p_{X\dot{X}}\right]$$

蒙特卡洛模拟：通过大量随机抽样和确定性分析获得响应的统计特性

---

2. 随机抽样方法

2.1 蒙特卡洛模拟方法

蒙特卡洛模拟（Monte Carlo Simulation, MCS）是一种基于随机抽样的数值方法，通过大量重复计算来估计复杂系统的统计特性。其基本思想是：利用大数定律，当样本量足够大时，样本均值收敛于期望值。

蒙特卡洛方法的基本步骤：

1. 定义输入随机变量：确定所有不确定参数的概率分布
2. 生成随机样本：根据概率分布生成大量随机数
3. 执行确定性分析：对每个样本进行结构动力学分析
4. 统计分析：计算响应的统计矩、概率分布或失效概率

蒙特卡洛估计的收敛性：

设 $N$ 为样本数，${\hat{\mu }}_Y$ 为响应均值的估计值，则估计误差的标准差为：

$${\sigma }_{{\hat{\mu }}_Y}=\frac{{\sigma }_Y}{\sqrt{N}}$$

这意味着要将估计误差减半，需要将样本量增加4倍。蒙特卡洛方法的收敛速度与问题维度无关，这是其最大优势，但对于小概率事件（如 $P_f=10^{-6}$），需要极大的样本量才能获得可靠的估计。

方差缩减技术是提高蒙特卡洛效率的重要手段，包括：

• 重要性抽样（Importance Sampling）：在重要区域增加抽样密度
• 分层抽样（Stratified Sampling）：将样本空间划分为若干层，每层独立抽样
• 控制变量法（Control Variates）：利用已知解析解的辅助变量减少方差
• 对偶变量法（Antithetic Variates）：利用负相关样本对减少方差

2.2 拉丁超立方抽样

拉丁超立方抽样（Latin Hypercube Sampling, LHS）是一种分层抽样技术，能够在保持随机性的同时确保样本在整个参数空间均匀分布。

LHS的基本原理：

对于 $n$ 维随机变量，将每个维度划分为 $N$ 个等概率区间，然后在每个区间中随机抽取一个样本，并确保在每个维度上，所有区间恰好有一个样本。

LHS的优势：

1. 更好的空间填充性：相比简单随机抽样，LHS能够更均匀地覆盖参数空间
2. 更快的收敛速度：对于光滑函数，LHS的收敛速度通常为 $O(N^{-1})$，优于简单随机抽样的 $O(N^{-1/2})$
3. 适合小样本分析：在样本量有限时，LHS能提供更稳定的估计

LHS的数学描述：

设 $X$ 为 $n$ 维随机变量，$F_i(x_i)$ 为第 $i$ 个分量的累积分布函数。LHS样本的生成过程为：

1. 生成 $N\times n$ 的随机排列矩阵 $P$，每列是 $(1,2,...,N)$ 的随机排列
2. 生成 $N\times n$ 的均匀随机数矩阵 $U$，元素在 $[0,1)$ 上均匀分布
3. 计算概率值：$P_{ij}=\frac{P_{ij}-U_{ij}}{N}$
4. 通过逆变换得到样本：$X_{ij}=F_i^{-1}(P_{ij})$

2.3 其他高级抽样技术

Sobol序列是一种低差异序列（Low-Discrepancy Sequence），也称为准随机序列。与伪随机数不同，低差异序列在空间中分布更加均匀，能够更有效地填充参数空间。

高斯过程抽样利用高斯过程回归模型，在已有样本的基础上，通过贝叶斯更新指导新样本的选择，实现自适应抽样。

子集模拟（Subset Simulation）是处理小失效概率问题的高效方法，通过引入中间失效事件，将 rare event 的估计转化为一系列条件概率的乘积。

---

3. 敏感性分析方法

3.1 局部敏感性分析

局部敏感性分析研究输入参数在基准值附近微小变化时，对输出响应的影响程度。最常用的方法是偏导数法。

对于模型 $Y=g(X_1,X_2,...,X_n)$，参数 $X_i$ 的局部敏感性系数定义为：

$$S_i=\frac{\partial g}{\partial X_i}{|}_{X={\mu }_X}$$

为了消除量纲影响，常采用标准化敏感性系数：

$$S_i^{\text{std}}=\frac{\partial g}{\partial X_i}{|}_{X={\mu }_X}\cdot \frac{{\sigma }_{X_i}}{{\sigma }_Y}$$

标准化敏感性系数表示当参数 $X_i$ 变化一个标准差时，响应 $Y$ 变化的标准差数量。

有限差分法是计算敏感性系数的数值方法：

$$\frac{\partial g}{\partial X_i}\approx \frac{g(X_i+\Delta X_i)-g(X_i)}{\Delta X_i}$$

前向差分、后向差分和中心差分是常用的差分格式，其中中心差分具有二阶精度：

$$\frac{\partial g}{\partial X_i}\approx \frac{g(X_i+\Delta X_i)-g(X_i-\Delta X_i)}{2\Delta X_i}$$

3.2 全局敏感性分析

全局敏感性分析考虑参数在其整个取值范围内的变化对响应的影响，能够识别参数间的交互作用。

Sobol敏感性指数是基于方差分解的全局敏感性指标。响应的总方差可以分解为：

$$\text{Var}(Y)=\sum _i{V}_i+\sum _{i<j}{V}_{ij}+\sum _{i<j<k}{V}_{ijk}+...+V_{\mathrm{12...}n}$$

其中：

• 一阶效应：$V_i={\text{Var}}_{X_i}(E_{X_{\sim i}}[Y|X_i])$
• 二阶交互效应：$V_{ij}={\text{Var}}_{X_i,X_j}(E_{X_{\sim ij}}[Y|X_i,X_j])-V_i-V_j$

一阶Sobol指数衡量单个参数对响应方差的贡献：

$$S_i=\frac{V_i}{\text{Var}(Y)}$$

总效应指数衡量参数及其与其他参数交互作用的总贡献：

$$S_{Ti}=1-\frac{{\text{Var}}_{X_{\sim i}}(E_{X_i}[Y|X_{\sim i}])}{\text{Var}(Y)}$$

Morris方法是一种基于"一次改变一个因素"（OAT）思想的筛选方法，通过计算基本效应（Elementary Effects）来识别重要参数。

3.3 基于回归的敏感性分析

当样本量足够大时，可以通过建立响应与输入参数之间的回归模型来分析敏感性。

线性回归模型：

$$Y={\beta }_0+\sum _{i=1}^n{\beta }_i{X}_i+ϵ$$

标准化回归系数（SRC）可以作为敏感性指标：

$${\text{SRC}}_i={\beta }_i\cdot \frac{{\sigma }_{X_i}}{{\sigma }_Y}$$

偏相关系数（Partial Correlation Coefficient, PCC）衡量在控制其他参数影响后，单个参数与响应的相关性。

秩相关系数（Spearman’s Rank Correlation Coefficient）对单调非线性关系敏感，适用于非线性模型。

---

4. 结构可靠性分析方法

4.1 一次二阶矩法（FORM）

一次二阶矩法（First Order Reliability Method, FORM）是结构可靠性分析中最常用的近似方法，通过将极限状态函数在验算点处线性化来估计失效概率。

FORM的基本步骤：

1. 标准化变换：将相关非正态随机变量变换为标准正态空间中的独立变量

   对于正态变量：$U_i=\frac{X_i-{\mu }_{X_i}}{{\sigma }_{X_i}}$

   对于非正态变量，使用Rackwitz-Fiessler变换或Nataf变换

2. 寻找验算点（Design Point）：在标准正态空间中找到极限状态曲面上距离原点最近的点

   验算点 $u^{\ast }$ 是以下约束优化问题的解：

   $$\min \Vert u\Vert \text{s.t.}g(u)=0$$

3. 计算可靠指标：$\beta =\Vert u^{\ast }\Vert$

4. 估计失效概率：$P_f\approx \Phi (-\beta )$

HL-RF算法是求解验算点的迭代算法：

$$u^{(k+1)}=\frac{1}{\Vert \nabla g(u^{(k)}){\Vert }^2}\left[\nabla g(u^{(k)}{)}^T{u}^{(k)}-g(u^{(k)})\right]\nabla g(u^{(k)})$$

4.2 二次二阶矩法（SORM）

二次二阶矩法（Second Order Reliability Method, SORM）在验算点处使用二次曲面近似极限状态函数，能够获得更精确的失效概率估计。

SORM的失效概率估计为：

$$P_f^{\text{SORM}}\approx \Phi (-\beta )\prod _{i=1}^{n-1}(1-\beta {\kappa }_i{)}^{-1/2}$$

其中，${\kappa }_i$ 为极限状态曲面在验算点处的主曲率。

4.3 响应面法

响应面法（Response Surface Method, RSM）通过构造一个简单的显式函数（通常是多项式）来近似复杂的隐式极限状态函数。

二次响应面模型：

$$\hat{g}(X)=a_0+\sum _{i=1}^n{a}_i{X}_i+\sum _{i=1}^n\sum _{j=i}^n{a}_{ij}{X}_i{X}_j$$

响应面系数的确定通常采用最小二乘法，通过在设计点附近选择适当的样本点来拟合响应面。

自适应响应面法通过迭代更新验算点和响应面，逐步逼近真实的极限状态函数。

---

5. Python实现：不确定性分析案例

5.1 案例1：单自由度系统的蒙特卡洛分析

本案例通过蒙特卡洛模拟分析单自由度系统在随机参数下的响应统计特性。

# -*- coding: utf-8 -*-
"""
主题012：结构动力学中的不确定性分析
案例1：单自由度系统的蒙特卡洛分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.integrate import solve_ivp
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\结构动力学仿真\主题012'
os.makedirs(output_dir, exist_ok=True)

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

print("=" * 70)
print("主题012：结构动力学中的不确定性分析")
print("=" * 70)

# ============================================================
# 案例1：单自由度系统的蒙特卡洛分析
# ============================================================
print("\n案例1：单自由度系统的蒙特卡洛分析")
print("-" * 70)

# 定义系统参数（随机变量）
# 质量 m ~ N(1.0, 0.05) kg
# 刚度 k ~ Lognormal(ln(100), 0.1) N/m
# 阻尼比 zeta ~ Uniform(0.02, 0.08)

m_mean, m_std = 1.0, 0.05
k_mean, k_cov = 100.0, 0.1  # 对数正态分布参数
zeta_low, zeta_high = 0.02, 0.08

# 转换为对数正态分布参数
k_ln_mean = np.log(k_mean / np.sqrt(1 + k_cov**2))
k_ln_std = np.sqrt(np.log(1 + k_cov**2))

print("随机参数定义：")
print(f"  质量 m ~ N({m_mean}, {m_std}²) kg")
print(f"  刚度 k ~ Lognormal({k_ln_mean:.3f}, {k_ln_std:.3f}²) N/m")
print(f"  阻尼比 ζ ~ Uniform({zeta_low}, {zeta_high})")

# 蒙特卡洛模拟参数
n_samples = 5000
np.random.seed(42)

# 生成随机样本
m_samples = np.random.normal(m_mean, m_std, n_samples)
k_samples = np.random.lognormal(k_ln_mean, k_ln_std, n_samples)
zeta_samples = np.random.uniform(zeta_low, zeta_high, n_samples)

print(f"\n蒙特卡洛模拟参数：")
print(f"  样本数量: {n_samples}")

# 计算固有频率和周期
omega_n_samples = np.sqrt(k_samples / m_samples)
f_n_samples = omega_n_samples / (2 * np.pi)
T_n_samples = 1.0 / f_n_samples

# 计算阻尼系数
c_samples = 2 * zeta_samples * np.sqrt(k_samples * m_samples)

print("\n正在计算系统响应统计特性...")

# 定义简谐激励下的稳态响应振幅
def steady_state_amplitude(m, c, k, p0, omega):
    """计算稳态响应振幅"""
    omega_n = np.sqrt(k / m)
    zeta = c / (2 * np.sqrt(m * k))
    beta = omega / omega_n
    
    H = 1.0 / np.sqrt((1 - beta**2)**2 + (2 * zeta * beta)**2)
    return (p0 / k) * H

# 激励参数
p0 = 50.0  # 激励幅值
omega = 8.0  # 激励频率

# 计算响应振幅
A_samples = np.array([steady_state_amplitude(m, c, k, p0, omega) 
                      for m, c, k in zip(m_samples, c_samples, k_samples)])

# 计算动力放大系数
DAF_samples = A_samples * k_samples / p0

print("  响应统计特性计算完成")

# 绘制参数分布
fig1, axes = plt.subplots(2, 3, figsize=(15, 10))

# 质量分布
axes[0, 0].hist(m_samples, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
x_m = np.linspace(m_mean - 4*m_std, m_mean + 4*m_std, 100)
axes[0, 0].plot(x_m, stats.norm.pdf(x_m, m_mean, m_std), 'r-', linewidth=2, label='PDF')
axes[0, 0].axvline(m_mean, color='g', linestyle='--', linewidth=2, label=f'Mean={m_mean:.3f}')
axes[0, 0].set_xlabel('Mass (kg)', fontsize=11)
axes[0, 0].set_ylabel('PDF', fontsize=11)
axes[0, 0].set_title('Mass Distribution', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 刚度分布
axes[0, 1].hist(k_samples, bins=50, density=True, alpha=0.7, color='green', edgecolor='black')
x_k = np.linspace(60, 150, 100)
axes[0, 1].plot(x_k, stats.lognorm.pdf(x_k, k_ln_std, scale=np.exp(k_ln_mean)), 'r-', linewidth=2)
axes[0, 1].axvline(k_mean, color='g', linestyle='--', linewidth=2, label=f'Mean={k_mean:.1f}')
axes[0, 1].set_xlabel('Stiffness (N/m)', fontsize=11)
axes[0, 1].set_ylabel('PDF', fontsize=11)
axes[0, 1].set_title('Stiffness Distribution', fontsize=12, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 阻尼比分布
axes[0, 2].hist(zeta_samples, bins=50, density=True, alpha=0.7, color='orange', edgecolor='black')
axes[0, 2].axvline((zeta_low + zeta_high)/2, color='g', linestyle='--', linewidth=2, 
                   label=f'Mean={(zeta_low+zeta_high)/2:.3f}')
axes[0, 2].set_xlabel('Damping Ratio', fontsize=11)
axes[0, 2].set_ylabel('PDF', fontsize=11)
axes[0, 2].set_title('Damping Ratio Distribution', fontsize=12, fontweight='bold')
axes[0, 2].legend()
axes[0, 2].grid(True, alpha=0.3)

# 固有频率分布
axes[1, 0].hist(omega_n_samples, bins=50, density=True, alpha=0.7, color='purple', edgecolor='black')
axes[1, 0].axvline(np.mean(omega_n_samples), color='r', linestyle='--', linewidth=2, 
                   label=f'Mean={np.mean(omega_n_samples):.2f}')
axes[1, 0].set_xlabel('Natural Frequency (rad/s)', fontsize=11)
axes[1, 0].set_ylabel('PDF', fontsize=11)
axes[1, 0].set_title('Natural Frequency Distribution', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 响应振幅分布
axes[1, 1].hist(A_samples, bins=50, density=True, alpha=0.7, color='red', edgecolor='black')
axes[1, 1].axvline(np.mean(A_samples), color='b', linestyle='--', linewidth=2, 
                   label=f'Mean={np.mean(A_samples):.4f}')
axes[1, 1].set_xlabel('Response Amplitude (m)', fontsize=11)
axes[1, 1].set_ylabel('PDF', fontsize=11)
axes[1, 1].set_title('Response Amplitude Distribution', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

# 动力放大系数分布
axes[1, 2].hist(DAF_samples, bins=50, density=True, alpha=0.7, color='cyan', edgecolor='black')
axes[1, 2].axvline(np.mean(DAF_samples), color='r', linestyle='--', linewidth=2, 
                   label=f'Mean={np.mean(DAF_samples):.2f}')
axes[1, 2].set_xlabel('Dynamic Amplification Factor', fontsize=11)
axes[1, 2].set_ylabel('PDF', fontsize=11)
axes[1, 2].set_title('DAF Distribution', fontsize=12, fontweight='bold')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/mc_parameter_distributions.png', dpi=150, bbox_inches='tight')
plt.close()
print("  参数分布图已保存")

# 统计结果汇总
print("\n蒙特卡洛分析结果：")
print("-" * 50)
print(f"{'Parameter':<25} {'Mean':<12} {'Std':<12} {'COV (%)':<10}")
print("-" * 50)
print(f"{'Mass (kg)':<25} {np.mean(m_samples):<12.4f} {np.std(m_samples):<12.4f} {100*np.std(m_samples)/np.mean(m_samples):<10.2f}")
print(f"{'Stiffness (N/m)':<25} {np.mean(k_samples):<12.2f} {np.std(k_samples):<12.2f} {100*np.std(k_samples)/np.mean(k_samples):<10.2f}")
print(f"{'Damping Ratio':<25} {np.mean(zeta_samples):<12.4f} {np.std(zeta_samples):<12.4f} {100*np.std(zeta_samples)/np.mean(zeta_samples):<10.2f}")
print(f"{'Natural Freq (rad/s)':<25} {np.mean(omega_n_samples):<12.2f} {np.std(omega_n_samples):<12.2f} {100*np.std(omega_n_samples)/np.mean(omega_n_samples):<10.2f}")
print(f"{'Response Amp (m)':<25} {np.mean(A_samples):<12.4f} {np.std(A_samples):<12.4f} {100*np.std(A_samples)/np.mean(A_samples):<10.2f}")
print(f"{'DAF':<25} {np.mean(DAF_samples):<12.2f} {np.std(DAF_samples):<12.2f} {100*np.std(DAF_samples)/np.mean(DAF_samples):<10.2f}")
print("-" * 50)

# 绘制散点图矩阵
fig2, axes = plt.subplots(3, 3, figsize=(14, 14))

params = [m_samples, k_samples, zeta_samples]
param_names = ['Mass (kg)', 'Stiffness (N/m)', 'Damping Ratio']

for i in range(3):
    for j in range(3):
        if i == j:
            # 对角线：直方图
            axes[i, j].hist(params[i], bins=40, alpha=0.7, color='steelblue', edgecolor='black')
            axes[i, j].set_title(param_names[i], fontsize=10, fontweight='bold')
        else:
            # 非对角线：散点图
            axes[i, j].scatter(params[j], params[i], alpha=0.3, s=1, c='blue')
            axes[i, j].set_xlabel(param_names[j], fontsize=9)
            axes[i, j].set_ylabel(param_names[i], fontsize=9)
            # 计算相关系数
            corr = np.corrcoef(params[j], params[i])[0, 1]
            axes[i, j].text(0.05, 0.95, f'r={corr:.3f}', transform=axes[i, j].transAxes,
                           fontsize=9, verticalalignment='top',
                           bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
        axes[i, j].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/mc_scatter_matrix.png', dpi=150, bbox_inches='tight')
plt.close()
print("  散点图矩阵已保存")

# 收敛性分析
print("\n正在进行收敛性分析...")
sample_sizes = np.logspace(2, np.log10(n_samples), 20, dtype=int)
sample_sizes = np.unique(sample_sizes)

mean_convergence = []
std_convergence = []

for n in sample_sizes:
    mean_convergence.append(np.mean(A_samples[:n]))
    std_convergence.append(np.std(A_samples[:n]))

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

axes[0].semilogx(sample_sizes, mean_convergence, 'b-o', linewidth=2, markersize=6)
axes[0].axhline(np.mean(A_samples), color='r', linestyle='--', linewidth=2, label='Final Mean')
axes[0].set_xlabel('Sample Size', fontsize=11)
axes[0].set_ylabel('Mean of Response Amplitude', fontsize=11)
axes[0].set_title('Convergence of Mean', fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].semilogx(sample_sizes, std_convergence, 'g-s', linewidth=2, markersize=6)
axes[1].axhline(np.std(A_samples), color='r', linestyle='--', linewidth=2, label='Final Std')
axes[1].set_xlabel('Sample Size', fontsize=11)
axes[1].set_ylabel('Std of Response Amplitude', fontsize=11)
axes[1].set_title('Convergence of Standard Deviation', fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/mc_convergence.png', dpi=150, bbox_inches='tight')
plt.close()
print("  收敛性分析图已保存")

5.2 案例2：拉丁超立方抽样与对比分析

本案例比较蒙特卡洛模拟和拉丁超立方抽样的效率和精度。

# ============================================================
# 案例2：拉丁超立方抽样与对比分析
# ============================================================
print("\n案例2：拉丁超立方抽样与对比分析")
print("-" * 70)

def latin_hypercube_sampling(n_samples, distributions):
    """
    拉丁超立方抽样
    distributions: 列表，每个元素为(scipy.stats分布对象, 参数元组)
    """
    n_vars = len(distributions)
    samples = np.zeros((n_samples, n_vars))
    
    for i, (dist, params) in enumerate(distributions):
        # 生成分层
        strata = np.random.permutation(n_samples)
        # 在每个分层内随机抽样
        u = (strata + np.random.rand(n_samples)) / n_samples
        # 逆变换抽样
        samples[:, i] = dist.ppf(u, *params)
    
    return samples

# 定义分布
# 正态分布: (stats.norm, (mean, std))
# 对数正态分布: (stats.lognorm, (s, scale))
# 均匀分布: (stats.uniform, (low, high-low))
distributions = [
    (stats.norm, (m_mean, m_std)),
    (stats.lognorm, (k_ln_std, 0, np.exp(k_ln_mean))),
    (stats.uniform, (zeta_low, zeta_high - zeta_low))
]

# 生成LHS样本
lhs_samples = latin_hypercube_sampling(n_samples, distributions)
m_lhs = lhs_samples[:, 0]
k_lhs = lhs_samples[:, 1]
zeta_lhs = lhs_samples[:, 2]

# 计算响应
omega_n_lhs = np.sqrt(k_lhs / m_lhs)
c_lhs = 2 * zeta_lhs * np.sqrt(k_lhs * m_lhs)
A_lhs = np.array([steady_state_amplitude(m, c, k, p0, omega) 
                  for m, c, k in zip(m_lhs, c_lhs, k_lhs)])

print("拉丁超立方抽样完成")

# 对比分析
fig4, axes = plt.subplots(2, 3, figsize=(15, 10))

# 质量对比
axes[0, 0].hist(m_samples, bins=50, density=True, alpha=0.5, label='MCS', color='blue')
axes[0, 0].hist(m_lhs, bins=50, density=True, alpha=0.5, label='LHS', color='red')
axes[0, 0].set_xlabel('Mass (kg)', fontsize=11)
axes[0, 0].set_ylabel('PDF', fontsize=11)
axes[0, 0].set_title('Mass: MCS vs LHS', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 刚度对比
axes[0, 1].hist(k_samples, bins=50, density=True, alpha=0.5, label='MCS', color='blue')
axes[0, 1].hist(k_lhs, bins=50, density=True, alpha=0.5, label='LHS', color='red')
axes[0, 1].set_xlabel('Stiffness (N/m)', fontsize=11)
axes[0, 1].set_ylabel('PDF', fontsize=11)
axes[0, 1].set_title('Stiffness: MCS vs LHS', fontsize=12, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 阻尼比对比
axes[0, 2].hist(zeta_samples, bins=50, density=True, alpha=0.5, label='MCS', color='blue')
axes[0, 2].hist(zeta_lhs, bins=50, density=True, alpha=0.5, label='LHS', color='red')
axes[0, 2].set_xlabel('Damping Ratio', fontsize=11)
axes[0, 2].set_ylabel('PDF', fontsize=11)
axes[0, 2].set_title('Damping Ratio: MCS vs LHS', fontsize=12, fontweight='bold')
axes[0, 2].legend()
axes[0, 2].grid(True, alpha=0.3)

# 响应振幅对比
axes[1, 0].hist(A_samples, bins=50, density=True, alpha=0.5, label='MCS', color='blue')
axes[1, 0].hist(A_lhs, bins=50, density=True, alpha=0.5, label='LHS', color='red')
axes[1, 0].set_xlabel('Response Amplitude (m)', fontsize=11)
axes[1, 0].set_ylabel('PDF', fontsize=11)
axes[1, 0].set_title('Response Amplitude: MCS vs LHS', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 二维散点对比：质量-刚度
axes[1, 1].scatter(m_samples[::10], k_samples[::10], alpha=0.3, s=5, label='MCS', c='blue')
axes[1, 1].scatter(m_lhs[::10], k_lhs[::10], alpha=0.3, s=5, label='LHS', c='red')
axes[1, 1].set_xlabel('Mass (kg)', fontsize=11)
axes[1, 1].set_ylabel('Stiffness (N/m)', fontsize=11)
axes[1, 1].set_title('Mass vs Stiffness', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

# 二维散点对比：阻尼比-响应
axes[1, 2].scatter(zeta_samples[::10], A_samples[::10], alpha=0.3, s=5, label='MCS', c='blue')
axes[1, 2].scatter(zeta_lhs[::10], A_lhs[::10], alpha=0.3, s=5, label='LHS', c='red')
axes[1, 2].set_xlabel('Damping Ratio', fontsize=11)
axes[1, 2].set_ylabel('Response Amplitude (m)', fontsize=11)
axes[1, 2].set_title('Damping Ratio vs Response', fontsize=12, fontweight='bold')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/lhs_comparison.png', dpi=150, bbox_inches='tight')
plt.close()
print("  LHS对比图已保存")

# 统计对比
print("\n抽样方法对比：")
print("-" * 70)
print(f"{'Method':<10} {'Response Mean':<15} {'Response Std':<15} {'95% CI Width':<15}")
print("-" * 70)

# MCS统计
mcs_mean = np.mean(A_samples)
mcs_std = np.std(A_samples)
mcs_ci = 1.96 * mcs_std / np.sqrt(n_samples)
print(f"{'MCS':<10} {mcs_mean:<15.6f} {mcs_std:<15.6f} {2*mcs_ci:<15.6f}")

# LHS统计
lhs_mean = np.mean(A_lhs)
lhs_std = np.std(A_lhs)
lhs_ci = 1.96 * lhs_std / np.sqrt(n_samples)
print(f"{'LHS':<10} {lhs_mean:<15.6f} {lhs_std:<15.6f} {2*lhs_ci:<15.6f}")
print("-" * 70)

5.3 案例3：敏感性分析

本案例通过多种方法分析各参数对系统响应的敏感性。

# ============================================================
# 案例3：敏感性分析
# ============================================================
print("\n案例3：敏感性分析")
print("-" * 70)

# 基于蒙特卡洛样本的敏感性分析
print("\n基于蒙特卡洛样本的敏感性分析：")

# 计算相关系数
corr_m = np.corrcoef(m_samples, A_samples)[0, 1]
corr_k = np.corrcoef(k_samples, A_samples)[0, 1]
corr_zeta = np.corrcoef(zeta_samples, A_samples)[0, 1]

print(f"  质量-响应相关系数: {corr_m:.4f}")
print(f"  刚度-响应相关系数: {corr_k:.4f}")
print(f"  阻尼比-响应相关系数: {corr_zeta:.4f}")

# 计算偏相关系数（简化计算）
from scipy.stats import pearsonr

# 多元线性回归分析
X = np.column_stack([m_samples, k_samples, zeta_samples])
X_centered = X - np.mean(X, axis=0)
A_centered = A_samples - np.mean(A_samples)

# 正规方程求解回归系数
beta = np.linalg.lstsq(X_centered, A_centered, rcond=None)[0]

print(f"\n多元线性回归系数：")
print(f"  β_m = {beta[0]:.6f}")
print(f"  β_k = {beta[1]:.6f}")
print(f"  β_ζ = {beta[2]:.6f}")

# 标准化回归系数（SRC）
beta_std = beta * np.std(X, axis=0) / np.std(A_samples)
print(f"\n标准化回归系数（SRC）：")
print(f"  SRC_m = {beta_std[0]:.4f}")
print(f"  SRC_k = {beta_std[1]:.4f}")
print(f"  SRC_ζ = {beta_std[2]:.4f}")

# 绘制敏感性分析结果
fig5, axes = plt.subplots(2, 2, figsize=(14, 12))

# 相关系数条形图
params = ['Mass', 'Stiffness', 'Damping']
correlations = [corr_m, corr_k, corr_zeta]
colors = ['blue', 'green', 'orange']

axes[0, 0].barh(params, correlations, color=colors, alpha=0.7, edgecolor='black')
axes[0, 0].axvline(0, color='black', linewidth=0.8)
axes[0, 0].set_xlabel('Correlation Coefficient', fontsize=11)
axes[0, 0].set_title('Correlation with Response', fontsize=12, fontweight='bold')
axes[0, 0].grid(True, alpha=0.3, axis='x')
for i, v in enumerate(correlations):
    axes[0, 0].text(v, i, f' {v:.3f}', va='center', fontsize=10)

# 标准化回归系数条形图
axes[0, 1].barh(params, beta_std, color=colors, alpha=0.7, edgecolor='black')
axes[0, 1].axvline(0, color='black', linewidth=0.8)
axes[0, 1].set_xlabel('Standardized Regression Coefficient', fontsize=11)
axes[0, 1].set_title('Standardized Regression Coefficients', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3, axis='x')
for i, v in enumerate(beta_std):
    axes[0, 1].text(v, i, f' {v:.3f}', va='center', fontsize=10)

# 散点图：质量 vs 响应
axes[1, 0].scatter(m_samples, A_samples, alpha=0.3, s=1, c='blue')
z_m = np.polyfit(m_samples, A_samples, 1)
p_m = np.poly1d(z_m)
axes[1, 0].plot(sorted(m_samples), p_m(sorted(m_samples)), 'r-', linewidth=2, 
                label=f'Linear fit (r={corr_m:.3f})')
axes[1, 0].set_xlabel('Mass (kg)', fontsize=11)
axes[1, 0].set_ylabel('Response Amplitude (m)', fontsize=11)
axes[1, 0].set_title('Mass vs Response', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 散点图：刚度 vs 响应
axes[1, 1].scatter(k_samples, A_samples, alpha=0.3, s=1, c='green')
z_k = np.polyfit(k_samples, A_samples, 1)
p_k = np.poly1d(z_k)
axes[1, 1].plot(sorted(k_samples), p_k(sorted(k_samples)), 'r-', linewidth=2, 
                label=f'Linear fit (r={corr_k:.3f})')
axes[1, 1].set_xlabel('Stiffness (N/m)', fontsize=11)
axes[1, 1].set_ylabel('Response Amplitude (m)', fontsize=11)
axes[1, 1].set_title('Stiffness vs Response', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/sensitivity_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  敏感性分析图已保存")

# 局部敏感性分析（基于有限差分）
print("\n局部敏感性分析（有限差分法）：")

# 基准值
m_base = 1.0
k_base = 100.0
zeta_base = 0.05
c_base = 2 * zeta_base * np.sqrt(k_base * m_base)

# 计算基准响应
A_base = steady_state_amplitude(m_base, c_base, k_base, p0, omega)

# 扰动
h = 0.01

# 质量敏感性
m_pert = m_base * (1 + h)
c_pert = 2 * zeta_base * np.sqrt(k_base * m_pert)
A_m_pert = steady_state_amplitude(m_pert, c_pert, k_base, p0, omega)
S_m = (A_m_pert - A_base) / (m_pert - m_base)

# 刚度敏感性
k_pert = k_base * (1 + h)
c_pert = 2 * zeta_base * np.sqrt(k_pert * m_base)
A_k_pert = steady_state_amplitude(m_base, c_pert, k_pert, p0, omega)
S_k = (A_k_pert - A_base) / (k_pert - k_base)

# 阻尼比敏感性
zeta_pert = zeta_base * (1 + h)
c_pert = 2 * zeta_pert * np.sqrt(k_base * m_base)
A_zeta_pert = steady_state_amplitude(m_base, c_pert, k_base, p0, omega)
S_zeta = (A_zeta_pert - A_base) / (zeta_pert - zeta_base)

print(f"  ∂A/∂m = {S_m:.6f}")
print(f"  ∂A/∂k = {S_k:.6f}")
print(f"  ∂A/∂ζ = {S_zeta:.6f}")

# 标准化局部敏感性系数
S_m_std = S_m * m_base / A_base
S_k_std = S_k * k_base / A_base
S_zeta_std = S_zeta * zeta_base / A_base

print(f"\n标准化局部敏感性系数：")
print(f"  (∂A/∂m)·(m/A) = {S_m_std:.4f}")
print(f"  (∂A/∂k)·(k/A) = {S_k_std:.4f}")
print(f"  (∂A/∂ζ)·(ζ/A) = {S_zeta_std:.4f}")

5.4 案例4：结构可靠性分析

本案例通过FORM方法分析结构的可靠性和失效概率。

# ============================================================
# 案例4：结构可靠性分析
# ============================================================
print("\n案例4：结构可靠性分析")
print("-" * 70)

# 定义极限状态函数
# g(R, S) = R - S，其中R为抗力，S为荷载效应
# 当 g < 0 时失效

# 假设：
# 抗力 R ~ Lognormal(μ_R, σ_R)
# 荷载效应 S ~ Normal(μ_S, σ_S)

mu_R = 200.0  # 抗力均值 (kN)
cov_R = 0.15  # 抗力变异系数
sigma_R = mu_R * cov_R

mu_S = 150.0  # 荷载效应均值 (kN)
cov_S = 0.20  # 荷载效应变异系数
sigma_S = mu_S * cov_S

print("可靠性分析参数：")
print(f"  抗力 R ~ Lognormal(μ={mu_R}, COV={cov_R})")
print(f"  荷载 S ~ Normal(μ={mu_S}, COV={cov_S})")

# 转换为对数正态参数
lambda_R = np.log(mu_R / np.sqrt(1 + cov_R**2))
zeta_R = np.sqrt(np.log(1 + cov_R**2))

# 蒙特卡洛估计失效概率
n_mc = 100000
np.random.seed(123)

R_samples = np.random.lognormal(lambda_R, zeta_R, n_mc)
S_samples = np.random.normal(mu_S, sigma_S, n_mc)

g_samples = R_samples - S_samples

# 失效概率估计
failure_indicator = (g_samples < 0).astype(int)
Pf_mc = np.mean(failure_indicator)
Pf_mc_std = np.std(failure_indicator) / np.sqrt(n_mc)

print(f"\n蒙特卡洛估计（N={n_mc}）：")
print(f"  失效概率 Pf = {Pf_mc:.6f}")
print(f"  95%置信区间 = [{Pf_mc - 1.96*Pf_mc_std:.6f}, {Pf_mc + 1.96*Pf_mc_std:.6f}]")

# 可靠指标
beta_mc = -stats.norm.ppf(Pf_mc)
print(f"  可靠指标 β = {beta_mc:.4f}")

# FORM方法（简化实现）
print("\nFORM分析：")

# 等效正态变换
# 对于对数正态变量R，在验算点处的等效正态参数
# 这里使用简化方法，直接在工作空间进行迭代

def limit_state(u):
    """
    标准正态空间中的极限状态函数
    u = [u_R, u_S]
    """
    # 逆变换到原始空间
    R = stats.lognorm.ppf(stats.norm.cdf(u[0]), zeta_R, scale=np.exp(lambda_R))
    S = stats.norm.ppf(stats.norm.cdf(u[1]), mu_S, sigma_S)
    return R - S

# 使用HL-RF算法寻找验算点
def hl_rf_algorithm(g_func, u0, max_iter=100, tol=1e-6):
    """HL-RF算法求解验算点"""
    u = u0.copy()
    
    for i in range(max_iter):
        # 计算函数值和梯度
        g_val = g_func(u)
        
        # 数值梯度
        eps = 1e-6
        grad = np.zeros_like(u)
        for j in range(len(u)):
            u_plus = u.copy()
            u_plus[j] += eps
            grad[j] = (g_func(u_plus) - g_val) / eps
        
        # 归一化梯度
        grad_norm = np.linalg.norm(grad)
        if grad_norm < 1e-10:
            break
        
        alpha = grad / grad_norm
        
        # 更新验算点
        u_new = (np.dot(grad, u) - g_val) / grad_norm * alpha
        
        # 检查收敛
        if np.linalg.norm(u_new - u) < tol:
            u = u_new
            break
        
        u = u_new
    
    return u, i+1

# 初始猜测
u0 = np.array([0.0, 0.0])
u_star, n_iter = hl_rf_algorithm(limit_state, u0)

beta_form = np.linalg.norm(u_star)
Pf_form = stats.norm.cdf(-beta_form)

print(f"  验算点 u* = [{u_star[0]:.4f}, {u_star[1]:.4f}]")
print(f"  可靠指标 β = {beta_form:.4f}")
print(f"  失效概率 Pf = {Pf_form:.6f}")
print(f"  迭代次数 = {n_iter}")

# 逆变换到原始空间
R_star = stats.lognorm.ppf(stats.norm.cdf(u_star[0]), zeta_R, scale=np.exp(lambda_R))
S_star = stats.norm.ppf(stats.norm.cdf(u_star[1]), mu_S, sigma_S)
print(f"  验算点（原始空间）：R* = {R_star:.2f}, S* = {S_star:.2f}")

# 绘制可靠性分析结果
fig6, axes = plt.subplots(2, 2, figsize=(14, 12))

# 散点图：R vs S
axes[0, 0].scatter(S_samples[::100], R_samples[::100], alpha=0.3, s=1, c='blue', label='Samples')
# 极限状态线
s_range = np.linspace(min(S_samples), max(S_samples), 100)
axes[0, 0].plot(s_range, s_range, 'r-', linewidth=2, label='g=0 (Limit State)')
axes[0, 0].scatter(S_star, R_star, c='red', s=100, marker='*', zorder=5, label='Design Point')
axes[0, 0].fill_between(s_range, s_range, max(R_samples), alpha=0.2, color='red', label='Failure Region')
axes[0, 0].set_xlabel('Load Effect S (kN)', fontsize=11)
axes[0, 0].set_ylabel('Resistance R (kN)', fontsize=11)
axes[0, 0].set_title('Reliability Analysis: R vs S', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 极限状态函数分布
axes[0, 1].hist(g_samples, bins=100, density=True, alpha=0.7, color='steelblue', edgecolor='black')
axes[0, 1].axvline(0, color='r', linestyle='--', linewidth=2, label='g=0')
axes[0, 1].axvline(np.mean(g_samples), color='g', linestyle='-', linewidth=2, 
                   label=f'Mean={np.mean(g_samples):.2f}')
axes[0, 1].set_xlabel('Limit State Function g(R,S)', fontsize=11)
axes[0, 1].set_ylabel('PDF', fontsize=11)
axes[0, 1].set_title('Distribution of g(R,S)', fontsize=12, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 标准正态空间
u1_range = np.linspace(-4, 4, 100)
u2_range = np.linspace(-4, 4, 100)
U1, U2 = np.meshgrid(u1_range, u2_range)

# 计算极限状态函数值
G = np.zeros_like(U1)
for i in range(len(u1_range)):
    for j in range(len(u2_range)):
        G[j, i] = limit_state([U1[j, i], U2[j, i]])

contour = axes[1, 0].contour(U1, U2, G, levels=[0], colors='red', linewidths=2)
axes[1, 0].clabel(contour, inline=True, fontsize=10)
axes[1, 0].scatter(u_star[0], u_star[1], c='red', s=100, marker='*', zorder=5, label='Design Point')
axes[1, 0].plot([0, u_star[0]], [0, u_star[1]], 'g--', linewidth=2, label=f'β={beta_form:.2f}')
# 绘制可靠度圆
theta = np.linspace(0, 2*np.pi, 100)
axes[1, 0].plot(beta_form*np.cos(theta), beta_form*np.sin(theta), 'k:', linewidth=1, alpha=0.5)
axes[1, 0].set_xlabel('u_R (Standard Normal)', fontsize=11)
axes[1, 0].set_ylabel('u_S (Standard Normal)', fontsize=11)
axes[1, 0].set_title('Standard Normal Space', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].axis('equal')

# 方法对比
methods = ['MCS', 'FORM']
betas = [beta_mc, beta_form]
Pfs = [Pf_mc, Pf_form]

x_pos = np.arange(len(methods))
width = 0.35

ax_twin = axes[1, 1].twinx()
bars1 = axes[1, 1].bar(x_pos - width/2, betas, width, label='Reliability Index β', 
                       color='steelblue', alpha=0.7)
bars2 = ax_twin.bar(x_pos + width/2, Pfs, width, label='Failure Probability Pf', 
                    color='coral', alpha=0.7)

axes[1, 1].set_xlabel('Method', fontsize=11)
axes[1, 1].set_ylabel('Reliability Index β', fontsize=11, color='steelblue')
ax_twin.set_ylabel('Failure Probability Pf', fontsize=11, color='coral')
axes[1, 1].set_title('Method Comparison', fontsize=12, fontweight='bold')
axes[1, 1].set_xticks(x_pos)
axes[1, 1].set_xticklabels(methods)
axes[1, 1].grid(True, alpha=0.3, axis='y')

# 添加数值标签
for bar, val in zip(bars1, betas):
    axes[1, 1].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.05,
                    f'{val:.3f}', ha='center', va='bottom', fontsize=9)
for bar, val in zip(bars2, Pfs):
    ax_twin.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.001,
                 f'{val:.4f}', ha='center', va='bottom', fontsize=9)

axes[1, 1].legend(loc='upper left')
ax_twin.legend(loc='upper right')

plt.tight_layout()
plt.savefig(f'{output_dir}/reliability_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  可靠性分析图已保存")

5.5 案例5：时变可靠性分析

本案例分析结构在随机时变荷载下的可靠性。

# ============================================================
# 案例5：时变可靠性分析
# ============================================================
print("\n案例5：时变可靠性分析")
print("-" * 70)

# 单自由度系统在随机地震作用下的可靠性分析
# 假设结构失效定义为位移超过阈值

# 系统参数（确定性）
m_tr = 1000.0  # kg
k_tr = 40000.0  # N/m
omega_n_tr = np.sqrt(k_tr / m_tr)
zeta_tr = 0.05
c_tr = 2 * zeta_tr * np.sqrt(k_tr * m_tr)

# 失效阈值
u_max = 0.15  # 最大允许位移 (m)

print("时变可靠性分析参数：")
print(f"  质量 m = {m_tr} kg")
print(f"  刚度 k = {k_tr} N/m")
print(f"  阻尼比 ζ = {zeta_tr}")
print(f"  失效阈值 u_max = {u_max} m")

# 生成随机地震动样本
def generate_random_earthquake(t, dt, seed=None):
    """生成随机地震动（白噪声通过滤波器）"""
    if seed is not None:
        np.random.seed(seed)
    
    n_points = len(t)
    
    # 生成白噪声
    white_noise = np.random.randn(n_points)
    
    # 简单滤波（模拟地震动频谱特性）
    from scipy import signal
    # 带通滤波 0.5-10 Hz
    sos = signal.butter(4, [0.5, 10], 'bandpass', fs=1/dt, output='sos')
    filtered = signal.sosfilt(sos, white_noise)
    
    # 包络函数
    envelope = np.ones(n_points)
    rise_idx = int(2.0 / dt)
    fall_start_idx = int(10.0 / dt)
    
    for i in range(n_points):
        if i < rise_idx:
            envelope[i] = (i / rise_idx) ** 2
        elif i > fall_start_idx:
            envelope[i] = np.exp(-0.5 * (i - fall_start_idx) * dt)
    
    earthquake = filtered * envelope
    
    # 归一化并缩放
    earthquake = earthquake / np.max(np.abs(earthquake)) * 3.0  # 峰值3 m/s²
    
    return earthquake

# 时间参数
t_tr = np.linspace(0, 20, 2000)
dt_tr = t_tr[1] - t_tr[0]

# 蒙特卡洛模拟时变可靠性
n_samples_tr = 1000
print(f"\n进行时变可靠性蒙特卡洛模拟（N={n_samples_tr}）...")

max_displacements = []

for i in range(n_samples_tr):
    # 生成随机地震动
    earthquake = generate_random_earthquake(t_tr, dt_tr, seed=i)
    
    # 等效荷载
    p_eq = -m_tr * earthquake
    
    # 求解系统响应（Newmark-beta法）
    def solve_sdof(m, c, k, p, dt):
        n = len(p)
        u = np.zeros(n)
        v = np.zeros(n)
        a = np.zeros(n)
        
        # 初始条件
        a[0] = p[0] / m
        
        # Newmark-beta参数（平均加速度法）
        beta = 0.25
        gamma = 0.5
        
        k_eff = k + gamma / (beta * dt) * c + 1 / (beta * dt**2) * m
        
        for j in range(n - 1):
            # 等效荷载
            p_eff = p[j+1] + m/(beta*dt**2)*u[j] + m/(beta*dt)*v[j] + (0.5/beta-1)*m*a[j] + \
                    gamma/(beta*dt)*c*u[j] + (gamma/beta-1)*c*v[j] + dt*(gamma/(2*beta)-1)*c*a[j]
            
            u[j+1] = p_eff / k_eff
            v[j+1] = gamma/(beta*dt)*(u[j+1]-u[j]) + (1-gamma/beta)*v[j] + dt*(1-0.5*gamma/beta)*a[j]
            a[j+1] = (p[j+1] - c*v[j+1] - k*u[j+1]) / m
        
        return u, v, a
    
    u_resp, _, _ = solve_sdof(m_tr, c_tr, k_tr, p_eq, dt_tr)
    max_displacements.append(np.max(np.abs(u_resp)))
    
    if (i + 1) % 200 == 0:
        print(f"  完成 {(i+1)/n_samples_tr*100:.0f}%")

max_displacements = np.array(max_displacements)

# 计算失效概率
failure_tr = max_displacements > u_max
Pf_tr = np.mean(failure_tr)
beta_tr = -stats.norm.ppf(Pf_tr)

print(f"\n时变可靠性分析结果：")
print(f"  最大位移均值 = {np.mean(max_displacements):.4f} m")
print(f"  最大位移标准差 = {np.std(max_displacements):.4f} m")
print(f"  失效概率 Pf = {Pf_tr:.6f}")
print(f"  可靠指标 β = {beta_tr:.4f}")

# 绘制时变可靠性分析结果
fig7, axes = plt.subplots(2, 2, figsize=(14, 12))

# 最大位移分布
axes[0, 0].hist(max_displacements, bins=50, density=True, alpha=0.7, color='steelblue', edgecolor='black')
axes[0, 0].axvline(u_max, color='r', linestyle='--', linewidth=2, label=f'Limit = {u_max} m')
axes[0, 0].axvline(np.mean(max_displacements), color='g', linestyle='-', linewidth=2, 
                   label=f'Mean = {np.mean(max_displacements):.4f} m')
axes[0, 0].set_xlabel('Max Displacement (m)', fontsize=11)
axes[0, 0].set_ylabel('PDF', fontsize=11)
axes[0, 0].set_title('Distribution of Max Displacement', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 累积分布函数
sorted_disp = np.sort(max_displacements)
cdf = np.arange(1, len(sorted_disp) + 1) / len(sorted_disp)
axes[0, 1].plot(sorted_disp, 1 - cdf, 'b-', linewidth=2)
axes[0, 1].axvline(u_max, color='r', linestyle='--', linewidth=2, label=f'Limit = {u_max} m')
axes[0, 1].axhline(Pf_tr, color='r', linestyle=':', linewidth=1, label=f'Pf = {Pf_tr:.4f}')
axes[0, 1].set_xlabel('Max Displacement (m)', fontsize=11)
axes[0, 1].set_ylabel('Exceedance Probability', fontsize=11)
axes[0, 1].set_title('Exceedance Probability Curve', fontsize=12, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 典型响应时程（选取几个样本）
for idx in range(5):
    earthquake = generate_random_earthquake(t_tr, dt_tr, seed=idx)
    p_eq = -m_tr * earthquake
    u_resp, _, _ = solve_sdof(m_tr, c_tr, k_tr, p_eq, dt_tr)
    axes[1, 0].plot(t_tr, u_resp, alpha=0.7, linewidth=0.8, label=f'Sample {idx+1}')
axes[1, 0].axhline(u_max, color='r', linestyle='--', linewidth=2, label=f'Limit = {u_max} m')
axes[1, 0].axhline(-u_max, color='r', linestyle='--', linewidth=2)
axes[1, 0].set_xlabel('Time (s)', fontsize=11)
axes[1, 0].set_ylabel('Displacement (m)', fontsize=11)
axes[1, 0].set_title('Typical Response Time Histories', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=8)
axes[1, 0].grid(True, alpha=0.3)

# 失效样本的响应
failure_indices = np.where(failure_tr)[0]
if len(failure_indices) > 0:
    for idx in failure_indices[:3]:
        earthquake = generate_random_earthquake(t_tr, dt_tr, seed=idx)
        p_eq = -m_tr * earthquake
        u_resp, _, _ = solve_sdof(m_tr, c_tr, k_tr, p_eq, dt_tr)
        axes[1, 1].plot(t_tr, u_resp, alpha=0.7, linewidth=0.8, label=f'Failure {idx+1}')
axes[1, 1].axhline(u_max, color='r', linestyle='--', linewidth=2, label=f'Limit = {u_max} m')
axes[1, 1].axhline(-u_max, color='r', linestyle='--', linewidth=2)
axes[1, 1].set_xlabel('Time (s)', fontsize=11)
axes[1, 1].set_ylabel('Displacement (m)', fontsize=11)
axes[1, 1].set_title('Failure Cases', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=8)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/time_variant_reliability.png', dpi=150, bbox_inches='tight')
plt.close()
print("  时变可靠性分析图已保存")

print("\n" + "=" * 70)
print("不确定性分析完成！")
print("=" * 70)

print("\n生成的文件：")
print("  1. mc_parameter_distributions.png - 参数分布图")
print("  2. mc_scatter_matrix.png - 散点图矩阵")
print("  3. mc_convergence.png - 收敛性分析")
print("  4. lhs_comparison.png - LHS对比分析")
print("  5. sensitivity_analysis.png - 敏感性分析")
print("  6. reliability_analysis.png - 可靠性分析")
print("  7. time_variant_reliability.png - 时变可靠性分析")

---

6. 结果分析与讨论

6.1 蒙特卡洛模拟结果分析

从蒙特卡洛模拟结果可以看出，输入参数的随机性通过系统传递后，响应的统计特性发生了显著变化。质量服从正态分布，刚度服从对数正态分布，而响应振幅的分布则呈现出复杂的非对称形态。这种不确定性传播的非线性特征在实际工程中具有重要意义。

收敛性分析表明，当样本量达到3000以上时，响应均值和标准差的估计基本稳定。这说明对于本案例，5000个样本足以获得可靠的统计估计。在实际应用中，应根据精度要求和计算资源合理选择样本量。

6.2 抽样方法对比

拉丁超立方抽样与简单随机抽样的对比显示，LHS在相同样本量下能够提供更稳定的估计，特别是对小样本情况。LHS的空间填充特性使其在参数敏感性分析和可靠性评估中具有优势。然而，LHS的优势随样本量增加而减弱，对于大样本蒙特卡洛模拟，两种方法的差异变得不明显。

6.3 敏感性分析结果

敏感性分析结果表明，在刚度变化范围内，刚度对响应的影响最为显著，其次是阻尼比，质量的影响相对较小。这一发现对工程设计具有指导意义：在需要控制结构响应时，应优先考虑调整刚度参数。

标准化回归系数（SRC）与相关系数的结果一致，验证了敏感性分析的可靠性。局部敏感性分析与全局敏感性分析的结果趋势一致，说明在本案例的参数范围内，系统响应近似线性依赖于输入参数。

6.4 可靠性分析结果

FORM方法成功找到了验算点，计算得到的可靠指标与蒙特卡洛结果吻合良好。验算点位于标准正态空间中距离原点最近的极限状态曲面上，这一几何解释直观地展示了可靠性的物理意义。

时变可靠性分析表明，在随机地震作用下，结构的失效概率可以通过蒙特卡洛模拟有效估计。失效样本的响应时程显示，失效通常由地震动的强震段引起，这与工程经验一致。

---

7. 工程应用与扩展

7.1 基于可靠性的设计优化

不确定性分析为基于可靠性的设计优化（RBDO）提供了理论基础。RBDO的目标是在满足可靠性约束的前提下，最小化结构成本或重量。数学表述为：

$$\underset{d}{\min }f(d)$$

$$\text{s.t.}P(g_i(X,d)\le 0)\le P_{f,i}^{target},i=1,2,...,m$$

其中，$d$ 为设计变量（如截面尺寸），$X$ 为随机变量，$P_{f,i}^{target}$ 为目标失效概率。

7.2 模型不确定性量化

除了参数不确定性，模型形式的不确定性也是工程实践中需要考虑的重要因素。贝叶斯模型平均（Bayesian Model Averaging）提供了一种系统的方法来量化和传播模型不确定性：

$$p(Y|D)=\sum _{k=1}^{K}p(Y|M_k,D)p(M_k|D)$$

其中，$M_k$ 为第 $k$ 个候选模型，$p(M_k|D)$ 为模型后验概率。

7.3 前沿研究方向

1. 深度学习在不确定性量化中的应用：利用神经网络代理模型加速蒙特卡洛模拟
2. 高维不确定性量化：处理数百甚至数千个随机变量的高效算法
3. 非概率不确定性分析：区间分析、模糊理论、证据理论等方法的发展
4. 数字孪生与实时不确定性更新：结合监测数据动态更新不确定性模型

---

8. 总结

本教程系统地介绍了结构动力学中的不确定性分析理论与方法。主要内容包括：

1. 不确定性理论基础：阐述了偶然不确定性和认知不确定性的概念，介绍了概率论基础和结构可靠性理论。

2. 随机抽样方法：详细讲解了蒙特卡洛模拟、拉丁超立方抽样等方法的原理和实现。

3. 敏感性分析方法：介绍了局部敏感性分析和全局敏感性分析的各种技术。

4. 结构可靠性分析：讲解了FORM、SORM、响应面法等可靠性评估方法。

5. Python实现案例：通过五个典型案例，展示了不确定性分析在实际问题中的应用。

通过本教程的学习，读者应能够：

• 识别和分类工程问题中的不确定性来源
• 选择合适的概率分布描述随机变量
• 运用蒙特卡洛模拟和LHS进行不确定性传播分析
• 执行敏感性分析识别关键参数
• 评估结构的可靠性和失效概率
• 使用Python实现不确定性分析算法

不确定性分析是现代结构动力学不可或缺的一部分，它为工程设计提供了更科学、更全面的决策依据。随着计算能力的提升和算法的发展，不确定性分析将在结构工程中发挥越来越重要的作用。