第六十四篇：结构动力学不确定性量化

摘要

不确定性量化（Uncertainty Quantification, UQ）是现代结构动力学分析的重要组成部分。本主题系统介绍结构动力学中不确定性的来源、分类和量化方法。重点讲解蒙特卡洛模拟、多项式混沌展开、随机有限元方法等核心算法，以及结构可靠性分析和敏感性分析方法。通过Python实现各种不确定性量化技术，帮助读者掌握处理工程实际中参数不确定性、模型不确定性和载荷不确定性的能力，为结构安全评估和优化设计提供理论基础。

关键词

不确定性量化；蒙特卡洛模拟；多项式混沌；随机有限元；可靠性分析；敏感性分析；结构安全

---

1. 不确定性量化基础理论

1.1 不确定性的来源与分类

在结构动力学分析中，不确定性主要来源于以下几个方面：

参数不确定性：材料属性（弹性模量、密度、阻尼系数）、几何尺寸、边界条件等物理参数的固有变异性。

模型不确定性：数学模型对真实物理过程的简化假设、模型形式的选择、以及模型参数的识别误差。

载荷不确定性：外部激励（地震、风载、波浪等）的随机性和不可预测性。

测量不确定性：实验数据采集过程中的噪声和系统误差。

1.2 不确定性表征方法

概率方法：将不确定性参数视为随机变量，用概率密度函数描述其统计特性。

$$X\sim f_X(x),x\in \mathbb{R}$$

区间方法：当概率信息不足时，用区间范围描述不确定性。

$$X\in [\underline{x},\overline{x}]$$

模糊方法：使用模糊集合理论处理认知不确定性。

1.3 不确定性传播

不确定性传播研究输入不确定性如何通过系统模型传递到输出响应。设系统模型为：

$$Y=g(X)$$

其中 $X$ 为输入随机变量，$Y$ 为输出响应。不确定性传播的目标是计算输出 $Y$ 的统计特性（均值、方差、概率分布等）。

---

2. 蒙特卡洛模拟方法

2.1 基本理论

蒙特卡洛模拟是最直观的不确定性量化方法，通过大量随机采样来估计输出统计特性。

算法步骤：

1. 从输入随机变量的概率分布中生成 $N$ 个样本
2. 对每个样本计算系统响应
3. 统计分析输出结果

统计估计：

均值估计：$${\hat{\mu }}_Y=\frac{1}{N}\sum _{i=1}^N{Y}_i$$

方差估计：$${\hat{\sigma }}_Y^2=\frac{1}{N-1}\sum _{i=1}^N(Y_i-{\hat{\mu }}_Y{)}^2$$

2.2 方差缩减技术

重要性采样：通过改变采样分布来减少方差。

分层采样：将样本空间划分为若干层，在每层内均匀采样。

拉丁超立方采样：确保样本在每一维度上均匀分布。

---

3. 多项式混沌展开

3.1 基本理论

多项式混沌展开（Polynomial Chaos Expansion, PCE）将随机输出表示为随机输入的正交多项式展开：

$$Y(\xi )=\sum _{i=0}^P{y}_i{\Psi }_i(\xi )$$

其中 ${\Psi }_i(\xi )$ 是关于随机变量 $\xi$ 的正交多项式，$y_i$ 是展开系数。

正交多项式选择：

• 高斯随机变量 $\to$ Hermite多项式
• 均匀随机变量 $\to$ Legendre多项式
• 伽马随机变量 $\to$ Laguerre多项式

3.2 系数计算方法

投影法：利用正交性计算系数

$$y_i=\frac{⟨Y,{\Psi }_i⟩}{⟨{\Psi }_i,{\Psi }_i⟩}$$

回归法：通过最小二乘拟合确定系数

---

4. 随机有限元方法

4.1 摄动法

摄动法将随机变量展开为确定性部分和小扰动部分：

$$X=X_0+ϵX_1$$

其中 $ϵ$ 是小参数。将展开式代入控制方程，按 $ϵ$ 的幂次整理，得到各阶摄动方程。

4.2 Neumann展开法

对于线性系统 $KU=F$，当刚度矩阵 $K$ 具有随机性时，使用Neumann级数展开：

$$K^{-1}=K_0^{-1}\sum _{i=0}^{\infty }(-K_0^{-1}\Delta K{)}^i$$

---

5. 可靠性分析

5.1 失效概率计算

结构失效概率定义为：

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

其中 $g(X)$ 是极限状态函数，$g(X)>0$ 表示安全，$g(X)\le 0$ 表示失效。

5.2 一阶可靠性方法（FORM）

FORM通过寻找设计点（最可能失效点）来近似计算失效概率。

可靠性指标：

$$\beta =\underset{g(u)=0}{\min }\sqrt{u^{T}u}$$

失效概率近似：

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

5.3 重要性测度

Sobol指数：衡量输入变量对输出方差的贡献程度。

$$S_i=\frac{V_{X_i}(E_{X_{\sim i}}[Y|X_i])}{V(Y)}$$

---

6. Python实现案例

6.1 案例1：蒙特卡洛模拟方法

问题描述：
分析简支梁在随机载荷作用下的响应不确定性。假设弹性模量和载荷幅值为随机变量，使用蒙特卡洛模拟计算梁中点挠度的统计特性。

Python代码实现：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
from scipy import stats
from scipy.stats import norm, uniform

