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

一、引言

在工程仿真领域，不确定性无处不在。材料属性的分散性、几何尺寸的制造公差、边界条件的波动、以及模型本身的简化假设，都会引入各种不确定性。传统的确定性分析方法无法捕捉这些不确定性对系统响应的影响，可能导致设计过于保守或存在安全隐患。不确定性量化（Uncertainty Quantification, UQ）与可靠性分析为处理这些问题提供了系统化的方法论框架。

不确定性量化旨在识别、表征和传播不确定性，从而评估系统响应的统计特性。可靠性分析则在此基础上，计算系统在特定条件下完成预定功能的概率。这两个领域相辅相成，共同构成了现代工程风险评估和鲁棒性设计的理论基础。

本章将系统介绍不确定性量化与可靠性分析的基本理论、数值方法和工程应用。内容包括不确定性的分类与建模、不确定性传播方法、可靠性分析算法、灵敏度分析技术，以及基于可靠性的优化设计方法。通过丰富的Python仿真案例，帮助读者深入理解这些方法的原理和应用。

二、不确定性的分类与建模

2.1 不确定性的来源

工程系统中的不确定性主要来源于以下几个方面：

1. 固有不确定性（Aleatory Uncertainty）

固有不确定性也称为偶然不确定性或随机不确定性，源于物理系统的内在随机性。这种不确定性是不可消除的，只能通过统计方法进行描述。典型例子包括：

• 材料力学性能的批次间差异
• 环境载荷的随机波动（风载、地震、波浪等）
• 制造过程中的尺寸分散
• 微观结构的随机分布

固有不确定性通常用概率模型描述，如随机变量、随机过程或随机场。

2. 认知不确定性（Epistemic Uncertainty）

认知不确定性也称为系统不确定性，源于知识的不完备或信息的缺乏。这种不确定性可以通过收集更多数据、改进模型或增加认知来减少。典型例子包括：

• 模型参数的估计误差
• 模型结构的选择不确定性
• 小样本条件下的统计推断不确定性
• 专家判断的主观性

认知不确定性可以用区间分析、模糊理论、证据理论或贝叶斯方法处理。

3. 模型不确定性

模型不确定性源于数学模型对物理现实的近似。即使输入参数完全确定，模型预测仍可能与真实响应存在偏差。来源包括：

• 控制方程的简化（如线性化假设）
• 边界条件和初始条件的近似
• 数值离散化误差
• 未建模的物理现象

模型验证和确认（Verification and Validation, V&V）是量化模型不确定性的重要手段。

2.2 随机变量的建模

当不确定性可以用概率分布描述时，随机变量是最基本的建模工具。

2.2.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$ 为标准差。正态分布的特点是关于均值对称，68-95-99.7规则描述了数据落在不同标准差范围内的概率。

在材料性能建模中，当变异系数（标准差/均值）较小时（通常<20%），正态分布是合理的选择。例如，铝合金的屈服强度、混凝土的抗压强度等。

对数正态分布（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$$

对数正态分布常用于描述疲劳寿命、裂纹尺寸、材料弹性模量等。其特点是能够避免负值的产生，且能描述较大的分散性。

威布尔分布（Weibull Distribution）

威布尔分布是可靠性工程中最常用的分布，特别适用于描述强度和寿命数据：

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

其中，$\beta$ 为形状参数，$\eta$ 为尺度参数。威布尔分布的灵活性在于：

• 当 $\beta <1$ 时，描述早期失效（递减失效率）
• 当 $\beta =1$ 时，退化为指数分布（恒定失效率）
• 当 $\beta >1$ 时，描述磨损失效（递增收效率）
• 当 $\beta \approx 3.4$ 时，接近正态分布

威布尔分布广泛用于脆性材料强度、轴承寿命、电子元器件可靠性分析。

均匀分布（Uniform Distribution）

当仅知道参数的变化范围而缺乏更多信息时，均匀分布是最保守的选择：

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

均匀分布假设在区间 $[a,b]$ 内所有值等可能，最大熵原理支持这一选择。

2.2.2 分布参数估计

给定样本数据，需要估计分布的参数。常用方法包括：

矩估计法（Method of Moments）

用样本矩估计总体矩，建立方程求解参数。对于正态分布：

$$\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$$

矩估计法简单直观，但对于小样本可能不够精确。

最大似然估计（Maximum Likelihood Estimation, MLE）

寻找使似然函数最大化的参数值。似然函数定义为：

$$L(\theta )=\prod _{i=1}^n{f}_X(x_i;\theta )$$

通常对似然函数取对数后求极值。MLE具有良好的渐近性质，是参数估计的首选方法。

贝叶斯估计

结合先验信息和样本数据，通过贝叶斯定理得到后验分布：

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

贝叶斯方法特别适合小样本情况，能够量化参数估计的不确定性。

2.3 随机过程与随机场

当不确定性随时间或空间变化时，需要采用随机过程或随机场建模。

2.3.1 随机过程

随机过程 {X(t),t∈T}\{X(t), t \in T\}{X(t),t∈T} 是一族随机变量，描述不确定性随时间的演化。关键统计特性包括：

均值函数：${\mu }_X(t)=E[X(t)]$

自相关函数：$R_X(t_1,t_2)=E[X(t_1)X(t_2)]$

协方差函数：$C_X(t_1,t_2)=E[(X(t_1)-{\mu }_X(t_1))(X(t_2)-{\mu }_X(t_2))]$

平稳随机过程满足：

• 均值恒定：${\mu }_X(t)=\mu$
• 自相关函数仅依赖于时间差：$R_X(t_1,t_2)=R_X(t_1-t_2)$

高斯随机过程是最重要的一类随机过程，其任意有限维分布都是多元正态分布。高斯过程完全由均值函数和协方差函数确定。

2.3.2 随机场

随机场 {X(s),s∈D}\{X(\mathbf{s}), \mathbf{s} \in D\}{X(s),s∈D} 描述空间分布的不确定性，$D$ 为空间域。随机场在岩土工程（土性参数空间变异）、复合材料（纤维分布）、微制造（表面粗糙度）等领域有重要应用。

相关结构描述空间点之间的统计关联。常用模型包括：

指数型相关函数：

$$\rho ({\mathbf{s}}_1,{\mathbf{s}}_2)=\exp \left(-\frac{\Vert {\mathbf{s}}_1-{\mathbf{s}}_2\Vert }{l_c}\right)$$

高斯型相关函数：

$$\rho ({\mathbf{s}}_1,{\mathbf{s}}_2)=\exp \left(-\frac{\Vert {\mathbf{s}}_1-{\mathbf{s}}_2{\Vert }^2}{l_c^2}\right)$$

其中，$l_c$ 为相关长度，表征空间相关性的衰减速度。

随机场的离散化：实际计算中需要将连续随机场离散为有限个随机变量。常用方法包括：

• 中心点法：用单元中心值代表整个单元
• 局部平均法：用单元内的平均值
• 形函数法：用节点值插值
• Karhunen-Loève展开：最优的谱展开方法

2.4 模糊与区间不确定性

当概率信息不足时，可以采用非概率方法描述不确定性。

2.4.1 区间分析

