主题045：不确定性量化与可靠性分析

1. 引言

1.1 工程中的不确定性

在实际工程问题中，不确定性无处不在：

材料参数的不确定性

• 材料属性的分散性（弹性模量、密度、热导率等）
• 制造公差导致的几何尺寸变化
• 材料缺陷和微观结构不均匀性

载荷条件的不确定性

• 环境载荷的随机性（风载、地震、波浪等）
• 使用载荷的变异性
• 边界条件的近似性

模型不确定性

• 数学模型的简化假设
• 数值离散误差
• 模型参数识别误差
  

1.2 不确定性量化的重要性

设计决策支持

• 评估设计的安全裕度
• 优化设计参数
• 权衡性能与可靠性

风险评估

• 计算失效概率
• 识别关键失效模式
• 制定风险缓解策略

认证与合规

• 满足安全标准
• 提供统计置信度
• 支持监管决策

1.3 应用领域

• 航空航天：飞行器结构可靠性、气动不确定性
• 土木工程：桥梁安全评估、地震风险分析
• 核工程：核反应堆安全分析、放射性物质扩散
• 机械工程：疲劳寿命预测、磨损分析
• 电子工程：电路参数容差分析、热可靠性

---

2. 不确定性来源与分类

2.1 不确定性的类型

2.1.1 偶然不确定性（Aleatory Uncertainty）

定义：固有的、不可减少的随机性

特征：

• 源于物理系统的内在随机性
• 无法通过增加数据消除
• 可用概率模型描述

例子：

• 材料微观结构的随机分布
• 环境载荷的随机波动
• 制造过程中的随机误差

2.1.2 认知不确定性（Epistemic Uncertainty）

定义：由于知识不足导致的不确定性

特征：

• 可通过增加数据或改进模型减少
• 反映信息的不完整性
• 可用区间分析或证据理论描述

例子：

• 模型参数估计误差
• 数学模型的简化假设
• 缺乏足够实验数据

2.2 不确定性表征方法

2.2.1 概率方法

概率密度函数（PDF）

$$p_X(x)\ge 0,{\int }_{-\infty }^{\infty }{p}_X(x)dx=1$$

累积分布函数（CDF）

$$F_X(x)=P(X\le x)={\int }_{-\infty }^x{p}_X(t)dt$$

常用分布类型

分布类型	概率密度函数	应用场景
正态分布	$p(x)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$	测量误差、材料属性
对数正态分布	$p(x)=\frac{1}{x\sigma \sqrt{2\pi }}{e}^{-\frac{(\ln x-\mu {)}^2}{2{\sigma }^2}}$	疲劳寿命、颗粒尺寸
均匀分布	$p(x)=\frac{1}{b-a},a\le x\le b$	公差范围、设计参数
威布尔分布	$p(x)=\frac{k}{\lambda }(\frac{x}{\lambda }{)}^{k-1}{e}^{-(x/\lambda {)}^k}$	强度、可靠性分析
指数分布	$p(x)=\lambda e^{-\lambda x}$	泊松过程、等待时间

2.2.2 区间方法

区间变量

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

区间运算

$$[a,b]+[c,d]=[a+c,b+d]$$
$$[a,b]\times [c,d]=[\min (ac,ad,bc,bd),\max (ac,ad,bc,bd)]$$

应用场景：

• 认知不确定性主导时
• 数据稀缺情况
• 保守安全评估

2.2.3 模糊方法

隶属度函数

$${\mu }_A(x)\in [0,1]$$

模糊运算

• 模糊加法、乘法
• 模糊微分方程
• 模糊优化

应用场景：

• 专家知识量化
• 语言变量描述
• 主观不确定性

2.3 不确定性建模流程

┌─────────────────┐
│  识别不确定性源  │
└────────┬────────┘
         ↓
┌─────────────────┐
│  选择表征方法   │
│ (概率/区间/模糊) │
└────────┬────────┘
         ↓
┌─────────────────┐
│  参数估计与验证  │
└────────┬────────┘
         ↓
┌─────────────────┐
│  不确定性传播   │
└────────┬────────┘
         ↓
┌─────────────────┐
│  结果分析与决策  │
└─────────────────┘

---

3. 概率统计基础

3.1 随机变量与分布

3.1.1 数字特征

期望值（均值）

$${\mu }_X=E[X]={\int }_{-\infty }^{\infty }x{p}_X(x)dx$$

方差

$${\sigma }_X^2=Var[X]=E[(X-{\mu }_X{)}^2]={\int }_{-\infty }^{\infty }(x-{\mu }_X{)}^2{p}_X(x)dx$$

标准差

$${\sigma }_X=\sqrt{Var[X]}$$

变异系数

$${\delta }_X=\frac{{\sigma }_X}{{\mu }_X}$$

3.1.2 多元随机变量

联合概率密度

$$p_{XY}(x,y)=p_X(x)p_{Y|X}(y|x)$$

协方差

$$Cov[X,Y]=E[(X-{\mu }_X)(Y-{\mu }_Y)]$$

相关系数

$${\rho }_{XY}=\frac{Cov[X,Y]}{{\sigma }_X{\sigma }_Y}$$

协方差矩阵

$$\mathbf{C}=\left[\begin{array}{ccc}{\sigma }_{X_1}^2& {\rho }_{12}{\sigma }_{X_1}{\sigma }_{X_2}& \cdots \\ {\rho }_{21}{\sigma }_{X_2}{\sigma }_{X_1}& {\sigma }_{X_2}^2& \cdots \\ ⋮& ⋮& \ddots \end{array}\right]$$

3.2 统计推断

3.2.1 参数估计

矩估计法

用样本矩估计总体矩：

$$\hat{\mu }=\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i$$
$${\hat{\sigma }}^2=s^2=\frac{1}{n-1}\sum _{i=1}^n(x_i-\bar{x}{)}^2$$

最大似然估计

$$\hat{\theta }=\arg \underset{\theta }{\max }L(\theta ;x_1,...,x_n)=\arg \underset{\theta }{\max }\prod _{i=1}^{n}p(x_i;\theta )$$

贝叶斯估计

$$p(\theta |D)=\frac{p(D|\theta )p(\theta )}{p(D)}$$

3.2.2 假设检验

正态性检验

• Kolmogorov-Smirnov检验
• Shapiro-Wilk检验
• Anderson-Darling检验

拟合优度检验

• 卡方检验
• 似然比检验

3.3 随机过程

3.3.1 平稳随机过程

定义：统计特性不随时间平移改变

自相关函数

$$R_X(\tau )=E[X(t)X(t+\tau )]$$

功率谱密度

$$S_X(\omega )={\int }_{-\infty }^{\infty }{R}_X(\tau )e^{-j\omega \tau }d\tau$$

3.3.2 高斯随机过程

性质：

• 任意有限维分布都是联合高斯分布
• 完全由均值函数和协方差函数确定
• 线性变换保持高斯性

Karhunen-Loève展开

$$X(t)=\mu (t)+\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}{\xi }_i{\varphi }_i(t)$$

其中 ${\xi }_i$ 是标准正态随机变量，${\lambda }_i$ 和 ${\varphi }_i(t)$ 是协方差核的特征值和特征函数。

---

4. 不确定性传播方法

4.1 蒙特卡洛模拟

4.1.1 基本方法

算法流程

def monte_carlo_simulation(n_samples, input_distributions, model):
    """
    蒙特卡洛模拟
    
    参数:
        n_samples: 样本数量
        input_distributions: 输入变量分布列表
        model: 仿真模型函数
    
    返回:
        output_samples: 输出样本
    """
    output_samples = []
    
    for i in range(n_samples):
        # 1. 从输入分布采样
        x = [dist.rvs() for dist in input_distributions]
        
        # 2. 运行模型
        y = model(x)
        
        # 3. 保存输出
        output_samples.append(y)
    
    return np.array(output_samples)