class MonteCarloSimulation:
    """蒙特卡洛模拟类"""
    
    def __init__(self, n_samples=10000):
        """
        初始化
        
        Parameters:
        -----------
        n_samples : int
            蒙特卡洛样本数
        """
        self.n_samples = n_samples
    
    def latin_hypercube_sampling(self, distributions, n_samples=None):
        """
        拉丁超立方采样
        
        Parameters:
        -----------
        distributions : list
            随机变量分布列表
        n_samples : int
            样本数
        
        Returns:
        --------
        samples : array
            采样矩阵 (n_samples x n_variables)
        """
        if n_samples is None:
            n_samples = self.n_samples
        
        n_vars = len(distributions)
        samples = np.zeros((n_samples, n_vars))
        
        for i, dist in enumerate(distributions):
            # 分层
            bins = np.linspace(0, 1, n_samples + 1)
            # 在每个层内随机采样
            u = np.random.rand(n_samples)
            u = (np.arange(n_samples) + u) / n_samples
            # 随机打乱
            np.random.shuffle(u)
            # 逆变换采样
            samples[:, i] = dist.ppf(u)
        
        return samples
    
    def simple_random_sampling(self, distributions, n_samples=None):
        """
        简单随机采样
        
        Parameters:
        -----------
        distributions : list
            随机变量分布列表
        n_samples : int
            样本数
        
        Returns:
        --------
        samples : array
            采样矩阵
        """
        if n_samples is None:
            n_samples = self.n_samples
        
        n_vars = len(distributions)
        samples = np.zeros((n_samples, n_vars))
        
        for i, dist in enumerate(distributions):
            samples[:, i] = dist.rvs(n_samples)
        
        return samples

def simply_supported_beam_deflection(E, P, L=2.0, b=0.05, h=0.1):
    """
    计算简支梁中点挠度
    
    Parameters:
    -----------
    E : float or array
        弹性模量 (Pa)
    P : float or array
        集中载荷 (N)
    L : float
        梁长度 (m)
    b : float
        截面宽度 (m)
    h : float
        截面高度 (m)
    
    Returns:
    --------
    delta : float or array
        中点挠度 (m)
    """
    I = b * h**3 / 12  # 截面惯性矩
    delta = P * L**3 / (48 * E * I)
    return delta