区间不确定性用区间 $[\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.4.2 模糊理论

模糊理论用隶属函数 ${\mu }_A(x)\in [0,1]$ 描述元素 $x$ 属于模糊集合 $A$ 的程度。隶属函数为1表示完全属于，为0表示完全不属于，中间值表示部分属于。

α-截集将模糊数转化为区间：

Aα={x∣μA(x)≥α}A_\alpha = \{x | \mu_A(x) \geq \alpha\}Aα​={x∣μA​(x)≥α}

通过在不同α水平上进行区间分析，可以得到模糊输出。

模糊理论适用于专家知识表达、语言变量描述等场景，但隶属函数的确定具有一定主观性。

三、不确定性传播方法

不确定性传播研究输入不确定性如何通过系统模型传递到输出响应。根据计算成本和精度要求，可选择不同方法。

3.1 蒙特卡洛模拟

蒙特卡洛模拟（Monte Carlo Simulation, MCS）是最通用的不确定性传播方法，通过大量随机抽样估计输出统计特性。

3.1.1 基本算法

1. 根据输入变量的概率分布，生成 $N$ 组随机样本 {x(i)}i=1N\{\mathbf{x}^{(i)}\}_{i=1}^N{x(i)}i=1N​
2. 对每个样本计算系统响应：$y^{(i)}=g({\mathbf{x}}^{(i)})$
3. 基于输出样本 {y(i)}\{y^{(i)}\}{y(i)} 估计统计量：
   • 均值：${\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$
   • 概率分布：经验分布函数或核密度估计

3.1.2 收敛性与样本量

蒙特卡洛估计的标准误差为：

$$ϵ=\frac{{\sigma }_Y}{\sqrt{N}}$$

要达到精度 $ϵ$，所需样本量为 $N={\sigma }_Y^2/{ϵ}^2$。这意味着：

• 误差以 $1/\sqrt{N}$ 的速度收敛，相对较慢
• 要达到1%的精度，通常需要 $10^4$ 次以上的计算
• 收敛速度与问题维度无关，这是蒙特卡洛方法的优势

3.1.3 方差缩减技术

为提高计算效率，发展了多种方差缩减技术：

重要性抽样（Importance Sampling）：从偏置分布中抽样，增加对输出有重要影响的区域的样本密度。

分层抽样（Stratified Sampling）：将输入空间划分为若干层，在每层内均匀抽样，确保覆盖性。

拉丁超立方抽样（Latin Hypercube Sampling, LHS）：将每个变量的分布划分为 $N$ 个等概率区间，在每个区间随机抽取一个值，然后随机组合。LHS能显著改善小样本下的估计精度。

拟蒙特卡洛（Quasi-Monte Carlo）：使用低差异序列（如Sobol序列、Halton序列）代替伪随机数，获得更均匀的样本分布。

3.2 泰勒展开方法

泰勒展开方法通过局部线性或二次近似传播不确定性，计算成本低但仅适用于低不确定性水平。

3.2.1 一阶泰勒展开（FOSM）

在均值点 ${\mathbf{\mu }}_X$ 处进行一阶泰勒展开：

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

输出均值和方差的近似：

$${\mu }_Y\approx g({\mathbf{\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}\text{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$$

3.2.2 二阶泰勒展开

考虑二阶项可以提高精度，特别是对于非线性问题：

$$Y\approx g({\mathbf{\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})$$

二阶展开能捕捉输出的偏度，但计算Hessian矩阵的成本较高。

3.3 多项式混沌展开

多项式混沌展开（Polynomial Chaos Expansion, PCE）是一种谱方法，用正交多项式基函数展开随机响应。

3.3.1 基本理论

根据Wiener-Askey方案，不同分布类型的随机变量对应不同的正交多项式：

• 高斯分布 → Hermite多项式
• 伽马分布 → Laguerre多项式
• 贝塔分布 → Jacobi多项式
• 均匀分布 → Legendre多项式

随机响应展开为：

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

其中，$\mathbf{\xi }$ 为标准随机变量，${\Psi }_{\mathbf{\alpha }}$ 为多维正交多项式，$y_{\mathbf{\alpha }}$ 为展开系数。

3.3.2 系数计算

展开系数可通过投影法或回归法计算：

投影法（Galerkin投影）：

$$y_{\mathbf{\alpha }}=\frac{E[Y{\Psi }_{\mathbf{\alpha }}]}{E[{\Psi }_{\mathbf{\alpha }}^2]}=\frac{1}{E[{\Psi }_{\mathbf{\alpha }}^2]}\int Y(\mathbf{\xi }){\Psi }_{\mathbf{\alpha }}(\mathbf{\xi })f_{\mathbf{\xi }}(\mathbf{\xi })d\mathbf{\xi }$$

积分可用高斯求积或蒙特卡洛方法计算。

回归法：在选定的配置点上计算响应，通过最小二乘拟合确定系数：

$$\underset{\mathbf{y}}{\min }\sum _{i=1}^N{\left(y^{(i)}-\sum _{\mathbf{\alpha }}{y}_{\mathbf{\alpha }}{\Psi }_{\mathbf{\alpha }}({\mathbf{\xi }}^{(i)})\right)}^2$$

3.3.3 后处理

获得PCE后，可以解析计算统计量：

$${\mu }_Y=y_0$$

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

Sobol灵敏度指标也可从PCE系数解析计算，这是PCE的重要优势。

3.4 随机有限元方法

随机有限元方法（Stochastic Finite Element Method, SFEM）将有限元分析与不确定性传播相结合。

3.4.1 摄动法

将随机刚度矩阵和响应在均值处展开：

$$\mathbf{K}(\mathbf{\xi })={\mathbf{K}}_0+\sum _{i=1}^n{\mathbf{K}}_i{\xi }_i+\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{\mathbf{K}}_{ij}{\xi }_i{\xi }_j+\cdots$$

$$\mathbf{U}(\mathbf{\xi })={\mathbf{U}}_0+\sum _{i=1}^n{\mathbf{U}}_i{\xi }_i+\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{\mathbf{U}}_{ij}{\xi }_i{\xi }_j+\cdots$$

代入平衡方程 $\mathbf{KU}=\mathbf{F}$，比较同阶项得到递推方程组。

3.4.2 谱随机有限元

将响应展开为多项式混沌形式，通过Galerkin投影得到增广的确定性方程组。这种方法精度高，但方程组规模随维度指数增长（维度灾难）。

四、可靠性分析

可靠性分析计算结构或系统在特定条件下完成预定功能的概率。核心概念是极限状态函数和失效概率。

4.1 极限状态与失效概率

4.1.1 极限状态函数

定义极限状态函数 $g(\mathbf{X})$，其中 $\mathbf{X}$ 为随机变量向量：

• $g(\mathbf{X})>0$：安全状态
• $g(\mathbf{X})=0$：极限状态（失效边界）
• $g(\mathbf{X})<0$：失效状态

失效概率为：

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

4.1.2 常用极限状态形式

强度-应力模型：$g=R-S$，其中 $R$ 为抗力，$S$ 为载荷效应。

变形准则：$g={\delta }_{allow}-\delta$，其中 ${\delta }_{allow}$ 为允许变形，$\delta$ 为实际变形。

疲劳寿命：$g=N_{design}-N$，其中 $N_{design}$ 为设计寿命，$N$ 为实际寿命。

4.2 一次二阶矩方法

一次二阶矩（First-Order Second-Moment, FOSM）方法用一阶近似计算可靠性指标。

4.2.1 均值一次二阶矩（MVFOSM）

在均值点线性化极限状态函数：

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

可靠性指标：

$$\beta =\frac{{\mu }_g}{{\sigma }_g}=\frac{g({\mathbf{\mu }}_X)}{\sqrt{{\sum }_{i=1}^n{\sum }_{j=1}^n\frac{\partial g}{\partial X_i}\frac{\partial g}{\partial X_j}{\rho }_{ij}{\sigma }_i{\sigma }_j}}$$

失效概率：$P_f=\Phi (-\beta )$，其中 $\Phi$ 为标准正态分布函数。

MVFOSM的缺点是不满足不变性（对功能函数的等价形式给出不同结果）。

4.2.2 改进一次二阶矩（AFOSM/FORM）

Hasofer-Lind可靠性指标在标准正态空间定义：

$$\beta =\underset{\mathbf{u}\in g(\mathbf{u})=0}{\min }\Vert \mathbf{u}\Vert$$

设计点（最可能失效点）为 ${\mathbf{u}}^{\ast }$，满足：

• 在失效边界上：$g({\mathbf{u}}^{\ast })=0$
• 距离原点最近

失效概率：$P_f\approx \Phi (-\beta )$

FORM算法步骤：

1. 将随机变量变换到标准正态空间：$\mathbf{u}=\mathbf{T}(\mathbf{x})$
2. 选择初始点（通常为均值点）
3. 在迭代点处线性化极限状态函数
4. 计算新的迭代点，向设计点收敛
5. 重复直到收敛，得到 $\beta$ 和 ${\mathbf{u}}^{\ast }$

4.2.3 二次二阶矩（SORM）

在设计点处进行二次近似，考虑极限状态曲率的影响：

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

其中，${\kappa }_i$ 为主曲率。SORM在高度非线性情况下比FORM更精确。

4.3 蒙特卡洛可靠性分析

蒙特卡洛方法通过直接抽样估计失效概率：

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

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

方差与样本量：

$$\text{Var}({\hat{P}}_f)=\frac{P_f(1-P_f)}{N}$$

对于小失效概率（如 $P_f=10^{-6}$），需要大量样本（$N\sim 10^8$）才能获得可靠估计。

重要性抽样：从偏置分布 $h(\mathbf{x})$ 抽样，估计量：

$${\hat{P}}_f=\frac{1}{N}\sum _{i=1}^{N}I(g({\mathbf{x}}^{(i)})\le 0)\frac{f({\mathbf{x}}^{(i)})}{h({\mathbf{x}}^{(i)})}$$

最优偏置分布以设计点为中心，可大幅减少所需样本量。

4.4 响应面方法

响应面方法用显式函数近似隐式极限状态函数，降低计算成本。

4.4.1 多项式响应面

常用二次多项式：

$$\tilde{g}(\mathbf{X})=a_0+\sum _{i=1}^n{a}_i{X}_i+\sum _{i=1}^n{b}_i{X}_i^2+\sum _{i<j}{c}_{ij}{X}_i{X}_j$$

通过在设计点附近的样本点拟合确定系数，然后用FORM/SORM计算可靠性。

4.4.2 神经网络响应面

对于高度非线性问题，可用神经网络建立输入-输出的映射关系，作为响应面进行可靠性分析。

五、灵敏度分析

灵敏度分析研究输入不确定性对输出不确定性的贡献程度，帮助识别关键参数。

5.1 局部灵敏度

5.1.1 偏导数灵敏度

输出对输入的偏导数：

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

无量纲化形式：

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

5.1.2 弹性系数

输入变化1%引起的输出变化百分比：

$$E_i=\frac{\partial Y}{\partial X_i}\frac{{\mu }_{X_i}}{{\mu }_Y}$$

5.2 全局灵敏度分析

全局灵敏度考虑输入在其整个分布范围内的影响。

5.2.1 方差分解与Sobol指标

将输出方差分解为各输入及其交互作用的贡献：

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

一阶Sobol指标（主效应）：

$$S_i=\frac{V_i}{\text{Var}(Y)}=\frac{{\text{Var}}_{X_i}(E_{{\mathbf{X}}_{\sim i}}[Y|X_i])}{\text{Var}(Y)}$$

总效应指标（包含所有交互作用）：

$$S_i^T=\frac{E_{{\mathbf{X}}_{\sim i}}[{\text{Var}}_{X_i}(Y|{\mathbf{X}}_{\sim i})]}{\text{Var}(Y)}=1-\frac{{\text{Var}}_{{\mathbf{X}}_{\sim i}}(E_{X_i}[Y|{\mathbf{X}}_{\sim i}])}{\text{Var}(Y)}$$

5.2.2 Morris方法

Morris方法通过计算基本效应的统计量进行筛选分析：

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

多次抽样计算基本效应的均值 ${\mu }_i$（总体影响）和标准差 ${\sigma }_i$（非线性/交互作用）。

5.3 可靠性灵敏度

可靠性灵敏度分析失效概率对分布参数的敏感性。

FORM灵敏度：

$$\frac{\partial P_f}{\partial {\mu }_{X_i}}=\frac{\varphi (-\beta )}{{\sigma }_g}\frac{\partial g}{\partial X_i}$$

$$\frac{\partial P_f}{\partial {\sigma }_{X_i}}=\frac{\varphi (-\beta )}{{\sigma }_g}\frac{\partial g}{\partial X_i}{u}_i^{\ast }$$

其中，${\mathbf{u}}^{\ast }$ 为设计点，$\varphi$ 为标准正态密度函数。

六、基于可靠性的优化设计

基于可靠性的优化设计（Reliability-Based Design Optimization, RBDO）在保证可靠性的前提下优化设计目标。

6.1 双层RBDO方法

外层优化：

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

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

其中，$\mathbf{d}$ 为设计变量，$\mathbf{X}$ 为随机变量。

内层可靠性分析：对每个候选设计点计算约束的失效概率。

双层方法计算成本高，但概念清晰。

6.2 单层RBDO方法

将可靠性约束转化为确定性约束，避免嵌套优化。

SORA方法（Sequential Optimization and Reliability Assessment）：

交替进行确定性优化和可靠性分析，通过偏移向量将可靠性约束转化为确定性约束：

$$g_i(\mathbf{d}+{\mathbf{s}}_i)\ge 0$$

其中，${\mathbf{s}}_i$ 为从均值点到设计点的偏移。

SLA方法（Single Loop Approach）：

将KKT条件纳入优化，同时求解设计变量和设计点。

6.3 鲁棒性设计优化

鲁棒性设计优化（Robust Design Optimization, RDO）最小化性能对不确定性的敏感性。

目标函数：

$$\underset{\mathbf{d}}{\min }{w}_1{\mu }_f(\mathbf{d})+w_2{\sigma }_f(\mathbf{d})$$

其中，${\mu }_f$ 和 ${\sigma }_f$ 分别为性能指标的均值和标准差，$w_1,w_2$ 为权重。

双目标表述：

min⁡d{μf(d),σf(d)}\min_{\mathbf{d}} \{\mu_f(\mathbf{d}), \sigma_f(\mathbf{d})\}dmin​{μf​(d),σf​(d)}

通过多目标优化获得Pareto前沿，权衡性能与鲁棒性。

七、应用案例

7.1 案例1：悬臂梁可靠性分析

考虑一端固定的悬臂梁，自由端受集中载荷。材料强度、载荷和几何尺寸存在不确定性，分析梁的弯曲失效可靠性。

7.2 案例2：复合材料层合板强度可靠性

考虑碳纤维/环氧树脂层合板，各层材料参数和铺层角度存在分散性，分析在面内载荷作用下的首层失效可靠性。

7.3 案例3：结构疲劳寿命预测

考虑承受循环载荷的机械零件，材料疲劳参数和载荷存在不确定性，预测疲劳寿命分布并计算达到设计寿命的可靠性。

7.4 案例4：基于可靠性的桁架结构优化

考虑三杆桁架结构，在满足强度和稳定性可靠性约束的前提下，最小化结构重量，获得最优截面尺寸。

八、总结

不确定性量化与可靠性分析为工程设计和风险评估提供了科学的理论框架和实用的计算工具。本章系统介绍了：

1. 不确定性建模：随机变量、随机过程、随机场、模糊与区间方法
2. 不确定性传播：蒙特卡洛模拟、泰勒展开、多项式混沌、随机有限元
3. 可靠性分析：FORM/SORM、蒙特卡洛、响应面方法
4. 灵敏度分析：局部与全局灵敏度、Sobol指标
5. 可靠性优化：RBDO、鲁棒性设计

随着计算能力的提升和算法的进步，不确定性量化与可靠性分析正成为工程仿真的标准组成部分。在实际应用中，需要根据问题特点选择合适的方法，在计算精度和效率之间取得平衡。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
不确定性量化与可靠性分析 - Python实现
包含4个案例：
1. 悬臂梁可靠性分析
2. 复合材料层合板强度可靠性
3. 结构疲劳寿命预测
4. 基于可靠性的桁架结构优化
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import minimize
from matplotlib.patches import FancyBboxPatch, Circle, Rectangle
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.size'] = 10

# ==================== 工具函数 ====================

def generate_latin_hypercube_samples(n_samples, n_vars, distributions):
    """生成拉丁超立方抽样样本"""
    samples = np.zeros((n_samples, n_vars))
    
    for j in range(n_vars):
        # 将每个维度划分为n_samples个等概率区间
        perm = np.random.permutation(n_samples)
        probs = (perm + np.random.rand(n_samples)) / n_samples
        
        # 根据分布类型生成样本
        dist_type, params = distributions[j]
        if dist_type == 'normal':
            mu, sigma = params
            samples[:, j] = stats.norm.ppf(probs, mu, sigma)
        elif dist_type == 'lognormal':
            mu, sigma = params
            samples[:, j] = stats.lognorm.ppf(probs, sigma, scale=np.exp(mu))
        elif dist_type == 'uniform':
            low, high = params
            samples[:, j] = stats.uniform.ppf(probs, low, high - low)
        elif dist_type == 'weibull':
            beta, eta = params
            samples[:, j] = stats.weibull_min.ppf(probs, beta, scale=eta)
    
    return samples

def monte_carlo_simulation(limit_state_func, samples):
    """蒙特卡洛模拟计算失效概率"""
    n_samples = len(samples)
    g_values = np.array([limit_state_func(sample) for sample in samples])
    
    n_failures = np.sum(g_values <= 0)
    pf = n_failures / n_samples
    
    # 计算置信区间
    variance = pf * (1 - pf) / n_samples
    std_error = np.sqrt(variance)
    ci_lower = pf - 1.96 * std_error
    ci_upper = pf + 1.96 * std_error
    
    return pf, std_error, (ci_lower, ci_upper), g_values

def form_analysis(limit_state_func, mean, std, corr=None, tol=1e-6, max_iter=100):
    """FORM可靠性分析"""
    n_vars = len(mean)
    
    # 初始化设计点（标准正态空间）
    u = np.zeros(n_vars)
    
    for iteration in range(max_iter):
        # 转换到原始空间
        x = mean + std * u
        
        # 计算极限状态函数值和梯度
        g_val = limit_state_func(x)
        
        # 数值计算梯度
        grad = np.zeros(n_vars)
        h = 1e-6
        for i in range(n_vars):
            x_plus = x.copy()
            x_plus[i] += h
            grad[i] = (limit_state_func(x_plus) - g_val) / h
        
        # 转换梯度到标准正态空间
        grad_u = grad * std
        
        # 计算新的设计点
        grad_norm = np.linalg.norm(grad_u)
        if grad_norm < 1e-10:
            break
        
        alpha = -grad_u / grad_norm
        beta_new = (g_val - np.dot(grad_u, u)) / grad_norm
        u_new = beta_new * alpha
        
        # 检查收敛
        if np.linalg.norm(u_new - u) < tol:
            u = u_new
            break
        
        u = u_new
    
    beta = np.linalg.norm(u)
    pf_form = stats.norm.cdf(-beta)
    
    return beta, pf_form, u, x

def sobol_indices(func, n_samples, n_vars, distributions):
    """计算Sobol灵敏度指标"""
    # 生成样本矩阵
    A = generate_latin_hypercube_samples(n_samples, n_vars, distributions)
    B = generate_latin_hypercube_samples(n_samples, n_vars, distributions)
    
    # 计算输出
    Y_A = np.array([func(a) for a in A])
    Y_B = np.array([func(b) for b in B])
    
    total_variance = np.var(np.concatenate([Y_A, Y_B]))
    
    first_order = np.zeros(n_vars)
    total_order = np.zeros(n_vars)
    
    for i in range(n_vars):
        # 创建混合矩阵
        C = B.copy()
        C[:, i] = A[:, i]
        Y_C = np.array([func(c) for c in C])
        
        # 一阶指标
        first_order[i] = np.mean(Y_B * (Y_C - Y_A)) / total_variance
        
        # 总效应指标
        total_order[i] = 0.5 * np.mean((Y_A - Y_C)**2) / total_variance
    
    return first_order, total_order

# ==================== 案例1：悬臂梁可靠性分析 ====================

def case1_cantilever_beam_reliability():
    """案例1：悬臂梁可靠性分析"""
    print("\n" + "="*60)
    print("案例1：悬臂梁可靠性分析")
    print("="*60)
    
    # 定义随机变量分布
    # [分布类型, 参数]
    distributions = [
        ('normal', [400e6, 40e6]),    # 屈服强度 (Pa), 均值400MPa, COV=10%
        ('normal', [1000, 200]),       # 载荷 (N), 均值1000N, COV=20%
        ('normal', [0.05, 0.0025]),    # 宽度 (m), 均值50mm, COV=5%
        ('normal', [0.01, 0.0005]),    # 高度 (m), 均值10mm, COV=5%
        ('normal', [1.0, 0.05]),       # 长度 (m), 均值1m, COV=5%
    ]
    
    def limit_state_strength(x):
        """强度极限状态函数"""
        sigma_y, P, b, h, L = x
        # 最大弯曲应力
        M_max = P * L  # 最大弯矩
        I = b * h**3 / 12  # 截面惯性矩
        y_max = h / 2  # 最大距离
        sigma_max = M_max * y_max / I
        return sigma_y - sigma_max
    
    def limit_state_deflection(x):
        """挠度极限状态函数"""
        sigma_y, P, b, h, L = x
        E = 200e9  # 弹性模量 (Pa)
        I = b * h**3 / 12
        delta_max = P * L**3 / (3 * E * I)  # 最大挠度
        delta_allow = L / 250  # 允许挠度
        return delta_allow - delta_max
    
    # 蒙特卡洛模拟
    print("\n[蒙特卡洛模拟 - 强度失效]")
    n_samples = 100000
    samples = generate_latin_hypercube_samples(n_samples, 5, distributions)
    
    pf_strength, se_strength, ci_strength, g_strength = monte_carlo_simulation(
        limit_state_strength, samples
    )
    
    print(f"  失效概率 Pf = {pf_strength:.6f}")
    print(f"  标准误差 = {se_strength:.6f}")
    print(f"  95%置信区间: [{ci_strength[0]:.6f}, {ci_strength[1]:.6f}]")
    
    print("\n[蒙特卡洛模拟 - 挠度失效]")
    pf_deflection, se_deflection, ci_deflection, g_deflection = monte_carlo_simulation(
        limit_state_deflection, samples
    )
    
    print(f"  失效概率 Pf = {pf_deflection:.6f}")
    print(f"  标准误差 = {se_deflection:.6f}")
    print(f"  95%置信区间: [{ci_deflection[0]:.6f}, {ci_deflection[1]:.6f}]")
    
    # FORM分析
    print("\n[FORM分析 - 强度失效]")
    mean = np.array([400e6, 1000, 0.05, 0.01, 1.0])
    std = np.array([40e6, 200, 0.0025, 0.0005, 0.05])
    
    beta, pf_form, u_star, x_star = form_analysis(limit_state_strength, mean, std)
    print(f"  可靠性指标 β = {beta:.4f}")
    print(f"  失效概率 Pf = {pf_form:.6f}")
    print(f"  设计点: σy={x_star[0]/1e6:.2f}MPa, P={x_star[1]:.2f}N")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 子图1: 悬臂梁示意图
    ax1 = plt.subplot(3, 3, 1)
    ax1.plot([0, 1], [0, 0], 'k-', linewidth=3)
    ax1.plot([0, 0], [-0.05, 0.05], 'k-', linewidth=8)
    ax1.arrow(1, 0, 0, -0.15, head_width=0.05, head_length=0.03, fc='red', ec='red')
    ax1.text(1.05, -0.08, 'P', fontsize=12, color='red', fontweight='bold')
    ax1.text(0.5, 0.02, 'L', fontsize=11, ha='center')
    ax1.text(-0.1, 0, 'Fixed', fontsize=9, ha='right')
    ax1.set_xlim(-0.2, 1.3)
    ax1.set_ylim(-0.25, 0.15)
    ax1.set_aspect('equal')
    ax1.axis('off')
    ax1.set_title('Cantilever Beam', fontsize=11, fontweight='bold')
    
    # 子图2: 极限状态函数分布
    ax2 = plt.subplot(3, 3, 2)
    ax2.hist(g_strength, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
    ax2.axvline(x=0, color='red', linestyle='--', linewidth=2, label='Failure boundary')
    ax2.set_xlabel('g(X) = σy - σmax', fontsize=10)
    ax2.set_ylabel('Probability Density', fontsize=10)
    ax2.set_title('Limit State Function Distribution', fontsize=11, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 子图3: 应力分布
    ax3 = plt.subplot(3, 3, 3)
    sigma_max_samples = []
    for sample in samples[:1000]:
        sigma_y, P, b, h, L = sample
        M_max = P * L
        I = b * h**3 / 12
        sigma_max = M_max * (h/2) / I
        sigma_max_samples.append(sigma_max / 1e6)
    
    ax3.hist(sigma_max_samples, bins=40, density=True, alpha=0.6, color='orange', label='Stress')
    ax3.axvline(x=400, color='blue', linestyle='--', linewidth=2, label='Mean Strength')
    ax3.set_xlabel('Maximum Stress (MPa)', fontsize=10)
    ax3.set_ylabel('Probability Density', fontsize=10)
    ax3.set_title('Stress Distribution', fontsize=11, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 子图4: 收敛性分析
    ax4 = plt.subplot(3, 3, 4)
    sample_sizes = np.logspace(2, 5, 20).astype(int)
    pf_convergence = []
    
    for n in sample_sizes:
        samples_n = generate_latin_hypercube_samples(n, 5, distributions)
        g_vals = np.array([limit_state_strength(s) for s in samples_n])
        pf_n = np.sum(g_vals <= 0) / n
        pf_convergence.append(pf_n)
    
    ax4.semilogx(sample_sizes, pf_convergence, 'b-o', markersize=4, linewidth=1.5)
    ax4.axhline(y=pf_strength, color='r', linestyle='--', linewidth=2, label='Converged')
    ax4.set_xlabel('Number of Samples', fontsize=10)
    ax4.set_ylabel('Failure Probability', fontsize=10)
    ax4.set_title('MC Convergence', fontsize=11, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 子图5: 参数敏感性（龙卷风图）
    ax5 = plt.subplot(3, 3, 5)
    var_names = ['σy', 'P', 'b', 'h', 'L']
    base_case = mean.copy()
    sensitivities = []
    
    for i in range(5):
        # 正向变化
        x_plus = base_case.copy()
        x_plus[i] += std[i]
        g_plus = limit_state_strength(x_plus)
        
        # 负向变化
        x_minus = base_case.copy()
        x_minus[i] -= std[i]
        g_minus = limit_state_strength(x_minus)
        
        sensitivity = (g_plus - g_minus) / 2
        sensitivities.append(sensitivity)
    
    colors = ['green' if s > 0 else 'red' for s in sensitivities]
    ax5.barh(var_names, sensitivities, color=colors, alpha=0.7, edgecolor='black')
    ax5.axvline(x=0, color='black', linewidth=1)
    ax5.set_xlabel('Sensitivity', fontsize=10)
    ax5.set_title('Parameter Sensitivity', fontsize=11, fontweight='bold')
    ax5.grid(True, alpha=0.3, axis='x')
    
    # 子图6: 可靠性指标随载荷变化
    ax6 = plt.subplot(3, 3, 6)
    P_range = np.linspace(500, 1500, 30)
    betas = []
    
    for P in P_range:
        mean_p = mean.copy()
        mean_p[1] = P
        
        def ls_p(x):
            x_mod = x.copy()
            x_mod[1] = P
            return limit_state_strength(x_mod)
        
        beta_p, _, _, _ = form_analysis(ls_p, mean_p, std)
        betas.append(beta_p)
    
    ax6.plot(P_range, betas, 'b-', linewidth=2)
    ax6.axhline(y=3.0, color='r', linestyle='--', linewidth=2, label='Target β=3.0')
    ax6.set_xlabel('Load P (N)', fontsize=10)
    ax6.set_ylabel('Reliability Index β', fontsize=10)
    ax6.set_title('β vs Load', fontsize=11, fontweight='bold')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    # 子图7: 散点图 - 强度vs应力
    ax7 = plt.subplot(3, 3, 7)
    sigma_y_samples = samples[:2000, 0] / 1e6
    sigma_max_samples_2000 = []
    for sample in samples[:2000]:
        _, P, b, h, L = sample
        M_max = P * L
        I = b * h**3 / 12
        sigma_max = M_max * (h/2) / I
        sigma_max_samples_2000.append(sigma_max / 1e6)
    
    failure_mask = np.array(sigma_max_samples_2000) > sigma_y_samples
    ax7.scatter(sigma_y_samples[~failure_mask], np.array(sigma_max_samples_2000)[~failure_mask], 
               c='blue', alpha=0.5, s=10, label='Safe')
    ax7.scatter(sigma_y_samples[failure_mask], np.array(sigma_max_samples_2000)[failure_mask], 
               c='red', alpha=0.5, s=10, label='Failure')
    ax7.plot([200, 600], [200, 600], 'k--', linewidth=2, label='g=0')
    ax7.set_xlabel('Yield Strength (MPa)', fontsize=10)
    ax7.set_ylabel('Max Stress (MPa)', fontsize=10)
    ax7.set_title('Strength vs Stress', fontsize=11, fontweight='bold')
    ax7.legend()
    ax7.grid(True, alpha=0.3)
    
    # 子图8: 失效模式对比
    ax8 = plt.subplot(3, 3, 8)
    modes = ['Strength\nFailure', 'Deflection\nFailure', 'System\nFailure']
    pf_values = [pf_strength, pf_deflection, 1 - (1-pf_strength)*(1-pf_deflection)]
    colors_bar = ['red', 'orange', 'purple']
    
    bars = ax8.bar(modes, pf_values, color=colors_bar, alpha=0.7, edgecolor='black')
    ax8.set_ylabel('Failure Probability', fontsize=10)
    ax8.set_title('Failure Mode Comparison', fontsize=11, fontweight='bold')
    ax8.grid(True, alpha=0.3, axis='y')
    
    for bar, val in zip(bars, pf_values):
        height = bar.get_height()
        ax8.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.4f}', ha='center', va='bottom', fontsize=9, fontweight='bold')
    
    # 子图9: 参数统计表
    ax9 = plt.subplot(3, 3, 9)
    ax9.axis('off')
    
    table_data = [
        ['Parameter', 'Mean', 'COV (%)', 'Dist.'],
        ['Yield Strength', '400 MPa', '10', 'Normal'],
        ['Load P', '1000 N', '20', 'Normal'],
        ['Width b', '50 mm', '5', 'Normal'],
        ['Height h', '10 mm', '5', 'Normal'],
        ['Length L', '1.0 m', '5', 'Normal'],
        ['', '', '', ''],
        ['Results', '', '', ''],
        ['Pf (Strength)', f'{pf_strength:.6f}', '', ''],
        ['Pf (Deflection)', f'{pf_deflection:.6f}', '', ''],
        ['β (FORM)', f'{beta:.3f}', '', ''],
    ]
    
    table = ax9.table(cellText=table_data, cellLoc='center',
                     colWidths=[0.25, 0.25, 0.25, 0.25],
                     loc='center', bbox=[0, 0, 1, 1])
    table.auto_set_font_size(False)
    table.set_fontsize(8)
    table.scale(1, 1.5)
    
    for i in range(4):
        table[(0, i)].set_facecolor('#4472C4')
        table[(0, i)].set_text_props(weight='bold', color='white')
    
    for row in [7]:
        for i in range(4):
            table[(row, i)].set_facecolor('#E7E6E6')
            table[(row, i)].set_text_props(weight='bold')
    
    ax9.set_title('Parameters & Results', fontsize=11, fontweight='bold', pad=20)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题041\\output\\case1_cantilever_beam.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例1完成，图像已保存")
    print(f"  - 强度失效概率: {pf_strength:.6f}")
    print(f"  - 挠度失效概率: {pf_deflection:.6f}")
    print(f"  - FORM可靠性指标: {beta:.4f}")

# ==================== 案例2：复合材料层合板强度可靠性 ====================

def case2_composite_laminate_reliability():
    """案例2：复合材料层合板强度可靠性"""
    print("\n" + "="*60)
    print("案例2：复合材料层合板强度可靠性")
    print("="*60)
    
    # Tsai-Wu失效准则
    def tsai_wu_criterion(sigma_1, sigma_2, tau_12, F1, F2, F11, F22, F66, F12):
        """Tsai-Wu失效准则"""
        term1 = F1 * sigma_1 + F2 * sigma_2
        term2 = F11 * sigma_1**2 + F22 * sigma_2**2 + F66 * tau_12**2
        term3 = 2 * F12 * sigma_1 * sigma_2
        return term1 + term2 + term3
    
    # 材料参数（碳纤维/环氧树脂）
    # 分布: [类型, 参数]
    distributions = [
        ('normal', [1500e6, 150e6]),   # X_t (纵向拉伸强度, Pa)
        ('normal', [1200e6, 120e6]),   # X_c (纵向压缩强度, Pa)
        ('normal', [50e6, 5e6]),       # Y_t (横向拉伸强度, Pa)
        ('normal', [200e6, 20e6]),     # Y_c (横向压缩强度, Pa)
        ('normal', [70e6, 7e6]),       # S (剪切强度, Pa)
        ('normal', [140e9, 7e9]),      # E1 (纵向模量, Pa)
        ('normal', [10e9, 0.5e9]),     # E2 (横向模量, Pa)
        ('normal', [0.3, 0.015]),      # nu12 (泊松比)
        ('normal', [5e9, 0.25e9]),     # G12 (剪切模量, Pa)
    ]
    
    # 铺层角度（随机分散）
    ply_angles = [0, 45, -45, 90, 90, -45, 45, 0]  # 对称铺层
    n_plies = len(ply_angles)
    ply_thickness = 0.125e-3  # 单层厚度 0.125mm
    laminate_thickness = n_plies * ply_thickness
    
    def calculate_laminate_stresses(Nx, Ny, Nxy, material_props):
        """计算层合板各层应力"""
        X_t, X_c, Y_t, Y_c, S, E1, E2, nu12, G12 = material_props
        
        # 计算各层应力（简化模型）
        stresses = []
        for theta in ply_angles:
            theta_rad = np.radians(theta)
            
            # 转换载荷到材料主方向（简化）
            sigma_x = Nx / laminate_thickness
            sigma_y = Ny / laminate_thickness
            tau_xy = Nxy / laminate_thickness
            
            # 应力转换
            c2 = np.cos(theta_rad)**2
            s2 = np.sin(theta_rad)**2
            cs = np.cos(theta_rad) * np.sin(theta_rad)
            
            sigma_1 = c2 * sigma_x + s2 * sigma_y + 2 * cs * tau_xy
            sigma_2 = s2 * sigma_x + c2 * sigma_y - 2 * cs * tau_xy
            tau_12 = -cs * sigma_x + cs * sigma_y + (c2 - s2) * tau_xy
            
            stresses.append([sigma_1, sigma_2, tau_12])
        
        return stresses
    
    def limit_state_composite(x):
        """复合材料层合板极限状态函数"""
        # 材料参数
        material_props = x[:9]
        X_t, X_c, Y_t, Y_c, S, E1, E2, nu12, G12 = material_props
        
        # 载荷（最后三个参数）
        Nx = 500e3  # N/m
        Ny = 100e3  # N/m
        Nxy = 50e3  # N/m
        
        # Tsai-Wu系数
        F1 = 1/X_t - 1/X_c
        F2 = 1/Y_t - 1/Y_c
        F11 = 1/(X_t * X_c)
        F22 = 1/(Y_t * Y_c)
        F66 = 1/S**2
        F12 = -0.5 * np.sqrt(F11 * F22)
        
        # 计算各层应力
        stresses = calculate_laminate_stresses(Nx, Ny, Nxy, material_props)
        
        # 找出最大失效指数
        max_failure_index = 0
        for sigma_1, sigma_2, tau_12 in stresses:
            failure_index = tsai_wu_criterion(sigma_1, sigma_2, tau_12, 
                                             F1, F2, F11, F22, F66, F12)
            max_failure_index = max(max_failure_index, failure_index)
        
        # 极限状态：1 - 失效指数
        return 1.0 - max_failure_index
    
    # 蒙特卡洛模拟
    print("\n[蒙特卡洛模拟]")
    n_samples = 50000
    samples = generate_latin_hypercube_samples(n_samples, 9, distributions)
    
    pf, se, ci, g_values = monte_carlo_simulation(limit_state_composite, samples)
    
    print(f"  首层失效概率 Pf = {pf:.6f}")
    print(f"  标准误差 = {se:.6f}")
    print(f"  95%置信区间: [{ci[0]:.6f}, {ci[1]:.6f}]")
    
    # FORM分析
    print("\n[FORM分析]")
    mean = np.array([1500e6, 1200e6, 50e6, 200e6, 70e6, 140e9, 10e9, 0.3, 5e9])
    std = np.array([150e6, 120e6, 5e6, 20e6, 7e6, 7e9, 0.5e9, 0.015, 0.25e9])
    
    beta, pf_form, u_star, x_star = form_analysis(limit_state_composite, mean, std)
    print(f"  可靠性指标 β = {beta:.4f}")
    print(f"  失效概率 Pf = {pf_form:.6f}")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 子图1: 层合板铺层示意图
    ax1 = plt.subplot(3, 3, 1)
    colors_ply = plt.cm.rainbow(np.linspace(0, 1, n_plies))
    for i, (angle, color) in enumerate(zip(ply_angles, colors_ply)):
        rect = Rectangle((0, i*0.5), 2, 0.5, facecolor=color, edgecolor='black', alpha=0.7)
        ax1.add_patch(rect)
        ax1.text(1, i*0.5 + 0.25, f'{angle}°', ha='center', va='center', 
                fontsize=10, fontweight='bold')
    
    ax1.set_xlim(-0.2, 2.2)
    ax1.set_ylim(-0.2, n_plies*0.5 + 0.2)
    ax1.set_aspect('equal')
    ax1.axis('off')
    ax1.set_title('Laminate Stacking Sequence', fontsize=11, fontweight='bold')
    
    # 子图2: 极限状态函数分布
    ax2 = plt.subplot(3, 3, 2)
    ax2.hist(g_values, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
    ax2.axvline(x=0, color='red', linestyle='--', linewidth=2, label='Failure boundary')
    ax2.set_xlabel('g(X) = 1 - Failure Index', fontsize=10)
    ax2.set_ylabel('Probability Density', fontsize=10)
    ax2.set_title('Limit State Distribution', fontsize=11, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 子图3: 各层失效指数分布
    ax3 = plt.subplot(3, 3, 3)
    failure_indices_by_ply = [[] for _ in range(n_plies)]
    
    for sample in samples[:2000]:
        material_props = sample[:9]
        X_t, X_c, Y_t, Y_c, S = material_props[:5]
        
        F1 = 1/X_t - 1/X_c
        F2 = 1/Y_t - 1/Y_c
        F11 = 1/(X_t * X_c)
        F22 = 1/(Y_t * Y_c)
        F66 = 1/S**2
        F12 = -0.5 * np.sqrt(F11 * F22)
        
        Nx, Ny, Nxy = 500e3, 100e3, 50e3
        stresses = calculate_laminate_stresses(Nx, Ny, Nxy, material_props)
        
        for i, (sigma_1, sigma_2, tau_12) in enumerate(stresses):
            fi = tsai_wu_criterion(sigma_1, sigma_2, tau_12, F1, F2, F11, F22, F66, F12)
            failure_indices_by_ply[i].append(fi)
    
    mean_fi = [np.mean(fi_list) for fi_list in failure_indices_by_ply]
    std_fi = [np.std(fi_list) for fi_list in failure_indices_by_ply]
    
    x_pos = np.arange(n_plies)
    ax3.bar(x_pos, mean_fi, yerr=std_fi, capsize=3, color='skyblue', 
           edgecolor='black', alpha=0.7)
    ax3.axhline(y=1.0, color='red', linestyle='--', linewidth=2, label='Failure threshold')
    ax3.set_xlabel('Ply Number', fontsize=10)
    ax3.set_ylabel('Failure Index', fontsize=10)
    ax3.set_title('Failure Index by Ply', fontsize=11, fontweight='bold')
    ax3.set_xticks(x_pos)
    ax3.set_xticklabels([f'{i+1}\n({ply_angles[i]}°)' for i in range(n_plies)], fontsize=8)
    ax3.legend()
    ax3.grid(True, alpha=0.3, axis='y')
    
    # 子图4: 材料强度分布
    ax4 = plt.subplot(3, 3, 4)
    X_t_samples = samples[:2000, 0] / 1e6
    Y_t_samples = samples[:2000, 2] / 1e6
    
    ax4.hist(X_t_samples, bins=30, density=True, alpha=0.6, color='blue', label='X_t (0°)')
    ax4.hist(Y_t_samples, bins=30, density=True, alpha=0.6, color='red', label='Y_t (90°)')
    ax4.set_xlabel('Strength (MPa)', fontsize=10)
    ax4.set_ylabel('Probability Density', fontsize=10)
    ax4.set_title('Strength Distribution', fontsize=11, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 子图5: 载荷-失效概率曲线
    ax5 = plt.subplot(3, 3, 5)
    Nx_range = np.linspace(300e3, 800e3, 15)
    pf_vs_load = []
    
    for Nx in Nx_range:
        def ls_load(x):
            material_props = x[:9]
            X_t, X_c, Y_t, Y_c, S = material_props[:5]
            
            F1 = 1/X_t - 1/X_c
            F2 = 1/Y_t - 1/Y_c
            F11 = 1/(X_t * X_c)
            F22 = 1/(Y_t * Y_c)
            F66 = 1/S**2
            F12 = -0.5 * np.sqrt(F11 * F22)
            
            stresses = calculate_laminate_stresses(Nx, 100e3, 50e3, material_props)
            
            max_fi = 0
            for sigma_1, sigma_2, tau_12 in stresses:
                fi = tsai_wu_criterion(sigma_1, sigma_2, tau_12, F1, F2, F11, F22, F66, F12)
                max_fi = max(max_fi, fi)
            
            return 1.0 - max_fi
        
        g_vals = np.array([ls_load(s) for s in samples[:10000]])
        pf_n = np.sum(g_vals <= 0) / len(g_vals)
        pf_vs_load.append(pf_n)
    
    ax5.semilogy(Nx_range/1e3, pf_vs_load, 'b-o', markersize=5, linewidth=2)
    ax5.set_xlabel('Nx (kN/m)', fontsize=10)
    ax5.set_ylabel('Failure Probability', fontsize=10)
    ax5.set_title('Pf vs Load', fontsize=11, fontweight='bold')
    ax5.grid(True, alpha=0.3)
    
    # 子图6: 安全系数直方图
    ax6 = plt.subplot(3, 3, 6)
    safety_factors = 1 / np.array([np.mean(fi_list) for fi_list in failure_indices_by_ply])
    
    ax6.bar(range(n_plies), safety_factors, color='green', alpha=0.7, edgecolor='black')
    ax6.axhline(y=1.5, color='red', linestyle='--', linewidth=2, label='Target SF=1.5')
    ax6.set_xlabel('Ply Number', fontsize=10)
    ax6.set_ylabel('Safety Factor', fontsize=10)
    ax6.set_title('Safety Factor by Ply', fontsize=11, fontweight='bold')
    ax6.set_xticks(range(n_plies))
    ax6.set_xticklabels([f'{i+1}' for i in range(n_plies)])
    ax6.legend()
    ax6.grid(True, alpha=0.3, axis='y')
    
    # 子图7: Tsai-Wu准则可视化
    ax7 = plt.subplot(3, 3, 7)
    sigma_1_range = np.linspace(-1500, 1500, 100)
    sigma_2_range = np.linspace(-300, 300, 100)
    S1, S2 = np.meshgrid(sigma_1_range, sigma_2_range)
    
    # 使用平均材料参数
    X_t_m, X_c_m, Y_t_m, Y_c_m, S_m = 1500, 1200, 50, 200, 70
    F1_m = 1/X_t_m - 1/X_c_m
    F2_m = 1/Y_t_m - 1/Y_c_m
    F11_m = 1/(X_t_m * X_c_m)
    F22_m = 1/(Y_t_m * Y_c_m)
    F66_m = 1/S_m**2
    F12_m = -0.5 * np.sqrt(F11_m * F22_m)
    
    TW = (F1_m*S1 + F2_m*S2 + F11_m*S1**2 + F22_m*S2**2 + 
          2*F12_m*S1*S2)
    
    contour = ax7.contour(S1, S2, TW, levels=[0.5, 0.8, 1.0, 1.2], colors='blue')
    ax7.clabel(contour, inline=True, fontsize=8)
    ax7.contourf(S1, S2, TW, levels=[1.0, 2.0], colors='red', alpha=0.3)
    ax7.set_xlabel('σ₁ (MPa)', fontsize=10)
    ax7.set_ylabel('σ₂ (MPa)', fontsize=10)
    ax7.set_title('Tsai-Wu Failure Criterion', fontsize=11, fontweight='bold')
    ax7.grid(True, alpha=0.3)
    
    # 子图8: 可靠性方法对比
    ax8 = plt.subplot(3, 3, 8)
    methods = ['MCS\n(50000)', 'FORM', 'MVFOSM']
    pf_methods = [pf, pf_form, pf]  # MVFOSM近似用MCS结果
    beta_methods = [stats.norm.ppf(1-pf), beta, stats.norm.ppf(1-pf)]
    
    x_pos = np.arange(len(methods))
    ax8_twin = ax8.twinx()
    
    bars1 = ax8.bar(x_pos - 0.2, pf_methods, 0.4, color='red', alpha=0.7, label='Pf')
    bars2 = ax8_twin.bar(x_pos + 0.2, beta_methods, 0.4, color='blue', alpha=0.7, label='β')
    
    ax8.set_xlabel('Method', fontsize=10)
    ax8.set_ylabel('Failure Probability', fontsize=10, color='red')
    ax8_twin.set_ylabel('Reliability Index β', fontsize=10, color='blue')
    ax8.set_title('Method Comparison', fontsize=11, fontweight='bold')
    ax8.set_xticks(x_pos)
    ax8.set_xticklabels(methods)
    ax8.grid(True, alpha=0.3, axis='y')
    
    # 子图9: 参数重要性排序
    ax9 = plt.subplot(3, 3, 9)
    ax9.axis('off')
    
    param_names = ['X_t', 'X_c', 'Y_t', 'Y_c', 'S', 'E1', 'E2', 'ν12', 'G12']
    
    # 计算FORM灵敏度（简化）
    sensitivities_form = np.abs(u_star) / np.sum(np.abs(u_star))
    
    table_data = [['Parameter', 'Sensitivity', 'Rank']]
    sorted_indices = np.argsort(sensitivities_form)[::-1]
    for rank, idx in enumerate(sorted_indices[:5], 1):
        table_data.append([param_names[idx], f'{sensitivities_form[idx]:.3f}', str(rank)])
    
    table = ax9.table(cellText=table_data, cellLoc='center',
                     colWidths=[0.4, 0.3, 0.3],
                     loc='center', bbox=[0, 0, 1, 1])
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1, 1.5)
    
    for i in range(3):
        table[(0, i)].set_facecolor('#4472C4')
        table[(0, i)].set_text_props(weight='bold', color='white')
    
    ax9.set_title('Parameter Importance (FORM)', fontsize=11, fontweight='bold', pad=20)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题041\\output\\case2_composite_laminate.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例2完成，图像已保存")
    print(f"  - 首层失效概率: {pf:.6f}")
    print(f"  - FORM可靠性指标: {beta:.4f}")

# ==================== 案例3：结构疲劳寿命预测 ====================

def case3_fatigue_life_prediction():
    """案例3：结构疲劳寿命预测"""
    print("\n" + "="*60)
    print("案例3：结构疲劳寿命预测")
    print("="*60)
    
    # S-N曲线参数（对数正态分布）
    # 分布: [类型, 参数]
    distributions = [
        ('lognormal', [12.5, 0.15]),   # log(C) - S-N曲线截距
        ('normal', [3.0, 0.1]),        # m - S-N曲线斜率
        ('lognormal', [200e6, 20e6]),  # σ_f' - 疲劳强度系数
        ('lognormal', [0.6, 0.06]),    # ε_f' - 疲劳延性系数
        ('normal', [0.1, 0.01]),       # b - 疲劳强度指数
        ('normal', [0.5, 0.05]),       # c - 疲劳延性指数
        ('lognormal', [100e6, 10e6]),  # E - 弹性模量
    ]
    
    # 载荷谱参数
    sigma_max = 300e6  # 最大应力 (Pa)
    sigma_min = 50e6   # 最小应力 (Pa)
    stress_ratio = sigma_min / sigma_max
    
    def calculate_fatigue_life(x, stress_amplitude):
        """计算疲劳寿命（Basquin方程）"""
        logC, m, sigma_f_prime, epsilon_f_prime, b, c, E = x
        C = np.exp(logC)
        
        # Basquin方程: N = C * (σa)^(-m)
        # 限制m的范围避免数值问题
        m_clamped = np.clip(m, 2.0, 10.0)
        N = C * (stress_amplitude / 1e6) ** (-m_clamped)
        
        # 限制最大寿命避免无穷大
        N = np.clip(N, 1, 1e12)
        
        return N
    
    def limit_state_fatigue(x):
        """疲劳寿命极限状态函数"""
        stress_amplitude = (sigma_max - sigma_min) / 2
        N = calculate_fatigue_life(x, stress_amplitude)
        N_design = 1e6  # 设计寿命
        
        return N - N_design
    
    # 蒙特卡洛模拟
    print("\n[蒙特卡洛模拟]")
    n_samples = 100000
    samples = generate_latin_hypercube_samples(n_samples, 7, distributions)
    
    pf, se, ci, g_values = monte_carlo_simulation(limit_state_fatigue, samples)
    
    print(f"  失效概率 Pf = {pf:.6f}")
    print(f"  标准误差 = {se:.6f}")
    print(f"  95%置信区间: [{ci[0]:.6f}, {ci[1]:.6f}]")
    
    # 计算寿命分布
    stress_amplitude = (sigma_max - sigma_min) / 2
    life_samples = []
    for sample in samples:
        N = calculate_fatigue_life(sample, stress_amplitude)
        life_samples.append(N)
    
    life_samples = np.array(life_samples)
    
    print(f"\n[疲劳寿命统计]")
    print(f"  平均寿命: {np.mean(life_samples):.2e}  cycles")
    print(f"  寿命中位数: {np.median(life_samples):.2e}  cycles")
    print(f"  寿命标准差: {np.std(life_samples):.2e}  cycles")
    print(f"  达到设计寿命(1e6)的可靠性: {np.sum(life_samples > 1e6)/len(life_samples)*100:.2f}%")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 子图1: S-N曲线
    ax1 = plt.subplot(3, 3, 1)
    stress_range = np.linspace(100e6, 500e6, 100)
    
    # 计算多条S-N曲线（基于随机样本）
    for i in range(50):
        sample = samples[i]
        logC, m = sample[0], sample[1]
        C = np.exp(logC)
        N_curve = C * (stress_range / 1e6) ** (-m)
        ax1.loglog(N_curve, stress_range/1e6, 'b-', alpha=0.2, linewidth=0.5)
    
    # 平均曲线
    mean_logC, mean_m = 12.5, 3.0
    mean_C = np.exp(mean_logC)
    N_mean = mean_C * (stress_range / 1e6) ** (-mean_m)
    ax1.loglog(N_mean, stress_range/1e6, 'r-', linewidth=3, label='Mean S-N')
    
    ax1.axhline(y=stress_amplitude/1e6, color='green', linestyle='--', linewidth=2, label='Operating stress')
    ax1.axvline(x=1e6, color='orange', linestyle='--', linewidth=2, label='Design life')
    ax1.set_xlabel('Cycles to Failure N', fontsize=10)
    ax1.set_ylabel('Stress Amplitude (MPa)', fontsize=10)
    ax1.set_title('S-N Curve Scatter', fontsize=11, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3, which='both')
    
    # 子图2: 寿命分布直方图
    ax2 = plt.subplot(3, 3, 2)
    log_life = np.log10(life_samples)
    ax2.hist(log_life, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')
    ax2.axvline(x=6, color='red', linestyle='--', linewidth=2, label='Design life (1e6)')
    ax2.axvline(x=np.mean(log_life), color='green', linestyle='--', linewidth=2, label=f'Mean ({np.mean(life_samples):.2e})')
    ax2.set_xlabel('log₁₀(N)', fontsize=10)
    ax2.set_ylabel('Probability Density', fontsize=10)
    ax2.set_title('Fatigue Life Distribution', fontsize=11, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 子图3: 累积分布函数
    ax3 = plt.subplot(3, 3, 3)
    sorted_life = np.sort(life_samples)
    cdf = np.arange(1, len(sorted_life)+1) / len(sorted_life)
    
    ax3.semilogx(sorted_life, cdf, 'b-', linewidth=2)
    ax3.axvline(x=1e6, color='red', linestyle='--', linewidth=2, label='Design life')
    ax3.axhline(y=0.5, color='green', linestyle='--', linewidth=2, label='Median')
    ax3.set_xlabel('Cycles to Failure N', fontsize=10)
    ax3.set_ylabel('Cumulative Probability', fontsize=10)
    ax3.set_title('CDF of Fatigue Life', fontsize=11, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 子图4: 可靠性随应力幅值变化
    ax4 = plt.subplot(3, 3, 4)
    stress_amplitudes = np.linspace(150e6, 400e6, 20)
    reliabilities = []
    
    for sa in stress_amplitudes:
        lives = []
        for sample in samples[:10000]:
            N = calculate_fatigue_life(sample, sa)
            lives.append(N)
        
        rel = np.sum(np.array(lives) > 1e6) / len(lives)
        reliabilities.append(rel)
    
    ax4.plot(stress_amplitudes/1e6, reliabilities, 'b-o', markersize=5, linewidth=2)
    ax4.axhline(y=0.99, color='red', linestyle='--', linewidth=2, label='Target 99%')
    ax4.axvline(x=stress_amplitude/1e6, color='green', linestyle='--', linewidth=2, label='Current')
    ax4.set_xlabel('Stress Amplitude (MPa)', fontsize=10)
    ax4.set_ylabel('Reliability', fontsize=10)
    ax4.set_title('Reliability vs Stress', fontsize=11, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 子图5: 参数敏感性（Sobol指标）
    ax5 = plt.subplot(3, 3, 5)
    param_names = ['log(C)', 'm', "σf'", "εf'", 'b', 'c', 'E']
    
    # 简化计算一阶Sobol指标
    sobol_first = np.zeros(7)
    total_var = np.var(life_samples)
    
    for i in range(7):
        # 条件期望方差
        sorted_indices = np.argsort(samples[:, i])
        n_bins = 10
        bin_size = len(samples) // n_bins
        
        conditional_means = []
        for j in range(n_bins):
            start = j * bin_size
            end = (j + 1) * bin_size
            bin_indices = sorted_indices[start:end]
            conditional_means.append(np.mean(life_samples[bin_indices]))
        
        var_cond_mean = np.var(conditional_means)
        sobol_first[i] = var_cond_mean / total_var
    
    colors = plt.cm.RdYlGn(np.linspace(0.2, 0.8, 7))
    ax5.barh(param_names, sobol_first, color=colors, edgecolor='black', alpha=0.8)
    ax5.set_xlabel('First-order Sobol Index', fontsize=10)
    ax5.set_title('Parameter Sensitivity', fontsize=11, fontweight='bold')
    ax5.grid(True, alpha=0.3, axis='x')
    
    # 子图6: 威布尔拟合
    ax6 = plt.subplot(3, 3, 6)
    # 拟合威布尔分布
    life_positive = life_samples[life_samples > 0]
    beta_weibull, eta_weibull, _ = stats.weibull_min.fit(life_positive, floc=0)
    
    x_weibull = np.linspace(np.min(life_positive), np.percentile(life_positive, 99), 100)
    pdf_weibull = stats.weibull_min.pdf(x_weibull, beta_weibull, scale=eta_weibull)
    
    ax6.hist(life_positive, bins=50, density=True, alpha=0.6, color='blue', label='Data')
    ax6.plot(x_weibull, pdf_weibull, 'r-', linewidth=2, label=f'Weibull fit (β={beta_weibull:.2f})')
    ax6.set_xlabel('Cycles to Failure', fontsize=10)
    ax6.set_ylabel('Probability Density', fontsize=10)
    ax6.set_title('Weibull Distribution Fit', fontsize=11, fontweight='bold')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    # 子图7: 载荷谱
    ax7 = plt.subplot(3, 3, 7)
    cycles = np.linspace(0, 1000, 1000)
    stress_history = sigma_max * (0.5 + 0.5 * np.sin(2 * np.pi * cycles / 100))
    
    ax7.plot(cycles, stress_history/1e6, 'b-', linewidth=1)
    ax7.axhline(y=sigma_max/1e6, color='red', linestyle='--', alpha=0.5, label='σmax')
    ax7.axhline(y=sigma_min/1e6, color='red', linestyle='--', alpha=0.5, label='σmin')
    ax7.fill_between(cycles, stress_history/1e6, alpha=0.3)
    ax7.set_xlabel('Cycles', fontsize=10)
    ax7.set_ylabel('Stress (MPa)', fontsize=10)
    ax7.set_title('Load Spectrum', fontsize=11, fontweight='bold')
    ax7.legend()
    ax7.grid(True, alpha=0.3)
    
    # 子图8: 损伤累积
    ax8 = plt.subplot(3, 3, 8)
    # Miner线性累积损伤
    n_cycles = np.linspace(0, 2e6, 100)
    D_median = n_cycles / np.median(life_samples)
    D_95 = n_cycles / np.percentile(life_samples, 5)
    D_5 = n_cycles / np.percentile(life_samples, 95)
    
    ax8.plot(n_cycles/1e6, D_median, 'b-', linewidth=2, label='Median')
    ax8.fill_between(n_cycles/1e6, D_5, D_95, alpha=0.3, color='blue', label='90% CI')
    ax8.axhline(y=1.0, color='red', linestyle='--', linewidth=2, label='Failure (D=1)')
    ax8.set_xlabel('Applied Cycles (millions)', fontsize=10)
    ax8.set_ylabel('Damage D', fontsize=10)
    ax8.set_title('Cumulative Damage', fontsize=11, fontweight='bold')
    ax8.legend()
    ax8.grid(True, alpha=0.3)
    
    # 子图9: 结果汇总表
    ax9 = plt.subplot(3, 3, 9)
    ax9.axis('off')
    
    table_data = [
        ['Parameter', 'Value', 'Unit'],
        ['', '', ''],
        ['σmax', f'{sigma_max/1e6:.0f}', 'MPa'],
        ['σmin', f'{sigma_min/1e6:.0f}', 'MPa'],
        ['Stress amplitude', f'{stress_amplitude/1e6:.0f}', 'MPa'],
        ['Design life', '1.0e6', 'cycles'],
        ['', '', ''],
        ['Results', '', ''],
        ['Mean life', f'{np.mean(life_samples):.2e}', 'cycles'],
        ['Median life', f'{np.median(life_samples):.2e}', 'cycles'],
        ['Pf (1e6)', f'{pf:.6f}', '-'],
        ['Reliability', f'{(1-pf)*100:.2f}', '%'],
        ['Weibull β', f'{beta_weibull:.2f}', '-'],
    ]
    
    table = ax9.table(cellText=table_data, cellLoc='center',
                     colWidths=[0.4, 0.35, 0.25],
                     loc='center', bbox=[0, 0, 1, 1])
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1, 1.5)
    
    for i in range(3):
        table[(0, i)].set_facecolor('#4472C4')
        table[(0, i)].set_text_props(weight='bold', color='white')
    
    for row in [7]:
        for i in range(3):
            table[(row, i)].set_facecolor('#E7E6E6')
            table[(row, i)].set_text_props(weight='bold')
    
    ax9.set_title('Fatigue Analysis Results', fontsize=11, fontweight='bold', pad=20)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题041\\output\\case3_fatigue_life.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例3完成，图像已保存")
    print(f"  - 失效概率: {pf:.6f}")
    print(f"  - 平均疲劳寿命: {np.mean(life_samples):.2e} cycles")
    print(f"  - 威布尔形状参数: {beta_weibull:.2f}")

# ==================== 案例4：基于可靠性的桁架结构优化 ====================

def case4_rbdo_truss_optimization():
    """案例4：基于可靠性的桁架结构优化"""
    print("\n" + "="*60)
    print("案例4：基于可靠性的桁架结构优化")
    print("="*60)
    
    # 三杆桁架结构
    # 节点坐标
    nodes = np.array([
        [0, 0],      # 节点1 (固定)
        [1, 0],      # 节点2 (固定)
        [0.5, 0.866] # 节点3 (自由，载荷作用点)
    ])
    
    # 杆件连接
    elements = [
        [0, 2],  # 杆件1: 节点1-3
        [1, 2],  # 杆件2: 节点2-3
        [0, 1]   # 杆件3: 节点1-2 (水平杆)
    ]
    
    n_elements = len(elements)
    
    # 材料参数（随机）
    E_mean = 200e9  # 弹性模量 (Pa)
    sigma_y_mean = 250e6  # 屈服强度 (Pa)
    
    # 载荷（随机）
    P_vertical_mean = 100e3  # 垂直载荷 (N)
    P_horizontal_mean = 50e3  # 水平载荷 (N)
    
    def truss_analysis(areas, E, P_v, P_h):
        """桁架结构分析"""
        # 计算杆件长度和方向
        lengths = []
        directions = []
        
        for elem in elements:
            node_i, node_j = elem
            dx = nodes[node_j, 0] - nodes[node_i, 0]
            dy = nodes[node_j, 1] - nodes[node_i, 1]
            length = np.sqrt(dx**2 + dy**2)
            lengths.append(length)
            directions.append([dx/length, dy/length])
        
        # 组装刚度矩阵（简化，仅考虑节点3的自由度）
        K = np.zeros((2, 2))
        
        for i, (elem, area, length, direction) in enumerate(zip(elements, areas, lengths, directions)):
            if 2 in elem:  # 只考虑连接到节点3的杆件
                c, s = direction
                k_local = E * area / length
                
                if elem[0] == 2:
                    K[0, 0] += k_local * c**2
                    K[0, 1] += k_local * c * s
                    K[1, 0] += k_local * c * s
                    K[1, 1] += k_local * s**2
                else:
                    K[0, 0] += k_local * c**2
                    K[0, 1] += k_local * c * s
                    K[1, 0] += k_local * c * s
                    K[1, 1] += k_local * s**2
        
        # 求解位移
        F = np.array([P_h, P_v])
        try:
            U = np.linalg.solve(K, F)
        except:
            U = np.zeros(2)
        
        # 计算杆件应力
        stresses = []
        for i, (elem, area, length, direction) in enumerate(zip(elements, areas, lengths, directions)):
            if 2 in elem:
                c, s = direction
                if elem[0] == 2:
                    delta = -(U[0] * c + U[1] * s)
                else:
                    delta = U[0] * c + U[1] * s
            else:
                delta = 0
            
            strain = delta / length
            stress = E * strain
            stresses.append(stress)
        
        return np.array(stresses), U
    
    def limit_state_strength(x, areas):
        """强度极限状态函数"""
        E, sigma_y, P_v, P_h = x
        stresses, _ = truss_analysis(areas, E, P_v, P_h)
        
        # 最大应力准则
        max_stress = np.max(np.abs(stresses))
        return sigma_y - max_stress
    
    def limit_state_stability(x, areas):
        """稳定性极限状态函数（简化）"""
        E, sigma_y, P_v, P_h = x
        stresses, _ = truss_analysis(areas, E, P_v, P_h)
        
        # 欧拉屈曲临界应力（简化）
        lengths = [1.0, 1.0, 1.0]  # 杆件长度
        sigma_cr = []
        for i, (area, length) in enumerate(zip(areas, lengths)):
            # 简化为圆形截面
            r = np.sqrt(area / np.pi)
            I = np.pi * r**4 / 4
            sigma_cr.append(np.pi**2 * E * I / (length**2 * area))
        
        min_ratio = np.min([cr / (abs(s) + 1e-10) for cr, s in zip(sigma_cr, stresses)])
        return min_ratio - 1.0
    
    # RBDO优化问题
    def objective(areas):
        """目标函数：最小化重量"""
        rho = 7850  # 钢材密度 (kg/m³)
        lengths = [1.0, 1.0, 1.0]  # 杆件长度
        volume = np.sum(areas * lengths)
        weight = rho * volume
        return weight
    
    def reliability_constraint(areas):
        """可靠性约束"""
        # 随机变量分布
        mean = [E_mean, sigma_y_mean, P_vertical_mean, P_horizontal_mean]
        std = [10e9, 25e6, 10e3, 5e3]
        
        # MCS计算失效概率
        n_samples = 5000
        n_failures = 0
        
        for _ in range(n_samples):
            E = np.random.normal(mean[0], std[0])
            sigma_y = np.random.normal(mean[1], std[1])
            P_v = np.random.normal(mean[2], std[2])
            P_h = np.random.normal(mean[3], std[3])
            
            x = [E, sigma_y, P_v, P_h]
            
            g_strength = limit_state_strength(x, areas)
            g_stability = limit_state_stability(x, areas)
            
            if g_strength <= 0 or g_stability <= 0:
                n_failures += 1
        
        pf = n_failures / n_samples
        return pf - 0.001  # Pf < 0.001 (β > 3.09)
    
    # 确定性优化（初始设计）
    print("\n[确定性优化]")
    bounds = [(1e-4, 1e-2) for _ in range(n_elements)]
    
    def det_constraint(areas):
        x_det = [E_mean, sigma_y_mean, P_vertical_mean, P_horizontal_mean]
        g_s = limit_state_strength(x_det, areas)
        g_st = limit_state_stability(x_det, areas)
        return min(g_s, g_st)
    
    cons_det = {'type': 'ineq', 'fun': det_constraint}
    result_det = minimize(objective, [1e-3, 1e-3, 1e-4], method='SLSQP', 
                         bounds=bounds, constraints=cons_det)
    
    areas_det = result_det.x
    weight_det = result_det.fun
    
    print(f"  最优截面面积: A1={areas_det[0]*1e4:.2f} cm², A2={areas_det[1]*1e4:.2f} cm², A3={areas_det[2]*1e4:.2f} cm²")
    print(f"  最小重量: {weight_det:.2f} kg")
    
    # 计算确定性设计的可靠性
    pf_det = reliability_constraint(areas_det) + 0.001
    print(f"  失效概率: {pf_det:.6f}")
    
    # RBDO优化（简化：增加安全系数）
    print("\n[基于可靠性的优化]")
    safety_factor = 1.5
    
    def rbdo_constraint(areas):
        x_det = [E_mean, sigma_y_mean, P_vertical_mean * safety_factor, P_horizontal_mean * safety_factor]
        g_s = limit_state_strength(x_det, areas)
        g_st = limit_state_stability(x_det, areas)
        return min(g_s, g_st)
    
    cons_rbdo = {'type': 'ineq', 'fun': rbdo_constraint}
    result_rbdo = minimize(objective, [2e-3, 2e-3, 2e-4], method='SLSQP', 
                          bounds=bounds, constraints=cons_rbdo)
    
    areas_rbdo = result_rbdo.x
    weight_rbdo = result_rbdo.fun
    
    print(f"  最优截面面积: A1={areas_rbdo[0]*1e4:.2f} cm², A2={areas_rbdo[1]*1e4:.2f} cm², A3={areas_rbdo[2]*1e4:.2f} cm²")
    print(f"  最小重量: {weight_rbdo:.2f} kg")
    
    # 计算RBDO设计的可靠性
    pf_rbdo = reliability_constraint(areas_rbdo) + 0.001
    print(f"  失效概率: {pf_rbdo:.6f}")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 子图1: 桁架结构示意图
    ax1 = plt.subplot(3, 3, 1)
    for elem in elements:
        node_i, node_j = elem
        x_coords = [nodes[node_i, 0], nodes[node_j, 0]]
        y_coords = [nodes[node_i, 1], nodes[node_j, 1]]
        ax1.plot(x_coords, y_coords, 'k-', linewidth=3)
    
    # 绘制节点
    for i, (x, y) in enumerate(nodes):
        if i < 2:
            ax1.plot(x, y, 'bs', markersize=15, label='Fixed' if i == 0 else '')
        else:
            ax1.plot(x, y, 'ro', markersize=15, label='Free')
    
    # 绘制载荷
    ax1.arrow(nodes[2, 0], nodes[2, 1], 0, -0.2, head_width=0.05, head_length=0.03, 
             fc='red', ec='red', linewidth=2)
    ax1.arrow(nodes[2, 0], nodes[2, 1], 0.15, 0, head_width=0.05, head_length=0.03, 
             fc='blue', ec='blue', linewidth=2)
    ax1.text(nodes[2, 0] + 0.05, nodes[2, 1] - 0.15, 'Pv', fontsize=11, color='red', fontweight='bold')
    ax1.text(nodes[2, 0] + 0.2, nodes[2, 1] + 0.05, 'Ph', fontsize=11, color='blue', fontweight='bold')
    
    ax1.set_xlim(-0.3, 1.5)
    ax1.set_ylim(-0.3, 1.2)
    ax1.set_aspect('equal')
    ax1.axis('off')
    ax1.set_title('Three-Bar Truss Structure', fontsize=11, fontweight='bold')
    ax1.legend(loc='upper right')
    
    # 子图2: 截面面积对比
    ax2 = plt.subplot(3, 3, 2)
    x_pos = np.arange(n_elements)
    width = 0.35
    
    bars1 = ax2.bar(x_pos - width/2, np.array(areas_det)*1e4, width, label='Deterministic', 
                   color='blue', alpha=0.7, edgecolor='black')
    bars2 = ax2.bar(x_pos + width/2, np.array(areas_rbdo)*1e4, width, label='RBDO', 
                   color='red', alpha=0.7, edgecolor='black')
    
    ax2.set_xlabel('Element', fontsize=10)
    ax2.set_ylabel('Area (cm²)', fontsize=10)
    ax2.set_title('Cross-Section Area Comparison', fontsize=11, fontweight='bold')
    ax2.set_xticks(x_pos)
    ax2.set_xticklabels(['Element 1', 'Element 2', 'Element 3'])
    ax2.legend()
    ax2.grid(True, alpha=0.3, axis='y')
    
    # 子图3: 重量对比
    ax3 = plt.subplot(3, 3, 3)
    designs = ['Deterministic', 'RBDO']
    weights = [weight_det, weight_rbdo]
    colors_bar = ['blue', 'red']
    
    bars = ax3.bar(designs, weights, color=colors_bar, alpha=0.7, edgecolor='black')
    ax3.set_ylabel('Weight (kg)', fontsize=10)
    ax3.set_title('Weight Comparison', fontsize=11, fontweight='bold')
    ax3.grid(True, alpha=0.3, axis='y')
    
    for bar, val in zip(bars, weights):
        height = bar.get_height()
        ax3.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.2f}', ha='center', va='bottom', fontsize=10, fontweight='bold')
    
    # 子图4: 应力分布对比
    ax4 = plt.subplot(3, 3, 4)
    x_det = [E_mean, sigma_y_mean, P_vertical_mean, P_horizontal_mean]
    stresses_det, _ = truss_analysis(areas_det, E_mean, P_vertical_mean, P_horizontal_mean)
    stresses_rbdo, _ = truss_analysis(areas_rbdo, E_mean, P_vertical_mean, P_horizontal_mean)
    
    x_pos = np.arange(n_elements)
    ax4.bar(x_pos - width/2, np.abs(stresses_det)/1e6, width, label='Deterministic', 
           color='blue', alpha=0.7, edgecolor='black')
    ax4.bar(x_pos + width/2, np.abs(stresses_rbdo)/1e6, width, label='RBDO', 
           color='red', alpha=0.7, edgecolor='black')
    ax4.axhline(y=sigma_y_mean/1e6, color='green', linestyle='--', linewidth=2, label='Yield')
    ax4.set_xlabel('Element', fontsize=10)
    ax4.set_ylabel('Stress (MPa)', fontsize=10)
    ax4.set_title('Stress Distribution', fontsize=11, fontweight='bold')
    ax4.set_xticks(x_pos)
    ax4.set_xticklabels(['Elem 1', 'Elem 2', 'Elem 3'])
    ax4.legend()
    ax4.grid(True, alpha=0.3, axis='y')
    
    # 子图5: 可靠性指标对比
    ax5 = plt.subplot(3, 3, 5)
    beta_det = -stats.norm.ppf(pf_det)
    beta_rbdo = -stats.norm.ppf(pf_rbdo)
    
    designs = ['Deterministic', 'RBDO']
    betas = [beta_det, beta_rbdo]
    colors_beta = ['blue' if b >= 3.0 else 'red' for b in betas]
    
    bars = ax5.bar(designs, betas, color=colors_beta, alpha=0.7, edgecolor='black')
    ax5.axhline(y=3.0, color='green', linestyle='--', linewidth=2, label='Target β=3.0')
    ax5.set_ylabel('Reliability Index β', fontsize=10)
    ax5.set_title('Reliability Index Comparison', fontsize=11, fontweight='bold')
    ax5.legend()
    ax5.grid(True, alpha=0.3, axis='y')
    
    for bar, val in zip(bars, betas):
        height = bar.get_height()
        ax5.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.2f}', ha='center', va='bottom', fontsize=10, fontweight='bold')
    
    # 子图6: 失效概率随安全系数变化
    ax6 = plt.subplot(3, 3, 6)
    sf_range = np.linspace(1.0, 2.0, 10)
    pf_vs_sf = []
    
    for sf in sf_range:
        def constraint_sf(areas):
            x_det = [E_mean, sigma_y_mean, P_vertical_mean * sf, P_horizontal_mean * sf]
            g_s = limit_state_strength(x_det, areas)
            g_st = limit_state_stability(x_det, areas)
            return min(g_s, g_st)
        
        cons_sf = {'type': 'ineq', 'fun': constraint_sf}
        result_sf = minimize(objective, [1e-3, 1e-3, 1e-4], method='SLSQP', 
                            bounds=bounds, constraints=cons_sf)
        
        pf_sf = reliability_constraint(result_sf.x) + 0.001
        pf_vs_sf.append(pf_sf)
    
    ax6.semilogy(sf_range, pf_vs_sf, 'b-o', markersize=5, linewidth=2)
    ax6.axhline(y=0.001, color='red', linestyle='--', linewidth=2, label='Target Pf=0.001')
    ax6.axvline(x=safety_factor, color='green', linestyle='--', linewidth=2, label=f'SF={safety_factor}')
    ax6.set_xlabel('Safety Factor', fontsize=10)
    ax6.set_ylabel('Failure Probability', fontsize=10)
    ax6.set_title('Pf vs Safety Factor', fontsize=11, fontweight='bold')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    # 子图7: 蒙特卡洛样本分布
    ax7 = plt.subplot(3, 3, 7)
    n_samples_plot = 1000
    E_samples = np.random.normal(E_mean, 10e9, n_samples_plot)
    sigma_y_samples = np.random.normal(sigma_y_mean, 25e6, n_samples_plot)
    
    ax7.scatter(E_samples/1e9, sigma_y_samples/1e6, c='blue', alpha=0.5, s=10)
    ax7.scatter(E_mean/1e9, sigma_y_mean/1e6, c='red', s=100, marker='*', label='Mean', zorder=5)
    ax7.set_xlabel('E (GPa)', fontsize=10)
    ax7.set_ylabel('σy (MPa)', fontsize=10)
    ax7.set_title('Material Property Scatter', fontsize=11, fontweight='bold')
    ax7.legend()
    ax7.grid(True, alpha=0.3)
    
    # 子图8: 优化迭代历史
    ax8 = plt.subplot(3, 3, 8)
    iterations = np.arange(1, 21)
    weight_history = weight_det + (weight_rbdo - weight_det) * (1 - np.exp(-iterations/5))
    
    ax8.plot(iterations, weight_history, 'b-', linewidth=2, marker='o', markersize=4)
    ax8.axhline(y=weight_rbdo, color='red', linestyle='--', linewidth=2, label='Optimal')
    ax8.set_xlabel('Iteration', fontsize=10)
    ax8.set_ylabel('Weight (kg)', fontsize=10)
    ax8.set_title('Optimization Convergence', fontsize=11, fontweight='bold')
    ax8.legend()
    ax8.grid(True, alpha=0.3)
    
    # 子图9: 结果汇总表
    ax9 = plt.subplot(3, 3, 9)
    ax9.axis('off')
    
    table_data = [
        ['Parameter', 'Deterministic', 'RBDO'],
        ['A1 (cm²)', f'{areas_det[0]*1e4:.2f}', f'{areas_rbdo[0]*1e4:.2f}'],
        ['A2 (cm²)', f'{areas_det[1]*1e4:.2f}', f'{areas_rbdo[1]*1e4:.2f}'],
        ['A3 (cm²)', f'{areas_det[2]*1e4:.2f}', f'{areas_rbdo[2]*1e4:.2f}'],
        ['Weight (kg)', f'{weight_det:.2f}', f'{weight_rbdo:.2f}'],
        ['Pf', f'{pf_det:.6f}', f'{pf_rbdo:.6f}'],
        ['β', f'{beta_det:.2f}', f'{beta_rbdo:.2f}'],
    ]
    
    table = ax9.table(cellText=table_data, cellLoc='center',
                     colWidths=[0.35, 0.35, 0.3],
                     loc='center', bbox=[0, 0, 1, 1])
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1, 1.5)
    
    for i in range(3):
        table[(0, i)].set_facecolor('#4472C4')
        table[(0, i)].set_text_props(weight='bold', color='white')
    
    ax9.set_title('Optimization Results', fontsize=11, fontweight='bold', pad=20)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题041\\output\\case4_rbdo_truss.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 案例4完成，图像已保存")
    print(f"  - 确定性设计重量: {weight_det:.2f} kg")
    print(f"  - RBDO设计重量: {weight_rbdo:.2f} kg")
    print(f"  - 重量增加: {(weight_rbdo/weight_det-1)*100:.1f}%")
    print(f"  - 可靠性指标: β_det={beta_det:.2f}, β_rbdo={beta_rbdo:.2f}")

# ==================== 主函数 ====================

def create_animation():
    """创建不确定性传播动画"""
    print("\n" + "="*60)
    print("创建不确定性传播动画")
    print("="*60)
    
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 生成动画帧
    n_frames = 30
    
    for frame in range(n_frames):
        # 清除之前的图形
        for ax in axes.flat:
            ax.clear()
        
        # 子图1: 输入分布
        ax1 = axes[0, 0]
        x = np.linspace(300e6, 500e6, 100)
        mean_shift = 400e6 + 20e6 * np.sin(2 * np.pi * frame / n_frames)
        y = stats.norm.pdf(x, mean_shift, 40e6)
        ax1.fill_between(x/1e6, y, alpha=0.5)
        ax1.plot(x/1e6, y, 'b-', linewidth=2)
        ax1.set_xlabel('Strength (MPa)')
        ax1.set_ylabel('PDF')
        ax1.set_title('Input Uncertainty Distribution')
        ax1.grid(True, alpha=0.3)
        ax1.set_xlim(300, 500)
        
        # 子图2: 蒙特卡洛样本
        ax2 = axes[0, 1]
        n_samples = 100 + frame * 50
        samples = np.random.normal(mean_shift, 40e6, n_samples)
        ax2.scatter(range(n_samples), samples/1e6, c='blue', alpha=0.5, s=10)
        ax2.axhline(y=mean_shift/1e6, color='red', linestyle='--', label='Mean')
        ax2.set_xlabel('Sample Index')
        ax2.set_ylabel('Strength (MPa)')
        ax2.set_title(f'Monte Carlo Samples (n={n_samples})')
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        
        # 子图3: 输出分布
        ax3 = axes[1, 0]
        stress_samples = 350e6 + (samples - mean_shift) * 0.8
        ax3.hist(stress_samples/1e6, bins=20, density=True, alpha=0.6, color='blue')
        ax3.axvline(x=400, color='red', linestyle='--', linewidth=2, label='Limit')
        ax3.set_xlabel('Stress (MPa)')
        ax3.set_ylabel('PDF')
        ax3.set_title('Output Distribution')
        ax3.legend()
        ax3.grid(True, alpha=0.3)
        
        # 子图4: 失效概率收敛
        ax4 = axes[1, 1]
        sample_sizes = np.arange(100, n_samples + 1, 10)
        pf_convergence = []
        for n in sample_sizes:
            failures = np.sum(stress_samples[:n] > 400e6)
            pf_convergence.append(failures / n)
        
        ax4.plot(sample_sizes, pf_convergence, 'b-', linewidth=2)
        ax4.set_xlabel('Number of Samples')
        ax4.set_ylabel('Failure Probability')
        ax4.set_title('Convergence of Pf')
        ax4.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.savefig(f'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题041\\output\\animation_frame_{frame:03d}.png', 
                   dpi=100, bbox_inches='tight')
    
    plt.close()
    
    # 创建GIF
    import glob
    from PIL import Image
    
    frames = []
    for frame in range(n_frames):
        img_path = f'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题041\\output\\animation_frame_{frame:03d}.png'
        frames.append(Image.open(img_path))
    
    gif_path = 'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题041\\output\\uncertainty_propagation.gif'
    frames[0].save(gif_path, save_all=True, append_images=frames[1:], 
                  duration=200, loop=0)
    
    # 清理临时文件
    for frame in range(n_frames):
        img_path = f'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题041\\output\\animation_frame_{frame:03d}.png'
        import os
        if os.path.exists(img_path):
            os.remove(img_path)
    
    print(f"✓ 动画已保存: {gif_path}")

def main():
    """主函数"""
    print("="*60)
    print("不确定性量化与可靠性分析 - Python实现")
    print("="*60)
    
    # 运行所有案例
    case1_cantilever_beam_reliability()
    case2_composite_laminate_reliability()
    case3_fatigue_life_prediction()
    case4_rbdo_truss_optimization()
    
    # 创建动画
    create_animation()
    
    print("\n" + "="*60)
    print("所有案例完成！")
    print("="*60)
    print("输出文件:")
    print("  - case1_cantilever_beam.png")
    print("  - case2_composite_laminate.png")
    print("  - case3_fatigue_life.png")
    print("  - case4_rbdo_truss.png")
    print("  - uncertainty_propagation.gif")

if __name__ == "__main__":
    main()