统计估计

均值估计：
$${\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$$

失效概率估计：
$${\hat{P}}_f=\frac{1}{N}\sum _{i=1}^{N}I(g(x_i)\le 0)$$

其中 $I(\cdot )$ 是指示函数。

4.1.2 方差缩减技术

重要性采样

$$P_f=\int I(g(x)\le 0)\frac{p_X(x)}{q_X(x)}{q}_X(x)dx$$

选择重要性密度 $q_X(x)$ 在失效域附近有更多采样。

分层采样

将输入空间划分为若干层，每层独立采样：

$${\hat{\mu }}_Y=\sum _{j=1}^L{w}_j{\bar{y}}_j$$

拉丁超立方采样（LHS）

确保样本在输入空间中均匀分布，减少聚类。

def latin_hypercube_sampling(n_samples, n_dims):
    """拉丁超立方采样"""
    # 生成排列
    samples = np.zeros((n_samples, n_dims))
    
    for i in range(n_dims):
        perm = np.random.permutation(n_samples)
        samples[:, i] = (perm + np.random.rand(n_samples)) / n_samples
    
    return samples

拟蒙特卡洛方法

使用低差异序列（如Sobol序列）代替伪随机数：

from scipy.stats import qmc

def sobol_sampling(n_samples, n_dims):
    """Sobol序列采样"""
    sampler = qmc.Sobol(d=n_dims, scramble=True)
    samples = sampler.random(n=n_samples)
    return samples

4.2 谱方法

4.2.1 多项式混沌展开（PCE）

Wiener-Askey方案

将随机输出展开为随机变量的多项式：

$$Y(\xi )=\sum _{\alpha \in \mathcal{A}}{y}_{\alpha }{\Psi }_{\alpha }(\xi )$$

其中：

• $\xi$：标准随机变量（高斯、均匀等）
• ${\Psi }_{\alpha }$：正交多项式（Hermite、Legendre等）
• $y_{\alpha }$：展开系数

正交多项式族

输入分布	多项式类型	正交性
高斯分布	Hermite	$\int H_m(x)H_n(x)e^{-x^2/2}dx=\sqrt{2\pi }n!{\delta }_{mn}$
均匀分布	Legendre	${\int }_{-1}^1{P}_m(x)P_n(x)dx=\frac{2}{2n+1}{\delta }_{mn}$
Gamma分布	Laguerre	${\int }_0^{\infty }{L}_m(x)L_n(x)e^{-x}dx={\delta }_{mn}$
Beta分布	Jacobi	${\int }_{-1}^1(1-x{)}^{\alpha }(1+x{)}^{\beta }{J}_m(x)J_n(x)dx={\delta }_{mn}$

系数计算

投影法：
$$y_{\alpha }=\frac{⟨Y,{\Psi }_{\alpha }⟩}{⟨{\Psi }_{\alpha },{\Psi }_{\alpha }⟩}=\frac{E[Y{\Psi }_{\alpha }]}{E[{\Psi }_{\alpha }^2]}$$

回归法：
$$\underset{y_{\alpha }}{\min }\sum _{i=1}^N{\left(y_i-\sum _{\alpha }{y}_{\alpha }{\Psi }_{\alpha }({\xi }_i)\right)}^2$$

统计矩计算

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

方差：${\sigma }_Y^2={\sum }_{\alpha \ne 0}{y}_{\alpha }^{2}E[{\Psi }_{\alpha }^2]$

4.2.2 随机配点法

基本思想

在配点上精确满足方程，通过插值获得全局近似。

配点选择

• 全张量积配点
• 稀疏网格配点（Smolyak算法）

Smolyak稀疏网格

$$A(q,d)=\sum _{q-d+1\le |i|\le q}(-1{)}^{q-|i|}(\genfrac{}{}{0ex}{}{d-1}{q-|i|})(U^{i_1}\otimes \cdots \otimes U^{i_d})$$

4.3 矩方法

4.3.1 一阶矩方法（FOSM）

泰勒展开

$$Y\approx g({\mu }_X)+\sum _{i=1}^n\frac{\partial g}{\partial X_i}{|}_{{\mu }_X}(X_i-{\mu }_{X_i})$$

均值估计

$${\mu }_Y\approx g({\mu }_X)$$

方差估计

$${\sigma }_Y^2\approx \sum _{i=1}^n\sum _{j=1}^n\frac{\partial g}{\partial X_i}\frac{\partial g}{\partial X_j}Cov[X_i,X_j]$$

对于独立变量：
$${\sigma }_Y^2\approx \sum _{i=1}^n{\left(\frac{\partial g}{\partial X_i}\right)}^2{\sigma }_{X_i}^2$$

4.3.2 二阶矩方法（SOSM）

包含二阶项：

$$Y\approx g({\mu }_X)+\sum _{i=1}^n\frac{\partial g}{\partial X_i}(X_i-{\mu }_{X_i})+\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n\frac{{\partial }^{2}g}{\partial X_i\partial X_j}(X_i-{\mu }_{X_i})(X_j-{\mu }_{X_j})$$

均值修正

$${\mu }_Y\approx g({\mu }_X)+\frac{1}{2}\sum _{i=1}^n\frac{{\partial }^{2}g}{\partial X_i^2}{\sigma }_{X_i}^2$$

4.4 代理模型方法

4.4.1 高斯过程回归

模型

$$Y(x)=f(x{)}^T\beta +Z(x)$$

其中 $Z(x)$ 是零均值高斯过程：

$$Cov[Z(x),Z(x^{\prime })]={\sigma }^{2}R(x,x^{\prime })$$

常用核函数

高斯核：
$$R(x,x^{\prime })=\exp \left(-\sum _{i=1}^d{\theta }_i(x_i-x_i^{\prime }{)}^2\right)$$

Matérn核：
$$R(x,x^{\prime })=\frac{1}{\Gamma (\nu )2^{\nu -1}}{\left(\frac{2\sqrt{\nu }|x-x^{\prime }|}{l}\right)}^{\nu }{K}_{\nu }\left(\frac{2\sqrt{\nu }|x-x^{\prime }|}{l}\right)$$

4.4.2 多项式响应面

一阶模型

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

二阶模型

$$\hat{Y}={\beta }_0+\sum _{i=1}^n{\beta }_i{X}_i+\sum _{i=1}^n{\beta }_{ii}{X}_i^2+\sum _{i<j}{\beta }_{ij}{X}_i{X}_j$$

---

5. 可靠性分析方法

5.1 失效概率计算

5.1.1 极限状态函数

定义

$$g(X)=\{\begin{array}{ll}>0& \text{安全}\\ =0& \text{极限状态}\\ <0& \text{失效}\end{array}$$

失效概率

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

可靠度指标

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

其中 ${\Phi }^{-1}$ 是标准正态CDF的逆函数。

5.1.2 蒙特卡洛可靠性分析

直接模拟

def mc_reliability_analysis(n_samples, input_dist, limit_state):
    """
    蒙特卡洛可靠性分析
    
    参数:
        n_samples: 样本数
        input_dist: 输入分布
        limit_state: 极限状态函数
    
    返回:
        pf: 失效概率估计
        cov: 变异系数
    """
    n_fail = 0
    
    for i in range(n_samples):
        x = input_dist.rvs()
        if limit_state(x) <= 0:
            n_fail += 1
    
    pf = n_fail / n_samples
    cov = np.sqrt((1 - pf) / (n_samples * pf))
    
    return pf, cov

样本量估计

对于目标变异系数 $CO{V}_{target}$：

$$N=\frac{1-P_f}{P_f\cdot CO{V}_{target}^2}$$

对于小失效概率（如 $P_f=10^{-6}$），需要大量样本。

5.2 一阶可靠性方法（FORM）

5.2.1 标准正态空间变换

Rosenblatt变换

对于独立随机变量：
$$u_i={\Phi }^{-1}(F_{X_i}(x_i))$$

Nataf变换

对于相关随机变量：
$$u=L^{-1}z$$

其中 $L$ 是相关矩阵的Cholesky分解，$z_i={\Phi }^{-1}(F_{X_i}(x_i))$。

5.2.2 最可能失效点（MPP）

优化问题

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

HL-RF算法

def hl_rf_algorithm(u0, g, dg, tol=1e-6, max_iter=100):
    """
    HL-RF算法求解MPP
    
    参数:
        u0: 初始点
        g: 极限状态函数
        dg: 梯度函数
        tol: 收敛容差
        max_iter: 最大迭代次数
    
    返回:
        u_star: MPP点
        beta: 可靠度指标
    """
    u = u0
    
    for i in range(max_iter):
        # 计算函数值和梯度
        g_val = g(u)
        grad = dg(u)
        grad_norm = np.linalg.norm(grad)
        
        # 更新方向
        alpha = grad / grad_norm
        
        # 更新MPP估计
        u_new = -alpha * (g_val / grad_norm - np.dot(alpha, u))
        
        # 检查收敛
        if np.linalg.norm(u_new - u) < tol:
            break
        
        u = u_new
    
    beta = np.linalg.norm(u)
    return u, beta

5.2.3 FORM近似

失效概率估计

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

误差来源

• 极限状态函数的非线性
• 高阶效应忽略

5.3 二阶可靠性方法（SORM）

5.3.1 曲率近似

Breitung公式

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

其中 ${\kappa }_i$ 是极限状态函数在MPP点的主曲率。

曲率计算

通过Hessian矩阵的特征值分解获得主曲率。

5.3.2 点拟合SORM

在MPP点周围拟合二次曲面：

$$g(u)\approx \beta -{\alpha }^{T}u+\frac{1}{2}(u-u^{\ast }{)}^{T}H(u-u^{\ast })$$

5.4 重要性抽样

5.4.1 最优重要性密度

$$q^{\ast }(x)=\frac{I(g(x)\le 0)p_X(x)}{P_f}$$

由于 $P_f$ 未知，使用近似：

$$q(x)\approx p_X(x|g(x)\le 0)$$

5.4.2 基于FORM的重要性抽样

以MPP为中心构建重要性密度：

$$q(x)={\varphi }_n(u-u^{\ast })\left|\frac{\partial u}{\partial x}\right|$$

失效概率估计

$${\hat{P}}_f=\frac{1}{N}\sum _{i=1}^{N}I(g(x_i)\le 0)\frac{p_X(x_i)}{q(x_i)}$$

---

6. 敏感性分析

6.1 局部敏感性分析

6.1.1 偏导数法

一阶敏感性系数

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

标准化敏感性系数

$$S_i^{std}=\frac{{\sigma }_{X_i}}{{\sigma }_Y}\frac{\partial Y}{\partial X_i}$$

6.1.2 基于方差的敏感性

FORM敏感性因子

$${\alpha }_i=-\frac{\partial \beta }{\partial u_i}=\frac{\partial g/\partial u_i}{\Vert {\nabla }_{u}g\Vert }$$

${\alpha }_i^2$ 表示第 $i$ 个变量对总方差的贡献比例。

6.2 全局敏感性分析

6.2.1 Sobol敏感性指数

方差分解

$$Var[Y]=\sum _i{V}_i+\sum _{i<j}{V}_{ij}+\cdots +V_{12\cdots n}$$

其中：
$$V_i=Var[E[Y|X_i]]$$
$$V_{ij}=Var[E[Y|X_i,X_j]]-V_i-V_j$$

一阶Sobol指数

$$S_i=\frac{V_i}{Var[Y]}$$

总效应指数

$$S_i^T=\frac{E[Var[Y|X_{\sim i}]]}{Var[Y]}=1-\frac{Var[E[Y|X_{\sim i}]]}{Var[Y]}$$

蒙特卡洛估计

def sobol_indices_mc(n_samples, model, input_dists):
    """
    基于蒙特卡洛的Sobol指数估计
    
    使用Saltelli采样方案
    """
    n_inputs = len(input_dists)
    
    # 生成两个独立样本矩阵
    A = np.array([dist.rvs(n_samples) for dist in input_dists]).T
    B = np.array([dist.rvs(n_samples) for dist in input_dists]).T
    
    # 计算模型输出
    Y_A = np.array([model(a) for a in A])
    Y_B = np.array([model(b) for b in B])
    
    # 总方差
    V_Y = np.var(Y_A)
    
    S1 = np.zeros(n_inputs)
    ST = np.zeros(n_inputs)
    
    for i in range(n_inputs):
        # 创建混合矩阵
        A_B = A.copy()
        A_B[:, i] = B[:, i]
        
        Y_AB = np.array([model(a_b) for a_b in A_B])
        
        # 一阶指数
        S1[i] = np.mean(Y_B * (Y_AB - Y_A)) / V_Y
        
        # 总效应指数
        ST[i] = 0.5 * np.mean((Y_A - Y_AB)**2) / V_Y
    
    return S1, ST

6.2.2 Morris方法

基本思想

通过计算"基本效应"来筛选重要参数：

$$E{E}_i=\frac{Y(x_1,...,x_i+\Delta ,...,x_n)-Y(x)}{\Delta }$$

统计量

• $\mu$：基本效应的均值（总体重要性）
• $\sigma$：基本效应的标准差（非线性/交互作用）

采样策略

在网格上构建轨迹，每个轨迹在 $n+1$ 个点上评估模型。

6.3 基于PCE的敏感性分析

** Sobol指数解析计算**

利用PCE的正交性：

$$S_i=\frac{{\sum }_{\alpha \in {\mathcal{A}}_i}{y}_{\alpha }^{2}E[{\Psi }_{\alpha }^2]}{{\sum }_{\alpha \ne 0}{y}_{\alpha }^{2}E[{\Psi }_{\alpha }^2]}$$

其中 Ai={α∣αi>0,αj≠i=0}\mathcal{A}_i = \{\alpha | \alpha_i > 0, \alpha_{j \neq i} = 0\}Ai​={α∣αi​>0,αj=i​=0}。

总效应指数

$$S_i^T=\frac{{\sum }_{\alpha \in {\mathcal{A}}_i^T}{y}_{\alpha }^{2}E[{\Psi }_{\alpha }^2]}{{\sum }_{\alpha \ne 0}{y}_{\alpha }^{2}E[{\Psi }_{\alpha }^2]}$$

其中 AiT={α∣αi>0}\mathcal{A}_i^T = \{\alpha | \alpha_i > 0\}AiT​={α∣αi​>0}。

---

7. 多物理场问题的不确定性量化

7.1 耦合系统的不确定性传播

7.1.1 分区耦合方法

松耦合

每个物理场独立进行不确定性量化，通过接口传递统计信息：

$$Y_1=M_1(X_1,Y_2),Y_2=M_2(X_2,Y_1)$$

紧耦合

同时考虑所有物理场的不确定性：

$$Y=M(X),X=[X_1,X_2,...,X_n]$$

7.1.2 多保真度方法

模型层次

• 高保真模型：CFD/FEA详细仿真
• 低保真模型：降阶模型、代理模型

自适应采样

def adaptive_multi_fidelity_uq(budget, hf_model, lf_model):
    """
    自适应多保真度UQ
    
    策略：
    1. 大量低保真样本
    2. 少量高保真样本修正
    3. 自适应分配计算资源
    """
    # 低保真模型建立响应面
    lf_samples = ...
    lf_responses = [lf_model(x) for x in lf_samples]
    
    # 高保真样本用于修正
    hf_samples = ...
    hf_responses = [hf_model(x) for x in hf_samples]
    
    # 构建修正模型
    correction = hf_responses - [lf_model(x) for x in hf_samples]
    
    # 组合预测
    return lf_prediction + correction_prediction

7.2 时变可靠性分析

7.2.1 首次穿越问题

定义

计算随机过程首次超过安全边界的概率：

$$P_f(0,T)=P(\exists t\in [0,T]:g(X(t))\le 0)$$

解析方法

对于高斯过程，使用Rice公式：

$${\nu }^{+}={\int }_0^{\infty }\dot{y}{p}_{Y\dot{Y}}(b,\dot{y})d\dot{y}$$

蒙特卡洛方法

def outcrossing_rate_mc(n_samples, process_generator, limit_state, dt):
    """
    穿越率蒙特卡洛估计
    """
    crossings = 0
    
    for _ in range(n_samples):
        # 生成样本路径
        trajectory = process_generator()
        
        # 检测穿越
        for i in range(len(trajectory)-1):
            if limit_state(trajectory[i]) > 0 and limit_state(trajectory[i+1]) <= 0:
                crossings += 1
                break
    
    return crossings / (n_samples * dt)

7.2.2 极值分布方法

极值统计

关注时间区间内的最大响应：

$$Y_{max}=\underset{t\in [0,T]}{\max }Y(t)$$

极值分布近似

使用Gumbel分布拟合：

$$F_{Y_{max}}(y)=\exp \left(-\exp \left(-\frac{y-\mu }{\sigma }\right)\right)$$

7.3 系统可靠性

7.3.1 系统模型

串联系统

$$P_f^{sys}=1-\prod _{i=1}^n(1-P_{f,i})$$

并联系统

$$P_f^{sys}=\prod _{i=1}^n{P}_{f,i}$$

混联系统

组合串联和并联结构。

7.3.2 故障树分析

基本事件

• 底事件：基本失效模式
• 顶事件：系统失效

逻辑门

• 与门：所有输入事件同时发生
• 或门：任一输入事件发生

最小割集

导致顶事件的最小底事件组合。

---

8. Python实现与案例分析

8.1 案例1：蒙特卡洛不确定性传播

"""
案例1：蒙特卡洛方法进行不确定性传播
悬臂梁挠度分析
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def cantilever_deflection(P, L, E, I, b, h):
    """
    计算悬臂梁端部挠度
    
    参数:
        P: 端部载荷 (N)
        L: 梁长度 (m)
        E: 弹性模量 (Pa)
        I: 截面惯性矩 (m^4)
        b: 截面宽度 (m)
        h: 截面高度 (m)
    
    返回:
        delta: 端部挠度 (m)
    """
    I = b * h**3 / 12
    delta = P * L**3 / (3 * E * I)
    return delta

def case1_mc_uncertainty_propagation():
    """案例1：蒙特卡洛不确定性传播"""
    print("="*60)
    print("案例1：蒙特卡洛不确定性传播 - 悬臂梁挠度分析")
    print("="*60)
    
    # 定义输入变量的不确定性
    # 载荷 P: 正态分布，均值1000N，标准差100N
    P_dist = stats.norm(loc=1000, scale=100)
    
    # 长度 L: 正态分布，均值2m，标准差0.02m
    L_dist = stats.norm(loc=2.0, scale=0.02)
    
    # 弹性模量 E: 对数正态分布，均值70GPa，变异系数5%
    E_dist = stats.lognorm(s=0.05, scale=np.exp(np.log(70e9)))
    
    # 截面宽度 b: 均匀分布，[0.09, 0.11]m
    b_dist = stats.uniform(loc=0.09, scale=0.02)
    
    # 截面高度 h: 均匀分布，[0.19, 0.21]m
    h_dist = stats.uniform(loc=0.19, scale=0.02)
    
    # 蒙特卡洛模拟
    n_samples = 10000
    print(f"\n运行蒙特卡洛模拟，样本数: {n_samples}")
    
    # 生成样本
    P_samples = P_dist.rvs(n_samples)
    L_samples = L_dist.rvs(n_samples)
    E_samples = E_dist.rvs(n_samples)
    b_samples = b_dist.rvs(n_samples)
    h_samples = h_dist.rvs(n_samples)
    
    # 计算输出
    delta_samples = np.array([
        cantilever_deflection(P, L, E, None, b, h)
        for P, L, E, b, h in zip(P_samples, L_samples, E_samples, b_samples, h_samples)
    ])
    
    # 统计结果
    delta_mean = np.mean(delta_samples)
    delta_std = np.std(delta_samples)
    delta_cov = delta_std / delta_mean
    
    print(f"\n挠度统计结果:")
    print(f"  均值: {delta_mean*1000:.4f} mm")
    print(f"  标准差: {delta_std*1000:.4f} mm")
    print(f"  变异系数: {delta_cov*100:.2f}%")
    print(f"  95%置信区间: [{np.percentile(delta_samples, 2.5)*1000:.4f}, {np.percentile(delta_samples, 97.5)*1000:.4f}] mm")
    
    # 可靠性分析（假设允许挠度为5mm）
    delta_allow = 0.005  # 5mm
    failure_samples = delta_samples > delta_allow
    pf_mc = np.sum(failure_samples) / n_samples
    beta_mc = -stats.norm.ppf(pf_mc)
    
    print(f"\n可靠性分析 (允许挠度 = {delta_allow*1000}mm):")
    print(f"  失效概率: {pf_mc:.6f}")
    print(f"  可靠度指标: {beta_mc:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 输入分布
    axes[0, 0].hist(P_samples, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
    x_P = np.linspace(P_dist.ppf(0.001), P_dist.ppf(0.999), 100)
    axes[0, 0].plot(x_P, P_dist.pdf(x_P), 'r-', linewidth=2)
    axes[0, 0].set_title('Load P Distribution')
    axes[0, 0].set_xlabel('P (N)')
    axes[0, 0].set_ylabel('PDF')
    axes[0, 0].grid(True, alpha=0.3)
    
    axes[0, 1].hist(L_samples, bins=50, density=True, alpha=0.7, color='green', edgecolor='black')
    axes[0, 1].set_title('Length L Distribution')
    axes[0, 1].set_xlabel('L (m)')
    axes[0, 1].set_ylabel('PDF')
    axes[0, 1].grid(True, alpha=0.3)
    
    axes[0, 2].hist(E_samples/1e9, bins=50, density=True, alpha=0.7, color='orange', edgecolor='black')
    axes[0, 2].set_title('Elastic Modulus E Distribution')
    axes[0, 2].set_xlabel('E (GPa)')
    axes[0, 2].set_ylabel('PDF')
    axes[0, 2].grid(True, alpha=0.3)
    
    # 输出分布
    axes[1, 0].hist(delta_samples*1000, bins=50, density=True, alpha=0.7, color='purple', edgecolor='black')
    axes[1, 0].axvline(delta_mean*1000, color='red', linestyle='--', linewidth=2, label=f'Mean: {delta_mean*1000:.2f}mm')
    axes[1, 0].axvline(delta_allow*1000, color='green', linestyle='--', linewidth=2, label=f'Limit: {delta_allow*1000}mm')
    axes[1, 0].set_title('Deflection Distribution')
    axes[1, 0].set_xlabel('Deflection (mm)')
    axes[1, 0].set_ylabel('PDF')
    axes[1, 0].legend()
    axes[1, 0].grid(True, alpha=0.3)
    
    # CDF
    sorted_delta = np.sort(delta_samples)
    cdf = np.arange(1, len(sorted_delta)+1) / len(sorted_delta)
    axes[1, 1].plot(sorted_delta*1000, cdf, 'b-', linewidth=2)
    axes[1, 1].axvline(delta_allow*1000, color='green', linestyle='--', linewidth=2)
    axes[1, 1].axhline(1-pf_mc, color='red', linestyle='--', linewidth=2, label=f'Reliability: {1-pf_mc:.4f}')
    axes[1, 1].set_title('Cumulative Distribution Function')
    axes[1, 1].set_xlabel('Deflection (mm)')
    axes[1, 1].set_ylabel('CDF')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    # 收敛分析
    sample_sizes = np.logspace(2, 4, 20).astype(int)
    means = []
    stds = []
    for n in sample_sizes:
        means.append(np.mean(delta_samples[:n]))
        stds.append(np.std(delta_samples[:n]))
    
    axes[1, 2].semilogx(sample_sizes, np.array(means)*1000, 'o-', label='Mean', linewidth=2)
    ax2 = axes[1, 2].twinx()
    ax2.semilogx(sample_sizes, np.array(stds)*1000, 's-', color='red', label='Std', linewidth=2)
    axes[1, 2].set_title('Monte Carlo Convergence')
    axes[1, 2].set_xlabel('Number of Samples')
    axes[1, 2].set_ylabel('Mean (mm)', color='blue')
    ax2.set_ylabel('Std (mm)', color='red')
    axes[1, 2].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题045\\output\\case1_mc_uncertainty.png', dpi=150)
    plt.close()
    
    print("\n案例1完成！")
    return delta_samples


# 运行案例1
if __name__ == "__main__":
    case1_mc_uncertainty_propagation()

8.2 案例2：可靠性分析与FORM

"""
案例2：一阶可靠性方法（FORM）
悬臂梁可靠性分析
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import minimize

def limit_state_function(x):
    """
    极限状态函数
    x = [P, L, E, b, h]
    g > 0: 安全
    g <= 0: 失效
    """
    P, L, E, b, h = x
    I = b * h**3 / 12
    delta = P * L**3 / (3 * E * I)
    delta_allow = 0.005  # 允许挠度5mm
    return delta_allow - delta

def transform_to_standard(x, dist_params):
    """变换到标准正态空间"""
    u = np.zeros_like(x)
    for i, (xi, params) in enumerate(zip(x, dist_params)):
        dist_type = params['type']
        if dist_type == 'normal':
            mean, std = params['mean'], params['std']
            u[i] = (xi - mean) / std
        elif dist_type == 'lognormal':
            mean, std = params['mean'], params['std']
            sigma_ln = np.sqrt(np.log(1 + (std/mean)**2))
            mu_ln = np.log(mean) - 0.5 * sigma_ln**2
            u[i] = (np.log(xi) - mu_ln) / sigma_ln
        elif dist_type == 'uniform':
            a, b = params['a'], params['b']
            # 使用Rosenblatt变换
            u[i] = stats.norm.ppf((xi - a) / (b - a))
    return u

def transform_from_standard(u, dist_params):
    """从标准正态空间变换回来"""
    x = np.zeros_like(u)
    for i, (ui, params) in enumerate(zip(u, dist_params)):
        dist_type = params['type']
        if dist_type == 'normal':
            mean, std = params['mean'], params['std']
            x[i] = mean + ui * std
        elif dist_type == 'lognormal':
            mean, std = params['mean'], params['std']
            sigma_ln = np.sqrt(np.log(1 + (std/mean)**2))
            mu_ln = np.log(mean) - 0.5 * sigma_ln**2
            x[i] = np.exp(mu_ln + ui * sigma_ln)
        elif dist_type == 'uniform':
            a, b = params['a'], params['params']['b']
            x[i] = a + stats.norm.cdf(ui) * (b - a)
    return x

def case2_form_reliability():
    """案例2：FORM可靠性分析"""
    print("="*60)
    print("案例2：一阶可靠性方法（FORM）")
    print("="*60)
    
    # 定义输入分布参数
    dist_params = [
        {'type': 'normal', 'mean': 1000, 'std': 100},      # P
        {'type': 'normal', 'mean': 2.0, 'std': 0.02},      # L
        {'type': 'lognormal', 'mean': 70e9, 'std': 3.5e9}, # E
        {'type': 'uniform', 'a': 0.09, 'b': 0.11},         # b
        {'type': 'uniform', 'a': 0.19, 'b': 0.21}          # h
    ]
    
    # 均值点
    x_mean = np.array([1000, 2.0, 70e9, 0.10, 0.20])
    
    # 在标准正态空间定义极限状态函数
    def g_u(u):
        x = transform_from_standard(u, dist_params)
        return limit_state_function(x)
    
    # 使用优化求解MPP
    print("\n求解最可能失效点（MPP）...")
    
    def objective(u):
        return np.linalg.norm(u)
    
    def constraint(u):
        return g_u(u)
    
    # 初始猜测
    u0 = np.zeros(5)
    
    # 约束优化
    result = minimize(objective, u0, method='SLSQP', 
                     constraints={'type': 'eq', 'fun': constraint})
    
    u_star = result.x
    beta = np.linalg.norm(u_star)
    x_star = transform_from_standard(u_star, dist_params)
    
    # FORM失效概率
    pf_form = stats.norm.cdf(-beta)
    
    print(f"\nFORM结果:")
    print(f"  可靠度指标 β = {beta:.4f}")
    print(f"  失效概率 Pf = {pf_form:.6f}")
    print(f"\nMPP点 (标准空间): {u_star}")
    print(f"MPP点 (原始空间):")
    print(f"  P = {x_star[0]:.2f} N")
    print(f"  L = {x_star[1]:.4f} m")
    print(f"  E = {x_star[2]/1e9:.2f} GPa")
    print(f"  b = {x_star[3]:.4f} m")
    print(f"  h = {x_star[4]:.4f} m")
    
    # 计算方向余弦（敏感性因子）
    alpha = u_star / beta
    print(f"\n方向余弦（敏感性因子）:")
    var_names = ['P', 'L', 'E', 'b', 'h']
    for name, a in zip(var_names, alpha):
        print(f"  {name}: {a:.4f} (贡献: {a**2*100:.2f}%)")
    
    # 蒙特卡洛验证
    print("\n蒙特卡洛验证...")
    n_mc = 100000
    P_samples = stats.norm.rvs(1000, 100, n_mc)
    L_samples = stats.norm.rvs(2.0, 0.02, n_mc)
    E_samples = stats.lognorm.rvs(0.05, scale=70e9, size=n_mc)
    b_samples = stats.uniform.rvs(0.09, 0.02, n_mc)
    h_samples = stats.uniform.rvs(0.19, 0.02, n_mc)
    
    g_samples = np.array([
        limit_state_function([P, L, E, b, h])
        for P, L, E, b, h in zip(P_samples, L_samples, E_samples, b_samples, h_samples)
    ])
    
    pf_mc = np.sum(g_samples <= 0) / n_mc
    
    print(f"  MC失效概率: {pf_mc:.6f}")
    print(f"  FORM失效概率: {pf_form:.6f}")
    print(f"  相对误差: {abs(pf_form - pf_mc)/pf_mc*100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 标准正态空间中的极限状态面（投影到2D）
    u1_range = np.linspace(-4, 4, 100)
    u2_range = np.linspace(-4, 4, 100)
    U1, U2 = np.meshgrid(u1_range, u2_range)
    
    # 固定其他变量为0，只变化P和L
    G_grid = np.zeros_like(U1)
    for i in range(len(u1_range)):
        for j in range(len(u2_range)):
            u_temp = np.array([U1[j,i], U2[j,i], 0, 0, 0])
            G_grid[j,i] = g_u(u_temp)
    
    contour = axes[0, 0].contour(U1, U2, G_grid, levels=[0], colors='red', linewidths=2)
    axes[0, 0].clabel(contour, inline=True, fontsize=10)
    axes[0, 0].plot(u_star[0], u_star[1], 'ro', markersize=10, label='MPP')
    axes[0, 0].plot(0, 0, 'go', markersize=8, label='Origin')
    axes[0, 0].arrow(0, 0, u_star[0], u_star[1], head_width=0.1, head_length=0.1, 
                    fc='blue', ec='blue', alpha=0.5)
    circle = plt.Circle((0, 0), beta, fill=False, color='green', linestyle='--', linewidth=2)
    axes[0, 0].add_patch(circle)
    axes[0, 0].set_xlabel('u_P (standardized P)')
    axes[0, 0].set_ylabel('u_L (standardized L)')
    axes[0, 0].set_title('Limit State Surface in Standard Normal Space')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].axis('equal')
    
    # 敏感性因子
    axes[0, 1].barh(var_names, alpha**2, color='steelblue', alpha=0.7)
    axes[0, 1].set_xlabel('Sensitivity Factor α²')
    axes[0, 1].set_title('Importance Factors')
    axes[0, 1].grid(True, alpha=0.3, axis='x')
    
    # 失效概率对比
    methods = ['FORM', 'Monte Carlo']
    pf_values = [pf_form, pf_mc]
    axes[1, 0].bar(methods, pf_values, color=['orange', 'green'], alpha=0.7)
    axes[1, 0].set_ylabel('Failure Probability')
    axes[1, 0].set_title('Failure Probability Comparison')
    for i, v in enumerate(pf_values):
        axes[1, 0].text(i, v + pf_mc*0.1, f'{v:.2e}', ha='center')
    axes[1, 0].grid(True, alpha=0.3, axis='y')
    
    # 极限状态函数值分布
    axes[1, 1].hist(g_samples, bins=50, density=True, alpha=0.7, color='purple', edgecolor='black')
    axes[1, 1].axvline(0, color='red', linestyle='--', linewidth=2, label='Limit State (g=0)')
    axes[1, 1].axvline(np.mean(g_samples), color='blue', linestyle='--', linewidth=2, 
                      label=f'Mean: {np.mean(g_samples):.4f}')
    axes[1, 1].set_xlabel('g(X)')
    axes[1, 1].set_ylabel('PDF')
    axes[1, 1].set_title('Limit State Function Distribution')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题045\\output\\case2_form_reliability.png', dpi=150)
    plt.close()
    
    print("\n案例2完成！")
    return beta, pf_form


# 运行案例2
if __name__ == "__main__":
    case2_form_reliability()

8.3 案例3：全局敏感性分析

"""
案例3：全局敏感性分析 - Sobol指数
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from matplotlib.patches import Rectangle

def ishigami_function(x):
    """
    Ishigami函数 - 常用于敏感性分析测试
    具有强非线性和交互作用
    """
    x1, x2, x3 = x
    return np.sin(x1) + 7*np.sin(x2)**2 + 0.1*x3**4*np.sin(x1)

def sobol_indices_mc(n_samples, model, n_inputs=3):
    """
    基于蒙特卡洛的Sobol指数估计（Saltelli方案）
    """
    # 生成两个独立样本矩阵
    A = np.random.rand(n_samples, n_inputs) * 2 * np.pi - np.pi
    B = np.random.rand(n_samples, n_inputs) * 2 * np.pi - np.pi
    
    # 计算模型输出
    Y_A = np.array([model(a) for a in A])
    Y_B = np.array([model(b) for b in B])
    
    # 总方差
    V_Y = np.var(Y_A)
    
    # 一阶和总效应指数
    S1 = np.zeros(n_inputs)
    ST = np.zeros(n_inputs)
    
    for i in range(n_inputs):
        # 创建混合矩阵 A_B
        A_B = A.copy()
        A_B[:, i] = B[:, i]
        Y_AB = np.array([model(a_b) for a_b in A_B])
        
        # 一阶指数
        S1[i] = np.mean(Y_B * (Y_AB - Y_A)) / V_Y
        
        # 总效应指数
        ST[i] = 0.5 * np.mean((Y_A - Y_AB)**2) / V_Y
    
    return S1, ST, V_Y

def case3_sobol_sensitivity():
    """案例3：Sobol敏感性分析"""
    print("="*60)
    print("案例3：全局敏感性分析 - Sobol指数")
    print("="*60)
    
    # 使用Ishigami函数进行演示
    print("\n使用Ishigami函数进行敏感性分析")
    print("Y = sin(X1) + 7*sin(X2)^2 + 0.1*X3^4*sin(X1)")
    print("其中 X1, X2, X3 ~ U(-π, π)")
    
    # 解析解（用于验证）
    print("\n解析解:")
    print("  S1 = 0.3139 (X1)")
    print("  S2 = 0.4424 (X2)")
    print("  S3 = 0.0000 (X3)")
    print("  S12 = 0.2437 (X1-X2交互)")
    print("  ST1 = 0.5576")
    print("  ST2 = 0.4424")
    print("  ST3 = 0.2437")
    
    # 蒙特卡洛估计
    n_samples = 10000
    print(f"\n蒙特卡洛估计（样本数: {n_samples}）...")
    
    S1, ST, V_Y = sobol_indices_mc(n_samples, ishigami_function, 3)
    
    print(f"\n一阶Sobol指数:")
    print(f"  S1 = {S1[0]:.4f} (X1)")
    print(f"  S2 = {S1[1]:.4f} (X2)")
    print(f"  S3 = {S1[2]:.4f} (X3)")
    print(f"  总和: {np.sum(S1):.4f}")
    
    print(f"\n总效应指数:")
    print(f"  ST1 = {ST[0]:.4f} (X1)")
    print(f"  ST2 = {ST[1]:.4f} (X2)")
    print(f"  ST3 = {ST[2]:.4f} (X3)")
    
    # 交互作用
    S_interaction = ST - S1
    print(f"\n交互作用:")
    print(f"  X1交互: {S_interaction[0]:.4f}")
    print(f"  X2交互: {S_interaction[1]:.4f}")
    print(f"  X3交互: {S_interaction[2]:.4f}")
    
    # Morris方法（筛选）
    print("\nMorris方法筛选...")
    
    def morris_screening(n_trajectories, n_inputs, delta=0.1):
        """Morris筛选方法"""
        EE = np.zeros((n_trajectories, n_inputs))
        
        for t in range(n_trajectories):
            # 随机起点
            x0 = np.random.rand(n_inputs)
            
            for i in range(n_inputs):
                # 创建轨迹
                x = x0.copy()
                y_base = ishigami_function(x * 2 * np.pi - np.pi)
                
                x[i] = (x[i] + delta) % 1
                y_perturbed = ishigami_function(x * 2 * np.pi - np.pi)
                
                EE[t, i] = (y_perturbed - y_base) / delta
        
        # 统计量
        mu = np.mean(EE, axis=0)
        sigma = np.std(EE, axis=0)
        mu_star = np.mean(np.abs(EE), axis=0)
        
        return mu, sigma, mu_star
    
    mu, sigma, mu_star = morris_screening(100, 3)
    
    print(f"\nMorris方法结果:")
    print(f"  μ* (平均绝对效应):")
    print(f"    X1: {mu_star[0]:.4f}")
    print(f"    X2: {mu_star[1]:.4f}")
    print(f"    X3: {mu_star[2]:.4f}")
    print(f"  σ (标准差，非线性指标):")
    print(f"    X1: {sigma[0]:.4f}")
    print(f"    X2: {sigma[1]:.4f}")
    print(f"    X3: {sigma[2]:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    var_names = ['X1', 'X2', 'X3']
    
    # Sobol指数对比
    x_pos = np.arange(len(var_names))
    width = 0.35
    axes[0, 0].bar(x_pos - width/2, S1, width, label='First-order', alpha=0.7, color='steelblue')
    axes[0, 0].bar(x_pos + width/2, ST, width, label='Total', alpha=0.7, color='coral')
    axes[0, 0].set_xlabel('Input Variable')
    axes[0, 0].set_ylabel('Sobol Index')
    axes[0, 0].set_title('Sobol Sensitivity Indices')
    axes[0, 0].set_xticks(x_pos)
    axes[0, 0].set_xticklabels(var_names)
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3, axis='y')
    
    # 一阶vs总效应
    axes[0, 1].scatter(S1, ST, s=100, alpha=0.7, color='green')
    for i, name in enumerate(var_names):
        axes[0, 1].annotate(name, (S1[i], ST[i]), xytext=(5, 5), 
                           textcoords='offset points', fontsize=12)
    axes[0, 1].plot([0, 1], [0, 1], 'k--', alpha=0.5)
    axes[0, 1].set_xlabel('First-order Index S1')
    axes[0, 1].set_ylabel('Total Index ST')
    axes[0, 1].set_title('First-order vs Total Effect')
    axes[0, 1].grid(True, alpha=0.3)
    
    # Morris筛选图
    axes[0, 2].scatter(mu_star, sigma, s=100, alpha=0.7, color='purple')
    for i, name in enumerate(var_names):
        axes[0, 2].annotate(name, (mu_star[i], sigma[i]), xytext=(5, 5), 
                           textcoords='offset points', fontsize=12)
    axes[0, 2].set_xlabel('μ* (Mean Absolute Effect)')
    axes[0, 2].set_ylabel('σ (Standard Deviation)')
    axes[0, 2].set_title('Morris Screening Plot')
    axes[0, 2].grid(True, alpha=0.3)
    
    # 输出分布
    n_samples_dist = 10000
    samples = np.random.rand(n_samples_dist, 3) * 2 * np.pi - np.pi
    Y_samples = np.array([ishigami_function(s) for s in samples])
    
    axes[1, 0].hist(Y_samples, bins=50, density=True, alpha=0.7, color='teal', edgecolor='black')
    axes[1, 0].set_xlabel('Output Y')
    axes[1, 0].set_ylabel('PDF')
    axes[1, 0].set_title('Output Distribution')
    axes[1, 0].grid(True, alpha=0.3)
    
    # 散点图矩阵
    sample_size = 500
    idx = np.random.choice(n_samples_dist, sample_size, replace=False)
    
    # X1 vs Y
    axes[1, 1].scatter(samples[idx, 0], Y_samples[idx], alpha=0.5, s=10, color='blue')
    axes[1, 1].set_xlabel('X1')
    axes[1, 1].set_ylabel('Y')
    axes[1, 1].set_title('Scatter: X1 vs Y')
    axes[1, 1].grid(True, alpha=0.3)
    
    # X2 vs Y
    axes[1, 2].scatter(samples[idx, 1], Y_samples[idx], alpha=0.5, s=10, color='red')
    axes[1, 2].set_xlabel('X2')
    axes[1, 2].set_ylabel('Y')
    axes[1, 2].set_title('Scatter: X2 vs Y')
    axes[1, 2].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题045\\output\\case3_sobol_sensitivity.png', dpi=150)
    plt.close()
    
    print("\n案例3完成！")
    return S1, ST


# 运行案例3
if __name__ == "__main__":
    case3_sobol_sensitivity()

8.4 案例4：时变可靠性分析

"""
案例4：时变可靠性分析
随机振动系统的首次穿越问题
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from matplotlib.animation import FuncAnimation, PillowWriter

def generate_random_process(t, dt, omega_0, zeta, S0):
    """
    生成高斯随机过程（白噪声激励的线性系统响应）
    
    参数:
        t: 时间数组
        dt: 时间步长
        omega_0: 固有频率
        zeta: 阻尼比
        S0: 白噪声谱密度
    
    返回:
        x: 位移响应
        v: 速度响应
    """
    n_steps = len(t)
    x = np.zeros(n_steps)
    v = np.zeros(n_steps)
    
    # 白噪声激励
    w = np.random.randn(n_steps) * np.sqrt(2 * np.pi * S0 / dt)
    
    # 数值积分（Newmark-beta方法）
    for i in range(n_steps - 1):
        # 加速度
        a = (w[i] - 2 * zeta * omega_0 * v[i] - omega_0**2 * x[i])
        
        # 速度和位移更新
        v[i+1] = v[i] + a * dt
        x[i+1] = x[i] + v[i+1] * dt
    
    return x, v

def case4_time_dependent_reliability():
    """案例4：时变可靠性分析"""
    print("="*60)
    print("案例4：时变可靠性分析 - 随机振动系统")
    print("="*60)
    
    # 系统参数
    omega_0 = 2 * np.pi * 5  # 固有频率 5Hz
    zeta = 0.05              # 阻尼比
    S0 = 0.1                 # 白噪声谱密度
    threshold = 0.5          # 位移阈值
    
    # 时间参数
    T = 10.0                 # 总时间
    dt = 0.01                # 时间步长
    t = np.arange(0, T, dt)
    
    # 蒙特卡洛模拟
    n_samples = 1000
    print(f"\n运行蒙特卡洛模拟，样本数: {n_samples}")
    print(f"时间范围: 0-{T}s，步长: {dt}s")
    
    trajectories = []
    first_passage_times = []
    n_failures = 0
    
    for i in range(n_samples):
        x, v = generate_random_process(t, dt, omega_0, zeta, S0)
        trajectories.append(x)
        
        # 检测首次穿越
        exceed_idx = np.where(np.abs(x) > threshold)[0]
        if len(exceed_idx) > 0:
            first_passage_times.append(t[exceed_idx[0]])
            n_failures += 1
        else:
            first_passage_times.append(None)
    
    # 时变失效概率
    failure_times = [fp for fp in first_passage_times if fp is not None]
    
    print(f"\n统计结果:")
    print(f"  失效样本数: {n_failures}/{n_samples} ({n_failures/n_samples*100:.2f}%)")
    if failure_times:
        print(f"  平均首次穿越时间: {np.mean(failure_times):.4f}s")
        print(f"  首次穿越时间标准差: {np.std(failure_times):.4f}s")
    
    # 计算时变失效概率
    pf_time = np.zeros(len(t))
    for i, ti in enumerate(t):
        count = sum(1 for fp in first_passage_times if fp is not None and fp <= ti)
        pf_time[i] = count / n_samples
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 样本轨迹
    n_plot = min(50, n_samples)
    for i in range(n_plot):
        axes[0, 0].plot(t, trajectories[i], alpha=0.3, linewidth=0.5)
    axes[0, 0].axhline(threshold, color='red', linestyle='--', linewidth=2, label=f'Threshold = {threshold}')
    axes[0, 0].axhline(-threshold, color='red', linestyle='--', linewidth=2)
    axes[0, 0].set_xlabel('Time (s)')
    axes[0, 0].set_ylabel('Displacement')
    axes[0, 0].set_title(f'Sample Trajectories (n={n_plot})')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)
    
    # 时变失效概率
    axes[0, 1].plot(t, pf_time, 'b-', linewidth=2)
    axes[0, 1].set_xlabel('Time (s)')
    axes[0, 1].set_ylabel('Failure Probability')
    axes[0, 1].set_title('Time-Dependent Failure Probability')
    axes[0, 1].grid(True, alpha=0.3)
    
    # 首次穿越时间分布
    if failure_times:
        axes[1, 0].hist(failure_times, bins=30, density=True, alpha=0.7, color='green', edgecolor='black')
        axes[1, 0].set_xlabel('First Passage Time (s)')
        axes[1, 0].set_ylabel('PDF')
        axes[1, 0].set_title('First Passage Time Distribution')
        axes[1, 0].grid(True, alpha=0.3)
    
    # 包络统计
    max_values = [np.max(np.abs(traj)) for traj in trajectories]
    axes[1, 1].hist(max_values, bins=30, density=True, alpha=0.7, color='purple', edgecolor='black')
    axes[1, 1].axvline(threshold, color='red', linestyle='--', linewidth=2, label=f'Threshold = {threshold}')
    axes[1, 1].set_xlabel('Maximum Absolute Displacement')
    axes[1, 1].set_ylabel('PDF')
    axes[1, 1].set_title('Maximum Response Distribution')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题045\\output\\case4_time_dependent_reliability.png', dpi=150)
    plt.close()
    
    print("\n案例4完成！")
    return t, pf_time, trajectories


def create_uncertainty_animation():
    """创建不确定性传播动画"""
    print("\n创建不确定性传播动画...")
    
    # 生成样本
    n_frames = 50
    n_samples_per_frame = 20
    
    # 输入分布
    P_dist = stats.norm(loc=1000, scale=100)
    L_dist = stats.norm(loc=2.0, scale=0.02)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    def update(frame):
        axes[0].clear()
        axes[1].clear()
        
        # 累积样本
        n_total = (frame + 1) * n_samples_per_frame
        P_samples = P_dist.rvs(n_total)
        L_samples = L_dist.rvs(n_total)
        
        # 计算挠度
        delta_samples = []
        for P, L in zip(P_samples, L_samples):
            I = 0.10 * 0.20**3 / 12
            delta = P * L**3 / (3 * 70e9 * I)
            delta_samples.append(delta * 1000)  # 转换为mm
        
        # 左图：输入空间
        axes[0].scatter(P_samples, L_samples, alpha=0.5, s=20, c='blue')
        axes[0].set_xlabel('Load P (N)')
        axes[0].set_ylabel('Length L (m)')
        axes[0].set_title(f'Input Space (n={n_total})')
        axes[0].grid(True, alpha=0.3)
        axes[0].set_xlim(700, 1300)
        axes[0].set_ylim(1.9, 2.1)
        
        # 右图：输出分布
        axes[1].hist(delta_samples, bins=30, density=True, alpha=0.7, 
                    color='purple', edgecolor='black')
        axes[1].axvline(np.mean(delta_samples), color='red', linestyle='--', 
                       linewidth=2, label=f'Mean: {np.mean(delta_samples):.2f}mm')
        axes[1].set_xlabel('Deflection (mm)')
        axes[1].set_ylabel('PDF')
        axes[1].set_title(f'Output Distribution (n={n_total})')
        axes[1].legend()
        axes[1].grid(True, alpha=0.3)
        
        plt.tight_layout()
        return axes
    
    anim = FuncAnimation(fig, update, frames=n_frames, interval=200, blit=False)
    writer = PillowWriter(fps=5)
    anim.save('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题045\\output\\uncertainty_propagation_animation.gif', 
              writer=writer)
    plt.close()
    
    print("动画创建完成！")


def main():
    """主函数"""
    print("="*60)
    print("主题045：不确定性量化与可靠性分析")
    print("="*60)
    
    # 运行所有案例
    case1_mc_uncertainty_propagation()
    case2_form_reliability()
    case3_sobol_sensitivity()
    case4_time_dependent_reliability()
    create_uncertainty_animation()
    
    print("\n" + "="*60)
    print("所有案例完成！")
    print("="*60)
    print("\n输出文件:")
    print("  - case1_mc_uncertainty.png")
    print("  - case2_form_reliability.png")
    print("  - case3_sobol_sensitivity.png")
    print("  - case4_time_dependent_reliability.png")
    print("  - uncertainty_propagation_animation.gif")


if __name__ == "__main__":
    main()

---

9. 总结与展望

9.1 方法对比与选择

方法	优点	缺点	适用场景
蒙特卡洛	通用、精确、易并行	计算量大	复杂模型、高精度要求
FORM/SORM	计算效率高	非线性误差	中等非线性、快速评估
PCE	解析矩计算	高维困难	光滑响应、多次评估
代理模型	加速计算	近似误差	昂贵仿真、优化设计

9.2 发展趋势

深度学习方法

• 神经网络代理模型
• 深度高斯过程
• 物理信息神经网络

高维问题处理

• 稀疏多项式混沌
• 主动子空间方法
• 流形学习技术

动态可靠性

• 时变可靠性分析
• 退化过程建模
• 剩余寿命预测

多尺度不确定性

• 微观-宏观关联
• 多尺度随机建模
• 异构不确定性融合

9.3 工程应用建议

1. 问题识别：明确不确定性来源和类型
2. 方法选择：根据问题特点选择合适方法
3. 验证确认：通过多种方法交叉验证
4. 灵敏度分析：识别关键不确定性因素
5. 决策支持：将UQ结果融入设计决策

---

6 深入理论分析

6.1 数学模型推导

流体力学问题的数学描述基于守恒定律和本构关系。控制方程的一般形式可以表示为：

$$\frac{\partial (\rho \varphi )}{\partial t}+\nabla \cdot (\rho \mathbf{u}\varphi )=\nabla \cdot (\Gamma \nabla \varphi )+S_{\varphi }$$

其中，$\varphi$表示通用变量，$\rho$是密度，$\mathbf{u}$是速度矢量，$\Gamma$是扩散系数，$S_{\varphi }$是源项。

对于稳态问题，控制方程简化为：

$$\nabla \cdot (\rho \mathbf{u}\varphi )=\nabla \cdot (\Gamma \nabla \varphi )+S_{\varphi }$$

边界条件的数学表达：

Dirichlet边界条件：
$$\varphi ={\varphi }_{specified}\text{在}{\Gamma }_D\text{上}$$

Neumann边界条件：
$$-\Gamma \frac{\partial \varphi }{\partial n}=q_{specified}\text{在}{\Gamma }_N\text{上}$$

Robin边界条件：
$$-\Gamma \frac{\partial \varphi }{\partial n}=h(\varphi -{\varphi }_{\infty })\text{在}{\Gamma }_R\text{上}$$

6.2 数值离散方法

有限差分法的离散格式：

时间离散：

• 显式格式：${\varphi }^{n+1}={\varphi }^n+\Delta t\cdot f({\varphi }^n)$
• 隐式格式：${\varphi }^{n+1}={\varphi }^n+\Delta t\cdot f({\varphi }^{n+1})$
• Crank-Nicolson格式：${\varphi }^{n+1}={\varphi }^n+\frac{\Delta t}{2}[f({\varphi }^n)+f({\varphi }^{n+1})]$

空间离散：

• 中心差分：$\frac{\partial \varphi }{\partial x}\approx \frac{{\varphi }_{i+1}-{\varphi }_{i-1}}{2\Delta x}$
• 迎风格式：$\frac{\partial \varphi }{\partial x}\approx \frac{{\varphi }_i-{\varphi }_{i-1}}{\Delta x}$（$u>0$时）

6.3 稳定性与收敛性分析

Von Neumann稳定性分析：

对于一维热传导方程的显式格式，稳定性条件为：

$$Fo=\frac{\alpha \Delta t}{(\Delta x{)}^2}\le \frac{1}{2}$$

其中$Fo$是傅里叶数，$\alpha$是热扩散系数。

收敛性条件：

根据Lax等价定理，对于适定的线性初值问题，一致性和稳定性是收敛性的充分必要条件。

误差分析：

截断误差：$\tau =O(\Delta t)+O((\Delta x{)}^2)$

累积误差：$ϵ=O(\Delta t)+O((\Delta x{)}^2)$