def case1_monte_carlo():
    """案例1：蒙特卡洛模拟方法"""
    print("="*60)
    print("主题064 - 案例1: 蒙特卡洛模拟方法")
    print("="*60)
    
    # 定义随机变量
    # 弹性模量：正态分布，均值210GPa，变异系数5%
    E_mean = 210e9
    E_std = 0.05 * E_mean
    E_dist = norm(E_mean, E_std)
    
    # 载荷：正态分布，均值1000N，变异系数10%
    P_mean = 1000
    P_std = 0.10 * P_mean
    P_dist = norm(P_mean, P_std)
    
    distributions = [E_dist, P_dist]
    
    # 确定性分析（参考值）
    delta_det = simply_supported_beam_deflection(E_mean, P_mean)
    print(f"\n确定性分析结果:")
    print(f"  中点挠度 = {delta_det*1000:.4f} mm")
    
    # 蒙特卡洛模拟
    print("\n蒙特卡洛模拟...")
    
    # 方法1：简单随机采样
    mc = MonteCarloSimulation(n_samples=10000)
    samples_srs = mc.simple_random_sampling(distributions)
    
    E_samples = samples_srs[:, 0]
    P_samples = samples_srs[:, 1]
    
    delta_samples_srs = simply_supported_beam_deflection(E_samples, P_samples)
    
    # 统计结果
    mean_srs = np.mean(delta_samples_srs)
    std_srs = np.std(delta_samples_srs)
    cov_srs = std_srs / mean_srs
    
    print(f"\n简单随机采样结果 (N=10000):")
    print(f"  均值 = {mean_srs*1000:.4f} mm")
    print(f"  标准差 = {std_srs*1000:.4f} mm")
    print(f"  变异系数 = {cov_srs:.4f}")
    
    # 方法2：拉丁超立方采样
    samples_lhs = mc.latin_hypercube_sampling(distributions)
    
    E_samples_lhs = samples_lhs[:, 0]
    P_samples_lhs = samples_lhs[:, 1]
    
    delta_samples_lhs = simply_supported_beam_deflection(E_samples_lhs, P_samples_lhs)
    
    mean_lhs = np.mean(delta_samples_lhs)
    std_lhs = np.std(delta_samples_lhs)
    cov_lhs = std_lhs / mean_lhs
    
    print(f"\n拉丁超立方采样结果 (N=10000):")
    print(f"  均值 = {mean_lhs*1000:.4f} mm")
    print(f"  标准差 = {std_lhs*1000:.4f} mm")
    print(f"  变异系数 = {cov_lhs:.4f}")
    
    # 收敛性分析
    sample_sizes = [100, 500, 1000, 2000, 5000, 10000]
    means_convergence = []
    stds_convergence = []
    
    for n in sample_sizes:
        samples = mc.simple_random_sampling(distributions, n)
        deltas = simply_supported_beam_deflection(samples[:, 0], samples[:, 1])
        means_convergence.append(np.mean(deltas))
        stds_convergence.append(np.std(deltas))
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 1. 样本散点图
    ax = axes[0, 0]
    ax.scatter(E_samples[:1000]/1e9, P_samples[:1000], alpha=0.5, s=10)
    ax.set_xlabel('弹性模量 E (GPa)')
    ax.set_ylabel('载荷 P (N)')
    ax.set_title('输入随机变量样本')
    ax.grid(True, alpha=0.3)
    
    # 2. 输出响应直方图
    ax = axes[0, 1]
    ax.hist(delta_samples_srs*1000, bins=50, density=True, alpha=0.7, 
            color='blue', edgecolor='black', label='SRS')
    ax.hist(delta_samples_lhs*1000, bins=50, density=True, alpha=0.5, 
            color='red', edgecolor='black', label='LHS')
    ax.axvline(delta_det*1000, color='green', linestyle='--', 
               linewidth=2, label='确定性值')
    ax.set_xlabel('中点挠度 (mm)')
    ax.set_ylabel('概率密度')
    ax.set_title('输出响应分布')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 3. 收敛性分析 - 均值
    ax = axes[0, 2]
    ax.plot(sample_sizes, np.array(means_convergence)*1000, 'bo-', linewidth=2)
    ax.axhline(mean_srs*1000, color='r', linestyle='--', alpha=0.7)
    ax.set_xlabel('样本数')
    ax.set_ylabel('均值 (mm)')
    ax.set_title('均值收敛性')
    ax.set_xscale('log')
    ax.grid(True, alpha=0.3)
    
    # 4. 收敛性分析 - 标准差
    ax = axes[1, 0]
    ax.plot(sample_sizes, np.array(stds_convergence)*1000, 'ro-', linewidth=2)
    ax.axhline(std_srs*1000, color='b', linestyle='--', alpha=0.7)
    ax.set_xlabel('样本数')
    ax.set_ylabel('标准差 (mm)')
    ax.set_title('标准差收敛性')
    ax.set_xscale('log')
    ax.grid(True, alpha=0.3)
    
    # 5. Q-Q图
    ax = axes[1, 1]
    stats.probplot(delta_samples_srs*1000, dist="norm", plot=ax)
    ax.set_title('正态Q-Q图')
    ax.grid(True, alpha=0.3)
    
    # 6. 统计结果对比
    ax = axes[1, 2]
    methods = ['SRS', 'LHS', '解析']
    means = [mean_srs*1000, mean_lhs*1000, delta_det*1000]
    stds = [std_srs*1000, std_lhs*1000, 0]
    
    x = np.arange(len(methods))
    width = 0.35
    
    ax.bar(x - width/2, means, width, label='均值', color='blue', alpha=0.8)
    ax.bar(x + width/2, stds, width, label='标准差', color='red', alpha=0.8)
    
    ax.set_ylabel('挠度 (mm)')
    ax.set_title('统计结果对比')
    ax.set_xticks(x)
    ax.set_xticklabels(methods)
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('monte_carlo_simulation.png', dpi=150, bbox_inches='tight')
    print("\n蒙特卡洛模拟结果图已保存: monte_carlo_simulation.png")
    plt.close()
    
    # 创建动画
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    ax1 = axes[0]
    scatter = ax1.scatter([], [], alpha=0.5, s=20)
    ax1.set_xlim(180, 240)
    ax1.set_ylim(700, 1300)
    ax1.set_xlabel('弹性模量 E (GPa)')
    ax1.set_ylabel('载荷 P (N)')
    ax1.set_title('样本累积过程')
    ax1.grid(True, alpha=0.3)
    
    ax2 = axes[1]
    line_mean, = ax2.plot([], [], 'b-', linewidth=2, label='均值')
    line_std, = ax2.plot([], [], 'r-', linewidth=2, label='标准差')
    ax2.set_xlim(0, 10000)
    ax2.set_ylim(0.8, 1.5)
    ax2.set_xlabel('样本数')
    ax2.set_ylabel('挠度 (mm)')
    ax2.set_title('统计量收敛过程')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    def animate(frame):
        n_show = min((frame + 1) * 200, 10000)
        
        # 更新散点
        scatter.set_offsets(np.column_stack([E_samples[:n_show]/1e9, P_samples[:n_show]]))
        
        # 更新收敛曲线
        if n_show >= 100:
            sample_sizes_anim = np.arange(100, n_show + 1, 100)
            means_anim = []
            stds_anim = []
            for n in sample_sizes_anim:
                means_anim.append(np.mean(delta_samples_srs[:n])*1000)
                stds_anim.append(np.std(delta_samples_srs[:n])*1000)
            
            line_mean.set_data(sample_sizes_anim, means_anim)
            line_std.set_data(sample_sizes_anim, stds_anim)
        
        return scatter, line_mean, line_std
    
    anim = FuncAnimation(fig, animate, frames=50, interval=100, blit=False)
    anim.save('monte_carlo_animation.gif', writer='pillow', fps=10, dpi=100)
    print("蒙特卡洛模拟动画已保存: monte_carlo_animation.gif")
    plt.close()
    
    return mc, delta_samples_srs

if __name__ == "__main__":
    mc, delta_samples = case1_monte_carlo()

6.2 案例2：多项式混沌展开方法

问题描述：
使用多项式混沌展开（PCE）方法分析单自由度振动系统的频率不确定性。系统刚度服从对数正态分布，通过PCE高效计算固有频率的统计矩。

理论基础：

多项式混沌展开将随机响应表示为随机输入的正交多项式级数：

$$Y(\xi )=\sum _{i=0}^P{y}_i{\Psi }_i(\xi )$$

其中 ${\Psi }_i(\xi )$ 是正交多项式基函数，$y_i$ 是展开系数。

对于高斯随机变量，使用Hermite多项式：

$$H_n(\xi )=(-1{)}^n{e}^{{\xi }^2/2}\frac{d^n}{d{\xi }^n}{e}^{-{\xi }^2/2}$$

统计矩计算：

均值：$${\mu }_Y=y_0$$

方差：$${\sigma }_Y^2=\sum _{i=1}^P{y}_i^2⟨{\Psi }_i^2⟩$$

Python代码实现：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import eval_hermitenorm
from scipy import stats

class PolynomialChaosExpansion:
    """多项式混沌展开类"""
    
    def __init__(self, order=3):
        """
        初始化
        
        Parameters:
        -----------
        order : int
            多项式阶数
        """
        self.order = order
        self.coefficients = None
        self.poly_type = None
    
    def hermite_polynomial(self, n, x):
        """计算归一化Hermite多项式"""
        return eval_hermitenorm(n, x)
    
    def generate_polynomial_basis(self, xi, dist_type='normal'):
        """
        生成多项式基函数矩阵
        
        Parameters:
        -----------
        xi : array
            标准随机变量样本
        dist_type : str
            分布类型
            
        Returns:
        --------
        basis : array
            多项式基函数矩阵
        """
        n_samples = len(xi)
        n_terms = self.order + 1
        basis = np.zeros((n_samples, n_terms))
        
        for i in range(n_terms):
            basis[:, i] = self.hermite_polynomial(i, xi)
        
        return basis
    
    def fit_regression(self, xi, y, dist_type='normal'):
        """
        使用回归法计算PCE系数
        
        Parameters:
        -----------
        xi : array
            标准随机变量
        y : array
            输出响应
        dist_type : str
            分布类型
        """
        basis = self.generate_polynomial_basis(xi, dist_type)
        coefficients = np.linalg.lstsq(basis, y, rcond=None)[0]
        self.coefficients = coefficients
        self.poly_type = dist_type
        return coefficients
    
    def predict(self, xi):
        """使用PCE模型预测"""
        if self.coefficients is None:
            raise ValueError("请先调用fit_regression方法拟合模型")
        basis = self.generate_polynomial_basis(xi, self.poly_type)
        return basis @ self.coefficients
    
    def calculate_statistics(self):
        """从PCE系数计算统计量"""
        if self.coefficients is None:
            raise ValueError("请先调用fit_regression方法拟合模型")
        
        # 均值（0阶系数）
        mean = self.coefficients[0]
        
        # 方差（高阶系数平方和）
        variance = np.sum(self.coefficients[1:]**2)
        std = np.sqrt(variance)
        cov = std / mean if mean != 0 else 0
        
        return {
            'mean': mean,
            'variance': variance,
            'std': std,
            'cov': cov
        }

def sdof_natural_frequency(k, m):
    """计算单自由度系统固有频率"""
    return np.sqrt(k / m) / (2 * np.pi)

def case2_polynomial_chaos():
    """案例2：多项式混沌展开"""
    print("="*60)
    print("主题064 - 案例2: 多项式混沌展开方法")
    print("="*60)
    
    # 系统参数
    m = 100.0  # 质量 (kg)
    
    # 刚度：对数正态分布
    k_mean = 10000.0  # 均值 (N/m)
    k_std = 2000.0    # 标准差
    
    # 对数正态分布参数
    sigma_ln = np.sqrt(np.log(1 + (k_std/k_mean)**2))
    mu_ln = np.log(k_mean) - 0.5 * sigma_ln**2
    
    # 确定性分析
    f_det = sdof_natural_frequency(k_mean, m)
    print(f"\n确定性分析结果:")
    print(f"  固有频率 = {f_det:.4f} Hz")
    
    # 解析解（对数正态分布的频率统计特性）
    f_mean_analytical = np.exp(mu_ln/2) * np.sqrt(1/m) / (2 * np.pi)
    f_var_analytical = (np.exp(sigma_ln**2) - 1) * f_mean_analytical**2
    
    print(f"\n解析解:")
    print(f"  均值 = {f_mean_analytical:.4f} Hz")
    print(f"  标准差 = {np.sqrt(f_var_analytical):.4f} Hz")
    print(f"  变异系数 = {np.sqrt(f_var_analytical)/f_mean_analytical:.4f}")
    
    # PCE方法
    print("\n多项式混沌展开分析...")
    
    # 生成训练样本
    n_train = 100
    xi_train = np.random.randn(n_train)
    k_train = np.exp(mu_ln + sigma_ln * xi_train)
    f_train = sdof_natural_frequency(k_train, m)
    
    # 拟合PCE模型（不同阶数）
    orders = [1, 2, 3, 4, 5]
    pce_results = []
    
    for order in orders:
        pce = PolynomialChaosExpansion(order=order)
        pce.fit_regression(xi_train, f_train, dist_type='normal')
        stats_dict = pce.calculate_statistics()
        pce_results.append({
            'order': order,
            'pce': pce,
            'stats': stats_dict
        })
        
        print(f"\nPCE阶数 = {order}:")
        print(f"  均值 = {stats_dict['mean']:.4f} Hz")
        print(f"  标准差 = {stats_dict['std']:.4f} Hz")
        print(f"  变异系数 = {stats_dict['cov']:.4f}")
    
    # 蒙特卡洛验证
    n_mc = 100000
    xi_mc = np.random.randn(n_mc)
    k_mc = np.exp(mu_ln + sigma_ln * xi_mc)
    f_mc = sdof_natural_frequency(k_mc, m)
    
    print(f"\n蒙特卡洛验证 (N={n_mc}):")
    print(f"  均值 = {np.mean(f_mc):.4f} Hz")
    print(f"  标准差 = {np.std(f_mc):.4f} Hz")
    print(f"  变异系数 = {np.std(f_mc)/np.mean(f_mc):.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 子图1：PCE模型预测 vs 真实值
    ax = axes[0, 0]
    xi_test = np.linspace(-3, 3, 100)
    k_test = np.exp(mu_ln + sigma_ln * xi_test)
    f_true = sdof_natural_frequency(k_test, m)
    
    for result in pce_results[::2]:  # 每隔一个显示
        order = result['order']
        pce = result['pce']
        f_pred = pce.predict(xi_test)
        ax.plot(xi_test, f_pred, '--', linewidth=2, label=f'PCE Order {order}')
    
    ax.plot(xi_test, f_true, 'k-', linewidth=2, label='True')
    ax.scatter(xi_train[:20], f_train[:20], c='red', s=50, alpha=0.5, label='Training Data')
    ax.set_xlabel('标准随机变量 ξ')
    ax.set_ylabel('固有频率 (Hz)')
    ax.set_title('PCE模型拟合')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 子图2：统计量收敛
    ax = axes[0, 1]
    means_pce = [r['stats']['mean'] for r in pce_results]
    stds_pce = [r['stats']['std'] for r in pce_results]
    
    ax.plot(orders, means_pce, 'bo-', linewidth=2, markersize=8, label='PCE Mean')
    ax.axhline(np.mean(f_mc), color='r', linestyle='--', linewidth=2, label='MC Mean')
    ax.set_xlabel('PCE阶数')
    ax.set_ylabel('均值 (Hz)')
    ax.set_title('均值收敛')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 子图3：标准差收敛
    ax = axes[1, 0]
    ax.plot(orders, stds_pce, 'ro-', linewidth=2, markersize=8, label='PCE Std')
    ax.axhline(np.std(f_mc), color='b', linestyle='--', linewidth=2, label='MC Std')
    ax.set_xlabel('PCE阶数')
    ax.set_ylabel('标准差 (Hz)')
    ax.set_title('标准差收敛')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 子图4：频率分布直方图
    ax = axes[1, 1]
    ax.hist(f_mc, bins=50, density=True, alpha=0.5, color='blue', 
            edgecolor='black', label='Monte Carlo')
    
    # PCE预测分布
    xi_pce = np.random.randn(10000)
    f_pce_best = pce_results[-1]['pce'].predict(xi_pce)
    ax.hist(f_pce_best, bins=50, density=True, alpha=0.5, color='red', 
            edgecolor='black', label='PCE Prediction')
    
    ax.axvline(f_det, color='green', linestyle='--', linewidth=2, label='Deterministic')
    ax.set_xlabel('固有频率 (Hz)')
    ax.set_ylabel('概率密度')
    ax.set_title('频率分布对比')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('pce_results.png', dpi=150, bbox_inches='tight')
    print("\n多项式混沌展开结果图已保存: pce_results.png")
    plt.close()

if __name__ == "__main__":
    case2_polynomial_chaos()

案例2结果分析：

多项式混沌展开方法展示了以下特点：

1. 计算效率：PCE方法仅需少量样本（如100个）即可建立代理模型，而蒙特卡洛方法需要数万个样本才能达到相同精度
2. 收敛特性：随着多项式阶数增加，PCE估计的统计量快速收敛到参考值
3. 分布拟合：高阶PCE能够准确捕捉输出响应的概率分布
4. 物理意义：本案例中频率与刚度的平方根成正比，PCE能够有效建模这种非线性关系

6.3 案例3：随机有限元方法

问题描述：
使用随机有限元方法（SFEM）分析悬臂梁的模态特性不确定性。考虑弹性模量和密度为空间随机场，通过Karhunen-Loève展开离散随机场，使用蒙特卡洛方法求解随机特征值问题。

理论基础：

Karhunen-Loève展开：

随机场 $H(x,\theta )$ 可以表示为：

$$H(x,\theta )=\bar{H}(x)+\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}{\xi }_i(\theta ){\varphi }_i(x)$$

其中 ${\lambda }_i$ 和 ${\varphi }_i(x)$ 是相关函数的特征值和特征函数，${\xi }_i$ 是标准正态随机变量。

相关函数（平方指数型）：

$$C(x_1,x_2)=\exp \left(-\frac{(x_1-x_2{)}^2}{L_c^2}\right)$$

其中 $L_c$ 是相关长度。

Python代码实现：

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg

class RandomField:
    """随机场类 - 使用Karhunen-Loève展开"""
    
    def __init__(self, mean, std, correlation_length, domain_length):
        self.mean = mean
        self.std = std
        self.correlation_length = correlation_length
        self.domain_length = domain_length
    
    def squared_exponential_correlation(self, x1, x2):
        """平方指数型（高斯）相关函数"""
        return np.exp(-((x1 - x2) / self.correlation_length)**2)
    
    def compute_kle_modes(self, n_modes, n_points=100):
        """计算Karhunen-Loève展开的特征值和特征函数"""
        x = np.linspace(0, self.domain_length, n_points)
        
        # 构建相关矩阵
        C = np.zeros((n_points, n_points))
        for i in range(n_points):
            for j in range(n_points):
                C[i, j] = self.squared_exponential_correlation(x[i], x[j])
        
        # 特征值分解
        eigenvalues, eigenvectors = np.linalg.eigh(C)
        
        # 按特征值降序排列
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx][:n_modes]
        eigenvectors = eigenvectors[:, idx][:, :n_modes]
        
        # 归一化特征向量
        for i in range(n_modes):
            eigenvectors[:, i] /= np.sqrt(np.trapezoid(eigenvectors[:, i]**2, x))
        
        return eigenvalues, eigenvectors, x
    
    def generate_realization(self, xi, eigenvalues, eigenvectors, x):
        """生成随机场的一个实现"""
        n_modes = len(xi)
        field = np.zeros_like(x)
        
        for i in range(n_modes):
            field += np.sqrt(eigenvalues[i]) * xi[i] * eigenvectors[:, i]
        
        # 转换为对数正态分布
        field = self.mean * np.exp(self.std * field - 0.5 * self.std**2)
        
        return field

class StochasticFEM:
    """随机有限元分析类"""
    
    def __init__(self, n_elements, length, E_mean, rho_mean, A=1.0, I=1.0):
        self.n_elements = n_elements
        self.length = length
        self.E_mean = E_mean
        self.rho_mean = rho_mean
        self.A = A
        self.I = I
        self.n_nodes = n_elements + 1
        self.h = length / n_elements
        self.x_nodes = np.linspace(0, length, self.n_nodes)
    
    def beam_element_matrices(self, E, rho, h):
        """计算欧拉-伯努利梁单元的刚度矩阵和质量矩阵"""
        # 刚度矩阵
        k = E * self.I / h**3 * np.array([
            [12, 6*h, -12, 6*h],
            [6*h, 4*h**2, -6*h, 2*h**2],
            [-12, -6*h, 12, -6*h],
            [6*h, 2*h**2, -6*h, 4*h**2]
        ])
        
        # 一致质量矩阵
        m = rho * self.A * h / 420 * np.array([
            [156, 22*h, 54, -13*h],
            [22*h, 4*h**2, 13*h, -3*h**2],
            [54, 13*h, 156, -22*h],
            [-13*h, -3*h**2, -22*h, 4*h**2]
        ])
        
        return k, m
    
    def assemble_system_matrices(self, E_field, rho_field):
        """组装全局刚度矩阵和质量矩阵"""
        ndof = 2 * self.n_nodes
        K = np.zeros((ndof, ndof))
        M = np.zeros((ndof, ndof))
        
        for e in range(self.n_elements):
            E_e = (E_field[e] + E_field[e+1]) / 2
            rho_e = (rho_field[e] + rho_field[e+1]) / 2
            
            k_e, m_e = self.beam_element_matrices(E_e, rho_e, self.h)
            dof = [2*e, 2*e+1, 2*(e+1), 2*(e+1)+1]
            
            for i in range(4):
                for j in range(4):
                    K[dof[i], dof[j]] += k_e[i, j]
                    M[dof[i], dof[j]] += m_e[i, j]
        
        return K, M
    
    def apply_boundary_conditions(self, K, M, fixed_dofs=[0, 1]):
        """应用边界条件"""
        ndof = K.shape[0]
        all_dofs = np.arange(ndof)
        free_dofs = np.setdiff1d(all_dofs, fixed_dofs)
        
        K_reduced = K[np.ix_(free_dofs, free_dofs)]
        M_reduced = M[np.ix_(free_dofs, free_dofs)]
        
        return K_reduced, M_reduced, free_dofs
    
    def solve_eigenvalue_problem(self, K, M, n_modes=5):
        """求解广义特征值问题"""
        eigenvalues, eigenvectors = linalg.eigh(K, M)
        frequencies = np.sqrt(np.abs(eigenvalues)) / (2 * np.pi)
        frequencies = frequencies[:n_modes]
        modes = eigenvectors[:, :n_modes]
        
        # 质量归一化
        for i in range(n_modes):
            mass_norm = np.sqrt(modes[:, i].T @ M @ modes[:, i])
            modes[:, i] /= mass_norm
        
        return frequencies, modes

def case3_stochastic_fem():
    """案例3：随机有限元方法"""
    print("="*60)
    print("主题064 - 案例3: 随机有限元方法")
    print("="*60)
    
    # 定义有限元模型
    sfem = StochasticFEM(
        n_elements=20,
        length=10.0,
        E_mean=200e9,
        rho_mean=7850.0,
        A=0.01,
        I=8.33e-6
    )
    
    # 定义随机场
    E_random_field = RandomField(
        mean=200e9,
        std=0.15,
        correlation_length=2.0,
        domain_length=10.0
    )
    
    rho_random_field = RandomField(
        mean=7850.0,
        std=0.10,
        correlation_length=3.0,
        domain_length=10.0
    )
    
    # 计算KL展开
    n_modes_kl = 8
    E_eigenvalues, E_eigenvectors, x_kl = E_random_field.compute_kle_modes(
        n_modes_kl, sfem.n_nodes)
    rho_eigenvalues, rho_eigenvectors, _ = rho_random_field.compute_kle_modes(
        n_modes_kl, sfem.n_nodes)
    
    # 确定性分析
    E_field = np.ones(sfem.n_nodes) * sfem.E_mean
    rho_field = np.ones(sfem.n_nodes) * sfem.rho_mean
    K, M = sfem.assemble_system_matrices(E_field, rho_field)
    K_reduced, M_reduced, _ = sfem.apply_boundary_conditions(K, M)
    det_freqs, _ = sfem.solve_eigenvalue_problem(K_reduced, M_reduced, 5)
    
    print("\n确定性分析结果:")
    for i, f in enumerate(det_freqs):
        print(f"  第{i+1}阶频率 = {f:.4f} Hz")
    
    # 蒙特卡洛模拟
    print("\n随机有限元分析...")
    n_samples = 500
    frequencies_samples = []
    
    for i in range(n_samples):
        if (i + 1) % 100 == 0:
            print(f"  进度: {i+1}/{n_samples}")
        
        # 生成随机变量
        xi_E = np.random.randn(n_modes_kl)
        xi_rho = np.random.randn(n_modes_kl)
        
        # 生成随机场
        E_field = E_random_field.generate_realization(
            xi_E, E_eigenvalues, E_eigenvectors, x_kl)
        rho_field = rho_random_field.generate_realization(
            xi_rho, rho_eigenvalues, rho_eigenvectors, x_kl)
        
        # 求解
        K, M = sfem.assemble_system_matrices(E_field, rho_field)
        K_reduced, M_reduced, _ = sfem.apply_boundary_conditions(K, M)
        frequencies, _ = sfem.solve_eigenvalue_problem(K_reduced, M_reduced, 5)
        
        frequencies_samples.append(frequencies)
    
    frequencies_samples = np.array(frequencies_samples)
    
    # 统计结果
    print("\n随机分析结果:")
    print(f"{'模态':<8} {'均值(Hz)':<15} {'标准差(Hz)':<15} {'变异系数(%)':<15}")
    print("-" * 55)
    for i in range(5):
        mean = np.mean(frequencies_samples[:, i])
        std = np.std(frequencies_samples[:, i])
        cov = std / mean * 100
        print(f"{i+1:<8} {mean:<15.4f} {std:<15.4f} {cov:<15.2f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 子图1：频率分布
    ax = axes[0, 0]
    for i in range(3):
        ax.hist(frequencies_samples[:, i], bins=30, alpha=0.5, 
                label=f'Mode {i+1}')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Count')
    ax.set_title('Distribution of Natural Frequencies')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 子图2：频率统计特性
    ax = axes[0, 1]
    mode_numbers = np.arange(1, 6)
    means = np.mean(frequencies_samples, axis=0)
    stds = np.std(frequencies_samples, axis=0)
    
    ax.bar(mode_numbers, means, yerr=2*stds, capsize=5, alpha=0.7)
    ax.set_xlabel('Mode Number')
    ax.set_ylabel('Frequency (Hz)')
    ax.set_title('Mean Natural Frequencies with 95% CI')
    ax.grid(True, alpha=0.3, axis='y')
    
    # 子图3：变异系数
    ax = axes[1, 0]
    cvs = stds / means * 100
    ax.bar(mode_numbers, cvs, alpha=0.7, color='coral')
    ax.set_xlabel('Mode Number')
    ax.set_ylabel('Coefficient of Variation (%)')
    ax.set_title('Frequency Variability')
    ax.grid(True, alpha=0.3, axis='y')
    
    # 子图4：确定性与随机对比
    ax = axes[1, 1]
    x_pos = np.arange(5)
    width = 0.35
    ax.bar(x_pos - width/2, det_freqs, width, label='Deterministic')
    ax.bar(x_pos + width/2, means, width, yerr=stds, capsize=3, label='Stochastic')
    ax.set_xlabel('Mode Number')
    ax.set_ylabel('Frequency (Hz)')
    ax.set_title('Deterministic vs Stochastic')
    ax.set_xticks(x_pos)
    ax.set_xticklabels([f'{i+1}' for i in range(5)])
    ax.legend()
    ax.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('sfem_results.png', dpi=150, bbox_inches='tight')
    print("\n随机有限元结果图已保存: sfem_results.png")
    plt.close()

if __name__ == "__main__":
    case3_stochastic_fem()

案例3结果分析：

随机有限元方法展示了以下特点：

1. 随机场建模：KL展开有效地将连续随机场离散为有限个随机变量，大幅降低了问题维度
2. 频率不确定性：各阶固有频率都表现出显著的不确定性，变异系数在7-10%之间
3. 模态相关性：高阶模态对材料参数的不确定性更加敏感
4. 工程应用：SFEM为复杂结构的可靠性分析提供了有效工具

6.4 案例4：可靠性分析方法

问题描述：
使用可靠性分析方法评估悬臂梁在随机载荷作用下的失效概率。考虑应力和挠度两种极限状态，使用FORM、蒙特卡洛和重要性采样方法计算失效概率。

理论基础：

极限状态函数：

应力极限状态：$$g_{\sigma }(X)={\sigma }_{allowable}-{\sigma }_{max}(X)$$

挠度极限状态：$$g_{\delta }(X)={\delta }_{allowable}-{\delta }_{max}(X)$$

失效概率：

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

FORM方法：

通过寻找设计点（最可能失效点）计算可靠性指标：

$$\beta =\underset{g(u)=0}{\min }\sqrt{u^{T}u}$$

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

Python代码实现：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm, gumbel_r

class RandomVariable:
    """随机变量类"""
    
    def __init__(self, name, distribution, params):
        self.name = name
        self.distribution = distribution
        self.params = params
        
        if distribution == 'normal':
            self.mean = params['mean']
            self.std = params['std']
            self.dist = norm(loc=self.mean, scale=self.std)
        elif distribution == 'lognormal':
            self.mean = params['mean']
            self.std = params['std']
            sigma_ln = np.sqrt(np.log(1 + (self.std/self.mean)**2))
            mu_ln = np.log(self.mean) - 0.5 * sigma_ln**2
            self.dist = lognorm(s=sigma_ln, scale=np.exp(mu_ln))
        elif distribution == 'gumbel':
            self.mean = params['mean']
            self.std = params['std']
            beta = self.std * np.sqrt(6) / np.pi
            mu = self.mean - 0.5772 * beta
            self.dist = gumbel_r(loc=mu, scale=beta)
    
    def pdf(self, x):
        return self.dist.pdf(x)
    
    def cdf(self, x):
        return self.dist.cdf(x)
    
    def sample(self, size=1):
        return self.dist.rvs(size=size)
    
    def transform_to_standard(self, x):
        u = self.dist.cdf(x)
        u = np.clip(u, 1e-10, 1 - 1e-10)
        return norm.ppf(u)
    
    def transform_from_standard(self, u):
        p = norm.cdf(u)
        return self.dist.ppf(p)

class FORM:
    """一阶可靠性方法"""
    
    def __init__(self, limit_state_func, random_vars):
        self.g = limit_state_func
        self.vars = random_vars
        self.n_vars = len(random_vars)
    
    def transform_to_standard(self, X):
        U = np.zeros(self.n_vars)
        for i, var in enumerate(self.vars):
            U[i] = var.transform_to_standard(X[i])
        return U
    
    def transform_from_standard(self, U):
        X = np.zeros(self.n_vars)
        for i, var in enumerate(self.vars):
            X[i] = var.transform_from_standard(U[i])
        return X
    
    def limit_state_in_standard(self, U):
        X = self.transform_from_standard(U)
        return self.g(X)
    
    def compute_reliability(self, tol=1e-6, max_iter=100):
        """计算可靠性指标"""
        X0 = np.array([var.mean for var in self.vars])
        U = self.transform_to_standard(X0)
        
        for iteration in range(max_iter):
            g_val = self.limit_state_in_standard(U)
            grad = self._numerical_gradient(U)
            grad_norm = np.linalg.norm(grad)
            
            if grad_norm < 1e-10:
                break
            
            alpha = grad / grad_norm
            beta_new = -g_val / grad_norm + np.dot(U, alpha)
            U_new = -beta_new * alpha
            
            if np.linalg.norm(U_new - U) < tol:
                U = U_new
                break
            
            U = U_new
        
        beta = np.linalg.norm(U)
        pf = norm.cdf(-beta)
        
        return {
            'beta': beta,
            'pf': pf,
            'design_point_X': self.transform_from_standard(U),
            'design_point_U': U
        }
    
    def _numerical_gradient(self, U, h=1e-6):
        grad = np.zeros(self.n_vars)
        g0 = self.limit_state_in_standard(U)
        
        for i in range(self.n_vars):
            U_plus = U.copy()
            U_plus[i] += h
            g_plus = self.limit_state_in_standard(U_plus)
            grad[i] = (g_plus - g0) / h
        
        return grad

def cantilever_stress(P, L, b, h, E):
    """悬臂梁最大应力"""
    I = b * h**3 / 12
    M_max = P * L
    stress = M_max * h / (2 * I)
    return stress

def case4_reliability():
    """案例4：可靠性分析"""
    print("="*60)
    print("主题064 - 案例4: 可靠性分析方法")
    print("="*60)
    
    # 定义随机变量
    random_vars = [
        RandomVariable('P', 'gumbel', {'mean': 5000, 'std': 1500}),
        RandomVariable('L', 'normal', {'mean': 5.0, 'std': 0.05}),
        RandomVariable('b', 'normal', {'mean': 0.2, 'std': 0.01}),
        RandomVariable('h', 'lognormal', {'mean': 0.3, 'std': 0.015}),
        RandomVariable('E', 'normal', {'mean': 200e9, 'std': 10e9})
    ]
    
    allowable_stress = 250e6
    
    # 定义极限状态函数
    def limit_state_stress(X):
        P, L, b, h, E = X
        stress = cantilever_stress(P, L, b, h, E)
        return allowable_stress - stress
    
    # FORM分析
    print("\nFORM分析...")
    form = FORM(limit_state_stress, random_vars)
    results = form.compute_reliability()
    
    print(f"\n应力极限状态分析结果:")
    print(f"  可靠性指标 β = {results['beta']:.4f}")
    print(f"  失效概率 Pf = {results['pf']:.4e}")
    print(f"\n设计点:")
    for i, var in enumerate(random_vars):
        print(f"  {var.name} = {results['design_point_X'][i]:.6f}")
    
    # 蒙特卡洛验证
    print("\n蒙特卡洛分析...")
    n_samples = 100000
    n_failures = 0
    
    for i in range(n_samples):
        X = np.array([var.sample() for var in random_vars])
        if limit_state_stress(X) < 0:
            n_failures += 1
    
    pf_mc = n_failures / n_samples
    beta_mc = -norm.ppf(pf_mc) if pf_mc > 0 else np.inf
    
    print(f"  失效概率 Pf = {pf_mc:.4e}")
    print(f"  可靠性指标 β = {beta_mc:.4f}")
    
    # 方法对比
    print("\n" + "-"*60)
    print("方法对比:")
    print(f"{'方法':<20} {'β':<15} {'Pf':<15}")
    print("-"*60)
    print(f"{'FORM':<20} {results['beta']:<15.4f} {results['pf']:<15.4e}")
    print(f"{'Monte Carlo':<20} {beta_mc:<15.4f} {pf_mc:<15.4e}")
    print("-"*60)

if __name__ == "__main__":
    case4_reliability()

案例4结果分析：

可靠性分析展示了以下特点：

1. FORM效率：FORM方法仅需少量迭代即可得到失效概率估计，计算效率远高于蒙特卡洛方法
2. 设计点信息：设计点指示了最可能导致失效的参数组合，为结构改进提供指导
3. 方法对比：对于小失效概率问题（Pf < 0.01），蒙特卡洛需要大量样本，而FORM保持高效
4. 工程决策：可靠性指标 $\beta$ 可直接用于结构设计规范和安全评估

---

7. 小结

本主题系统介绍了结构动力学中的不确定性量化方法，包括：

7.1 方法对比

方法	优点	缺点	适用场景
蒙特卡洛模拟	直观、通用、精度可控	计算成本高	验证其他方法、复杂非线性问题
多项式混沌展开	计算效率高、收敛快	需要光滑响应	光滑响应函数、参数不确定性
随机有限元方法	与FEM结合、适合复杂结构	实现复杂	复杂结构、空间随机场
可靠性分析	直接给出失效概率	对非线性敏感	安全评估、设计优化

7.2 关键要点

1. 不确定性表征：概率方法是工程中最常用的不确定性表征方式
2. 方差缩减：拉丁超立方采样、重要性采样等技术可显著提高蒙特卡洛效率
3. 代理模型：PCE等代理模型方法在保持精度的同时大幅降低计算成本
4. 工程应用：不确定性量化为结构安全评估和优化设计提供了更全面的视角

7.3 发展趋势

• 高维不确定性：发展适用于高维问题的UQ方法
• 模型不确定性：考虑模型形式误差和参数识别误差
• 多尺度分析：跨尺度的不确定性传播方法
• 机器学习结合：利用深度学习加速不确定性量化

不确定性量化为结构设计和安全评估提供了更全面的视角，是现代工程分析不可或缺的一部分。掌握这些方法，工程师可以更好地处理实际工程中的不确定性因素，设计出更安全、更可靠的结构系统。

---