主题059：不确定性优化

1. 引言

在实际工程设计中，不确定性无处不在。材料属性的波动、制造公差、载荷变化、环境条件等因素都会影响结构性能。传统的确定性优化往往产生"脆弱"的最优解——在名义条件下性能优异，但在实际工况中可能失效。不确定性优化（Uncertainty-based Optimization）正是为解决这一问题而诞生的。

1.1 什么是不确定性优化？

不确定性优化是一类考虑设计变量、参数或环境存在不确定性时的优化方法。其核心思想是：

不再追求理论上的最优，而是追求在实际工况中稳健可靠的性能

1.2 为什么需要不确定性优化？

考虑一个经典案例：某桥梁设计通过确定性优化，在名义载荷下应力恰好达到许用值。但实际中：

• 材料强度可能有±10%的波动
• 载荷可能因极端天气增加30%
• 制造误差导致截面尺寸偏差

结果：这座"最优"设计的桥梁在实际使用中频繁出现裂缝，维护成本远超预期。

1.3 本教程学习目标

• 理解不确定性的数学表征方法
• 掌握蒙特卡洛模拟技术
• 学会鲁棒优化（Robust Optimization）
• 理解可靠性优化（Reliability-based Optimization）
• 能够处理实际工程中的不确定性问题

---

2. 不确定性的数学表征

2.1 不确定性类型

工程中的不确定性主要分为以下几类：

2.1.1 概率不确定性（Aleatory Uncertainty）

本质：固有的随机性，无法通过增加数据消除

例子：

• 材料强度的自然波动
• 风载荷的随机变化
• 制造过程中的随机误差

数学表征：概率密度函数（PDF）

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

2.1.2 认知不确定性（Epistemic Uncertainty）

本质：由于知识不足导致的不确定性，可通过增加数据减少

例子：

• 新材料的性能数据不足
• 极端载荷的历史记录有限
• 复杂物理现象的模型简化

数学表征：区间分析、模糊集、证据理论

2.1.3 混合不确定性

实际工程中往往是多种不确定性的组合，需要综合处理方法。

2.2 不确定性量化方法

2.2.1 蒙特卡洛模拟（Monte Carlo Simulation）

原理：通过大量随机采样估计统计特性

算法步骤：

1. 定义输入变量的概率分布
2. 从分布中随机采样 N 次
3. 对每个样本执行确定性分析
4. 统计输出结果的分布特性

收敛性分析：

蒙特卡洛的标准误差与样本数的关系：

$${ϵ}_{MC}\propto \frac{1}{\sqrt{N}}$$

这意味着：

• 样本增加4倍，误差减半
• 达到1%精度通常需要10,000+样本

代码实现：

def monte_carlo_simulation(self, design_vars, n_samples=1000):
    """蒙特卡洛模拟"""
    results = {'stress': [], 'deflection': [], 'weight': []}
    
    for _ in range(n_samples):
        # 采样不确定参数
        E = np.random.normal(E_nominal, E_std)
        P = np.random.normal(P_nominal, P_std)
        sigma_allow = np.random.uniform(sigma_low, sigma_high)
        
        # 确定性分析
        result = self.analyze(design_vars, 
                             {'E': E, 'P': P, 'sigma_allow': sigma_allow})
        
        for key in results:
            results[key].append(result[key])
    
    # 计算统计量
    stats = {}
    for key in results:
        data = np.array(results[key])
        stats[key] = {
            'mean': np.mean(data),
            'std': np.std(data),
            'min': np.min(data),
            'max': np.max(data),
            'p95': np.percentile(data, 95),
            'p5': np.percentile(data, 5)
        }
    
    return stats

2.2.2 可靠性指数（Reliability Index）

极限状态函数（Limit State Function）：

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

其中：

• $R$：抗力（Resistance）
• $S$：载荷效应（Load Effect）

失效概率：

$$P_f=P(g(X)\le 0)$$

可靠性指数（Hasofer-Lind）：

$$\beta =\underset{g(X)=0}{\min }\sqrt{(X-\mu {)}^T{\Sigma }^{-1}(X-\mu )}$$

$\beta$ 与失效概率的关系：

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

可靠性指数 $\beta$	失效概率 $P_f$	可靠性 $R$
1.0	15.87%	84.13%
1.65	4.95%	95.05%
2.33	0.99%	99.01%
3.0	0.13%	99.87%
3.72	0.01%	99.99%

---

3. 鲁棒优化（Robust Optimization）

3.1 鲁棒性定义

鲁棒性是指设计对不确定性的不敏感性。一个鲁棒的设计应该：

1. 性能稳定：性能指标的方差小
2. 约束满足：在不确定性下仍能满足约束
3. 容错能力强：对参数扰动不敏感

3.2 鲁棒性度量方法

3.2.1 均值-方差法（Mean-Variance）

目标函数：

$$\min f_{robust}={\mu }_f+k\cdot {\sigma }_f$$

其中：

• ${\mu }_f$：目标函数的均值
• ${\sigma }_f$：目标函数的标准差
• $k$：风险厌恶系数（$k>0$）

物理意义：

• $k=0$：传统确定性优化
• $k=1$：平衡均值和方差
• $k=3$：保守设计，强调稳定性

适用场景：

• 质量要求高的产品
• 成本波动敏感的项目
• 需要稳定性能的系统

3.2.2 最坏情况法（Worst-Case）

目标函数：

$$\min f_{robust}=\underset{u\in U}{\max }f(x,u)$$

其中 $U$ 是不确定参数的集合。

特点：

• 极度保守
• 保证在所有可能情况下都可行
• 适用于安全关键系统

适用场景：

• 航空航天
• 核工程
• 医疗器械

3.2.3 田口法（Taguchi Method）

信噪比（Signal-to-Noise Ratio）：

对于"越小越好"的目标：

$$SN=-10{\log }_{10}\left(\frac{1}{n}\sum _{i=1}^n{y}_i^2\right)$$

对于"越接近目标越好"：

$$SN=-10{\log }_{10}\left(\frac{1}{n}\sum _{i=1}^n(y_i-m{)}^2\right)$$

特点：

• 同时考虑均值和方差
• 源自实验设计理论
• 适合离散参数优化

3.2.4 六西格玛法（Six Sigma）

目标：确保性能在 $\pm 6\sigma$ 范围内都满足要求

目标函数：

$$\min f_{robust}={\mu }_f+3{\sigma }_f$$

约束：

$${\mu }_g-6{\sigma }_g\ge 0$$

特点：

• 极高的可靠性要求
• 制造业质量控制标准
• 缺陷率低于3.4ppm

3.3 鲁棒约束处理

传统约束 $g(x)\le 0$ 在不确定性下需要修改：

3.3.1 机会约束（Chance Constraints）

$$P(g(x,\xi )\le 0)\ge 1-\alpha$$

其中 $\alpha$ 是允许的最大失效概率。

转换方法：

假设 $g$ 服从正态分布：

$${\mu }_g+{\Phi }^{-1}(1-\alpha )\cdot {\sigma }_g\le 0$$

3.3.2 鲁棒约束的蒙特卡洛实现

def robust_constraints(self, design_vars):
    """鲁棒约束"""
    stats = self.monte_carlo_simulation(design_vars, n_samples=300)
    
    constraints = []
    
    # 1. 应力约束（99%可靠性）
    stress_mean = stats['stress']['mean']
    stress_std = stats['stress']['std']
    beta = (sigma_allow - stress_mean) / stress_std
    constraints.append(beta - 2.33)  # 要求 beta > 2.33
    
    # 2. 挠度约束（95%分位数）
    delta_p95 = stats['deflection']['p95']
    delta_limit = L / 250
    constraints.append(delta_limit - delta_p95)
    
    # 3. 安全因子约束（5%分位数）
    sf_p5 = stats['safety_factor']['p5']
    constraints.append(sf_p5 - 1.5)
    
    return constraints

3.4 悬臂梁鲁棒优化案例

问题描述：

设计一个悬臂梁，考虑以下不确定性：

• 弹性模量 $E\sim N(200000,10000^2)$ MPa
• 载荷 $P\sim N(1000,100^2)$ N
• 许用应力 ${\sigma }_{allow}\in [200,300]$ MPa（区间不确定性）

设计变量：

• 宽度 $b\in [10,50]$ mm
• 高度 $h\in [20,100]$ mm

目标：最小化重量

约束：

• 应力可靠性 > 99%
• 挠度 < L/250
• 安全因子 > 1.5

优化结果对比：

方法	宽度(mm)	高度(mm)	重量(kg)	应力(MPa)	安全因子(5%)
均值-方差法	17.25	52.19	7.07	127.3±12.8	1.52
最坏情况法	24.61	49.79	9.62	98.3±9.7	1.97
田口法	10.00	100.00	7.85	59.9±6.1	3.24

结果分析：

1. 均值-方差法：平衡了重量和稳定性，是最经济的选择
2. 最坏情况法：最保守，重量增加36%，但可靠性最高
3. 田口法：充分利用高度方向的材料效率，应力最低

---

4. 可靠性优化（Reliability-based Design Optimization, RBDO）

4.1 RBDO vs 鲁棒优化

特性	鲁棒优化	可靠性优化
关注点	性能稳定性	失效概率
目标	减小方差	满足可靠性要求
约束	概率约束	可靠性约束
计算量	较大	大
适用场景	性能波动敏感	安全关键系统

4.2 可靠性分析方法

4.2.1 FORM（First Order Reliability Method）

基本思想：在标准正态空间中寻找极限状态曲面上距离原点最近的点（MPP）。

算法步骤：

1. 将随机变量转换为标准正态变量 U
2. 迭代寻找 MPP（Most Probable Point）
3. 计算可靠性指数 β = ||u*||
4. 计算失效概率 Pf = Φ(-β)

数学推导：

变换到标准正态空间：

$$U_i=\frac{X_i-{\mu }_i}{{\sigma }_i}$$

极限状态函数：

$$g(U)=0$$

优化问题：

$$\begin{gathered} \min \beta =\sqrt{U^{T}U} \\ s.t.g(U)=0 \end{gathered}$$

4.2.2 SORM（Second Order Reliability Method）

在FORM基础上，在MPP处用二次曲面近似极限状态函数：

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

其中 ${\kappa }_i$ 是主曲率。

4.3 RBDO的数学表述

双层优化：

$$\begin{gathered} \underset{x}{\min }f(x) \\ s.t.{\beta }_i(x)\ge {\beta }_{target},i=1,...,m \end{gathered}$$

其中内层优化计算可靠性指数：

$$\begin{gathered} \beta (x)=\underset{u}{\min }||u|| \\ s.t.g(x,u)=0 \end{gathered}$$

单层等价形式（SORA方法）：

将可靠性约束转化为确定性约束：

$$g(x+\Delta x)\le 0$$

其中 $\Delta x$ 是随机向量的平移量。

4.4 悬臂梁可靠性优化案例

目标：最小化重量

约束：可靠性指数 $\beta \ge 2.33$（对应99%可靠性）

优化结果：

初始设计: 宽度=30mm, 高度=60mm
初始可靠性指数: 6.614 (目标: 2.326)

优化完成!
最优设计: 宽度=10.00mm, 高度=59.09mm
可靠性指数: 2.326
蒙特卡洛验证可靠性: 100.00%

结果分析：

• 初始设计过于保守（$\beta =6.6$）
• 优化后设计刚好满足可靠性要求
• 重量从约14kg降至4.6kg，减重67%
• 蒙特卡洛验证确认可靠性达标

---

5. 机会约束优化（Chance-Constrained Optimization）

5.1 机会约束的概念

传统约束要求100%满足，但在不确定性下可能过于保守。机会约束允许小概率违反：

$$P(g(x,\xi )\le 0)\ge 1-\alpha$$

其中 $\alpha$ 是风险水平（如5%）。

5.2 机会约束的确定性等价形式

5.2.1 线性机会约束

对于线性约束 $a^{T}x\le b$，其中 $a$ 是随机向量：

$$P(a^{T}x\le b)\ge 1-\alpha$$

假设 $a\sim N({\mu }_a,{\Sigma }_a)$，则等价于：

$${\mu }_a^{T}x+{\Phi }^{-1}(1-\alpha )\sqrt{x^T{\Sigma }_{a}x}\le b$$

5.2.2 凸性分析

• 当 $\alpha \le 0.5$ 时，${\Phi }^{-1}(1-\alpha )\ge 0$
• 约束是凸的，可用凸优化求解

5.3 场景近似法

当解析形式不可得时，用蒙特卡洛近似：

$$\frac{1}{N}\sum _{i=1}^N\mathbb{I}(g(x,{\xi }_i)\le 0)\ge 1-\alpha$$

其中 $\mathbb{I}$ 是指示函数。

大M法转换：

引入二元变量 $z_i$：

g ( x , ξ i ) ≤ M z i 1 N ∑ i = 1 N ( 1 − z i ) ≥ 1 − α z i ∈ { 0 , 1 } g(x, \xi_i) \leq M z_i \\ \frac{1}{N} \sum_{i=1}^N (1 - z_i) \geq 1 - \alpha \\ z_i \in \{0, 1\} g(x,ξi​)≤Mzi​N1​i=1∑N​(1−zi​)≥1−αzi​∈{0,1}

---

6. 区间优化（Interval Optimization）

6.1 区间分析基础

区间数：$[a]=[\underline{a},\overline{a}]$

区间运算：

$$\begin{gathered} [a]+[b]=[\underline{a}+\underline{b},\overline{a}+\overline{b}] \\ [a]\cdot [b]=[\min (\underline{a}\underline{b},\underline{a}\overline{b},\overline{a}\underline{b},\overline{a}\overline{b}),\max (\underline{a}\underline{b},\underline{a}\overline{b},\overline{a}\underline{b},\overline{a}\overline{b})] \end{gathered}$$

6.2 区间优化方法

6.2.1 最坏情况优化

$$\underset{x}{\min }\underset{u\in [\underline{u},\overline{u}]}{\max }f(x,u)$$

6.2.2 区间可能度法

定义区间比较的可能度：

$$P([a]\le [b])=\frac{\max (0,\overline{b}-\overline{a})-\max (0,\underline{b}-\underline{a})}{\overline{a}-\underline{a}}$$

约束转化为：

$$P([g](x)\le 0)\ge \lambda$$

其中 $\lambda$ 是决策者给定的置信水平。

---

7. 实际应用案例

7.1 汽车碰撞安全性设计

不确定性来源：

• 碰撞速度：实际事故速度分布
• 材料性能：批次差异
• 制造工艺：焊接质量波动

优化策略：

• 乘员伤害指标（HIC）的95%分位数最小化
• 确保99%的碰撞场景下安全气囊正常展开

7.2 风力发电机叶片设计

不确定性来源：

• 风速分布：Weibull分布
• 材料疲劳性能
• 冰载荷

优化策略：

• 20年疲劳寿命的可靠性 > 95%
• 极端风载下的失效概率 < 0.1%

7.3 微电子封装热设计

不确定性来源：

• 芯片功耗波动
• 散热材料接触热阻
• 环境温度变化

优化策略：

• 结温的99%分位数 < 125°C
• 热循环疲劳寿命的变异系数 < 10%

---

8. 计算效率与近似方法

8.1 蒙特卡洛的替代方法

8.1.1 多项式混沌展开（PCE）

用正交多项式近似随机输出：

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

优点：

• 计算量小（通常几十次仿真）
• 可解析计算统计矩
• 适合平滑响应

8.1.2 克里金/高斯过程

用代理模型替代昂贵仿真：

$$f(x)\sim GP(\mu (x),k(x,x^{\prime }))$$

结合蒙特卡洛进行不确定性传播。

8.2 降维技术

8.2.1 主成分分析（PCA）

识别最重要的不确定性方向：

$$\xi \approx \mu +\sum _{i=1}^k\sqrt{{\lambda }_i}{\varphi }_i{Z}_i$$

8.2.2 主动子空间法

寻找对目标函数影响最大的方向：

$$C=\mathbb{E}[\nabla f\nabla f^T]=W\Lambda W^T$$

8.3 自适应采样

根据不确定性重要性自适应增加样本：

def adaptive_sampling():
    """自适应采样"""
    # 初始均匀采样
    samples = latin_hypercube_sampling(n_initial)
    
    for iteration in range(max_iter):
        # 计算当前统计量
        stats = compute_statistics(samples)
        
        # 识别高方差区域
        high_variance_regions = identify_regions(samples, stats)
        
        # 在高方差区域增加样本
        new_samples = sample_in_regions(high_variance_regions)
        samples = np.vstack([samples, new_samples])
        
        # 检查收敛
        if converged(stats):
            break
    
    return stats

---

9. 软件工具与实践建议

9.1 常用软件

软件	特点	适用场景
UQLab	MATLAB平台，PCE、Kriging	学术研究
Dakota	开源，多种UQ方法	工程应用
OpenCossan	可靠性分析专用	土木工程
ANSYS Prob	与CAE集成	工业仿真
Python + SciPy	灵活，开源	快速原型

9.2 实践建议

9.2.1 不确定性建模

1. 数据收集：尽可能收集历史数据
2. 分布选择：根据物理意义选择分布
   • 强度：Weibull或Lognormal
   • 尺寸误差：Normal
   • 计数：Poisson
3. 相关性处理：考虑参数间的相关性

9.2.2 优化策略选择

场景	推荐方法
快速评估	蒙特卡洛 + 代理模型
高精度要求	FORM/SORM
多失效模式	系统可靠性方法
数据稀缺	区间分析 + 贝叶斯更新

9.2.3 验证与确认

1. 交叉验证：用独立数据集验证模型
2. 敏感性分析：识别关键不确定性来源
3. 反向验证：用实测数据校准模型

---

10. 前沿研究方向

10.1 深度学习与UQ

• 深度高斯过程：处理高维非线性问题
• 神经网络不确定性：贝叶斯神经网络
• 生成模型：用VAE/GAN生成样本

10.2 多保真度方法

结合高低保真度模型：

$$f_{high}(x)\approx \rho f_{low}(x)+\delta (x)$$

大幅减少昂贵仿真次数。

10.3 实时不确定性量化

• 数字孪生中的在线UQ
• 流数据驱动的自适应更新
• 边缘计算实现

---

11. 总结

11.1 关键知识点回顾

1. 不确定性类型：概率型、区间型、混合型
2. 量化方法：蒙特卡洛、FORM/SORM、PCE
3. 鲁棒优化：均值-方差、最坏情况、田口法
4. 可靠性优化：RBDO、机会约束
5. 计算策略：代理模型、自适应采样、降维

11.2 方法选择指南

是否需要考虑不确定性？
├── 否 → 确定性优化
└── 是 → 性能是否关键？
    ├── 是 → 鲁棒优化
    └── 否 → 安全性是否关键？
        ├── 是 → 可靠性优化
        └── 否 → 机会约束优化

11.3 学习建议

1. 理论结合实践：从简单案例开始
2. 重视可视化：理解不确定性传播
3. 验证结果：用多种方法交叉验证
4. 关注前沿：跟踪深度学习方法

---

12. 思考题与练习

12.1 思考题

1. 为什么确定性优化的最优解往往对不确定性敏感？
2. 均值-方差法中，风险厌恶系数 $k$ 如何影响优化结果？
3. FORM方法的局限性是什么？什么情况下需要使用SORM？
4. 如何权衡鲁棒性和经济性？

12.2 编程练习

1. 修改代码，实现三杆桁架的鲁棒优化
2. 添加参数相关性处理（使用Copula）
3. 实现自适应蒙特卡洛方法
4. 对比不同代理模型在UQ中的精度

12.3 案例分析

选择一个你熟悉的工程问题：

1. 识别所有不确定性来源
2. 建立不确定性模型
3. 进行鲁棒或可靠性优化
4. 分析结果并给出设计建议

---

附录：代码清单

完整代码见：uncertainty_optimization.py

核心类说明

类名	功能
CantileverBeam	悬臂梁分析模型
RobustOptimizer	鲁棒优化器
ReliabilityBasedOptimizer	可靠性优化器
UncertainParameter	不确定参数定义

运行结果

运行代码将生成：

1. uncertainty_optimization.png - 鲁棒优化结果对比
2. mc_convergence.png - 蒙特卡洛收敛性分析

---

"""
主题059：不确定性优化
Uncertainty-based Optimization / Robust Design Optimization

本代码实现考虑不确定性的优化方法：
1. 鲁棒优化 (Robust Optimization)
2. 可靠性优化 (Reliability-based Optimization)
3. 机会约束优化 (Chance-constrained Optimization)
4. 区间优化 (Interval Optimization)

以悬臂梁设计为例，考虑材料属性、载荷的不确定性
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, differential_evolution
from scipy.stats import norm, uniform
from typing import Callable, Tuple, List, Dict, Optional
from dataclasses import dataclass
from enum import Enum
import warnings

warnings.filterwarnings('ignore')

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


class UncertaintyType(Enum):
    """不确定性类型"""
    PROBABILISTIC = "概率不确定性"
    INTERVAL = "区间不确定性"
    FUZZY = "模糊不确定性"
    MIXED = "混合不确定性"


class RobustnessMeasure(Enum):
    """鲁棒性度量方式"""
    MEAN_VARIANCE = "均值-方差法"
    WORST_CASE = "最坏情况法"
    TAGUCHI = "田口法"
    SIX_SIGMA = "六西格玛法"


@dataclass
class UncertainParameter:
    """
    不确定参数定义
    
    Attributes:
        name: 参数名称
        nominal: 名义值
        uncertainty_type: 不确定性类型
        distribution: 概率分布（如果是概率型）
        lower_bound: 下界（区间型）
        upper_bound: 上界（区间型）
        std: 标准差（概率型）
    """
    name: str
    nominal: float
    uncertainty_type: UncertaintyType
    distribution: Optional[str] = None
    lower_bound: Optional[float] = None
    upper_bound: Optional[float] = None
    std: Optional[float] = None


class CantileverBeam:
    """
    悬臂梁结构分析
    
    考虑以下不确定性：
    - 材料弹性模量 E
    - 载荷 P
    - 许用应力 sigma_allow
    """
    
    def __init__(self):
        # 确定性参数
        self.L = 1000  # 梁长度 (mm)
        
        # 不确定参数定义
        self.uncertain_params = {
            'E': UncertainParameter(
                name='弹性模量',
                nominal=200000,  # MPa
                uncertainty_type=UncertaintyType.PROBABILISTIC,
                distribution='normal',
                std=10000  # 5%变异系数
            ),
            'P': UncertainParameter(
                name='载荷',
                nominal=1000,  # N
                uncertainty_type=UncertaintyType.PROBABILISTIC,
                distribution='normal',
                std=100  # 10%变异系数
            ),
            'sigma_allow': UncertainParameter(
                name='许用应力',
                nominal=250,  # MPa
                uncertainty_type=UncertaintyType.INTERVAL,
                lower_bound=200,
                upper_bound=300
            )
        }
    
    def analyze(self, design_vars: np.ndarray, 
                uncertain_values: Dict[str, float]) -> Dict:
        """
        悬臂梁分析
        
        设计变量：
        - x[0]: 宽度 b (mm)
        - x[1]: 高度 h (mm)
        
        返回：应力、挠度、重量
        """
        b, h = design_vars
        E = uncertain_values.get('E', self.uncertain_params['E'].nominal)
        P = uncertain_values.get('P', self.uncertain_params['P'].nominal)
        
        # 截面惯性矩
        I = b * h**3 / 12
        
        # 最大应力（固定端）
        M_max = P * self.L  # 最大弯矩
        sigma_max = M_max * (h/2) / I  # 弯曲应力
        
        # 最大挠度（自由端）
        delta_max = P * self.L**3 / (3 * E * I)
        
        # 重量（简化）
        density = 7.85e-6  # kg/mm^3 (钢材)
        weight = density * b * h * self.L
        
        return {
            'stress': sigma_max,
            'deflection': delta_max,
            'weight': weight,
            'safety_factor': uncertain_values.get('sigma_allow', 250) / sigma_max
        }
    
    def monte_carlo_simulation(self, design_vars: np.ndarray, 
                               n_samples: int = 1000) -> Dict:
        """
        蒙特卡洛模拟
        
        评估设计在不确定性下的性能分布
        """
        results = {
            'stress': [],
            'deflection': [],
            'weight': [],
            'safety_factor': []
        }
        
        for _ in range(n_samples):
            # 采样不确定参数
            uncertain_values = {}
            
            # E: 正态分布
            E_param = self.uncertain_params['E']
            uncertain_values['E'] = np.random.normal(
                E_param.nominal, E_param.std
            )
            
            # P: 正态分布
            P_param = self.uncertain_params['P']
            uncertain_values['P'] = np.random.normal(
                P_param.nominal, P_param.std
            )
            
            # sigma_allow: 区间均匀分布
            sigma_param = self.uncertain_params['sigma_allow']
            uncertain_values['sigma_allow'] = np.random.uniform(
                sigma_param.lower_bound, sigma_param.upper_bound
            )
            
            # 分析
            result = self.analyze(design_vars, uncertain_values)
            
            for key in results:
                results[key].append(result[key])
        
        # 计算统计量
        stats = {}
        for key in results:
            data = np.array(results[key])
            stats[key] = {
                'mean': np.mean(data),
                'std': np.std(data),
                'min': np.min(data),
                'max': np.max(data),
                'p95': np.percentile(data, 95),
                'p5': np.percentile(data, 5)
            }
        
        return stats


class RobustOptimizer:
    """
    鲁棒优化器
    
    实现多种鲁棒优化方法
    """
    
    def __init__(self, beam: CantileverBeam, 
                 robustness_measure: RobustnessMeasure = RobustnessMeasure.MEAN_VARIANCE):
        self.beam = beam
        self.robustness_measure = robustness_measure
        self.bounds = [(10, 50), (20, 100)]  # 宽度、高度边界
        
    def robust_objective(self, design_vars: np.ndarray) -> float:
        """
        鲁棒目标函数
        
        最小化：均值 + 惩罚项
        """
        # 蒙特卡洛评估
        stats = self.beam.monte_carlo_simulation(design_vars, n_samples=500)
        
        if self.robustness_measure == RobustnessMeasure.MEAN_VARIANCE:
            # 均值-方差法
            mean_weight = stats['weight']['mean']
            std_weight = stats['weight']['std']
            
            # 鲁棒目标：均值 + k * 标准差
            k = 1.0  # 风险厌恶系数
            robust_obj = mean_weight + k * std_weight
            
        elif self.robustness_measure == RobustnessMeasure.WORST_CASE:
            # 最坏情况法
            robust_obj = stats['weight']['max']
            
        elif self.robustness_measure == RobustnessMeasure.TAGUCHI:
            # 田口法：信噪比
            mean_weight = stats['weight']['mean']
            std_weight = stats['weight']['std']
            
            # 信噪比（越大越好，所以取负）
            sn_ratio = -10 * np.log10(mean_weight**2 + std_weight**2)
            robust_obj = -sn_ratio
            
        else:  # SIX_SIGMA
            # 六西格玛法
            mean_weight = stats['weight']['mean']
            std_weight = stats['weight']['std']
            
            # 考虑6倍标准差的范围
            robust_obj = mean_weight + 3 * std_weight
        
        return robust_obj
    
    def robust_constraints(self, design_vars: np.ndarray) -> List[float]:
        """
        鲁棒约束
        
        约束需要以高概率满足
        """
        # 蒙特卡洛评估
        stats = self.beam.monte_carlo_simulation(design_vars, n_samples=300)
        
        constraints = []
        
        # 1. 应力约束（可靠性要求：P(stress < sigma_allow) > 0.99）
        stress_mean = stats['stress']['mean']
        stress_std = stats['stress']['std']
        
        # 计算可靠性指数 beta
        sigma_allow_mean = (self.beam.uncertain_params['sigma_allow'].lower_bound + 
                           self.beam.uncertain_params['sigma_allow'].upper_bound) / 2
        beta = (sigma_allow_mean - stress_mean) / stress_std
        
        # 要求 beta > 2.33 (对应99%可靠性)
        constraints.append(beta - 2.33)
        
        # 2. 挠度约束（95%分位数）
        delta_p95 = stats['deflection']['p95']
        delta_limit = self.beam.L / 250  # L/250
        constraints.append(delta_limit - delta_p95)
        
        # 3. 安全因子约束
        sf_p5 = stats['safety_factor']['p5']  # 5%分位数
        constraints.append(sf_p5 - 1.5)
        
        return constraints
    
    def optimize(self) -> Tuple[np.ndarray, Dict]:
        """执行鲁棒优化"""
        print("=" * 70)
        print(f"鲁棒优化 ({self.robustness_measure.value})")
        print("=" * 70)
        
        # 初始设计
        x0 = np.array([25, 50])
        
        print(f"初始设计: 宽度={x0[0]}mm, 高度={x0[1]}mm")
        
        # 评估初始设计
        init_stats = self.beam.monte_carlo_simulation(x0, n_samples=1000)
        print(f"\n初始设计性能:")
        print(f"  重量: {init_stats['weight']['mean']:.3f} ± {init_stats['weight']['std']:.3f} kg")
        print(f"  应力: {init_stats['stress']['mean']:.2f} ± {init_stats['stress']['std']:.2f} MPa")
        print(f"  挠度: {init_stats['deflection']['mean']:.2f} ± {init_stats['deflection']['std']:.2f} mm")
        print(f"  安全因子(5%): {init_stats['safety_factor']['p5']:.2f}")
        
        # 定义约束函数
        def constraint_func(x):
            return self.robust_constraints(x)
        
        # 包装约束
        constraints = []
        for i in range(3):
            constraints.append({
                'type': 'ineq',
                'fun': lambda x, idx=i: constraint_func(x)[idx]
            })
        
        # 执行优化
        print(f"\n开始优化...")
        result = minimize(
            self.robust_objective,
            x0,
            method='SLSQP',
            bounds=self.bounds,
            constraints=constraints,
            options={'maxiter': 50, 'ftol': 1e-4}
        )
        
        x_opt = result.x
        
        print(f"\n优化完成!")
        print(f"最优设计: 宽度={x_opt[0]:.2f}mm, 高度={x_opt[1]:.2f}mm")
        
        # 评估最优设计
        opt_stats = self.beam.monte_carlo_simulation(x_opt, n_samples=1000)
        print(f"\n最优设计性能:")
        print(f"  重量: {opt_stats['weight']['mean']:.3f} ± {opt_stats['weight']['std']:.3f} kg")
        print(f"  应力: {opt_stats['stress']['mean']:.2f} ± {opt_stats['stress']['std']:.2f} MPa")
        print(f"  挠度: {opt_stats['deflection']['mean']:.2f} ± {opt_stats['deflection']['std']:.2f} mm")
        print(f"  安全因子(5%): {opt_stats['safety_factor']['p5']:.2f}")
        
        return x_opt, {
            'x_opt': x_opt,
            'initial_stats': init_stats,
            'optimal_stats': opt_stats,
            'robustness_measure': self.robustness_measure
        }


class ReliabilityBasedOptimizer:
    """
    可靠性优化器
    
    基于FORM/SORM方法的可靠性优化
    """
    
    def __init__(self, beam: CantileverBeam, target_reliability: float = 0.99):
        self.beam = beam
        self.target_reliability = target_reliability
        self.beta_target = norm.ppf(target_reliability)  # 目标可靠性指数
        self.bounds = [(10, 50), (20, 100)]
        
    def form_analysis(self, design_vars: np.ndarray) -> float:
        """
        一次可靠性方法（FORM）
        
        计算可靠性指数 beta
        """
        b, h = design_vars
        
        # 名义值分析
        nominal_values = {
            'E': self.beam.uncertain_params['E'].nominal,
            'P': self.beam.uncertain_params['P'].nominal,
            'sigma_allow': self.beam.uncertain_params['sigma_allow'].nominal
        }
        
        result = self.beam.analyze(design_vars, nominal_values)
        
        # 极限状态函数：g = sigma_allow - sigma_max
        # 简化计算：假设线性关系
        
        # 灵敏度（解析导数）
        # sigma = 6*P*L / (b*h^2)
        P = nominal_values['P']
        L = self.beam.L
        sigma_nominal = 6 * P * L / (b * h**2)
        
        # 对P的灵敏度
        dsigma_dP = 6 * L / (b * h**2)
        
        # 对sigma_allow的灵敏度
        dsigma_dsigma_allow = 1.0
        
        # 标准差
        sigma_P = self.beam.uncertain_params['P'].std
        sigma_sigma_allow = (self.beam.uncertain_params['sigma_allow'].upper_bound - 
                            self.beam.uncertain_params['sigma_allow'].lower_bound) / (2 * np.sqrt(3))
        
        # 极限状态函数的标准差
        sigma_g = np.sqrt((dsigma_dP * sigma_P)**2 + 
                         (dsigma_dsigma_allow * sigma_sigma_allow)**2)
        
        # 可靠性指数
        g_mean = nominal_values['sigma_allow'] - sigma_nominal
        beta = g_mean / sigma_g if sigma_g > 0 else float('inf')
        
        return beta
    
    def optimize(self) -> Tuple[np.ndarray, Dict]:
        """执行可靠性优化"""
        print("\n" + "=" * 70)
        print(f"可靠性优化 (目标可靠性: {self.target_reliability*100:.1f}%)")
        print("=" * 70)
        
        def objective(x):
            """最小化重量"""
            stats = self.beam.monte_carlo_simulation(x, n_samples=100)
            return stats['weight']['mean']
        
        def reliability_constraint(x):
            """可靠性约束"""
            beta = self.form_analysis(x)
            return beta - self.beta_target
        
        # 初始设计
        x0 = np.array([30, 60])
        
        print(f"初始设计: 宽度={x0[0]}mm, 高度={x0[1]}mm")
        init_beta = self.form_analysis(x0)
        print(f"初始可靠性指数: {init_beta:.3f} (目标: {self.beta_target:.3f})")
        
        # 优化
        constraints = [{'type': 'ineq', 'fun': reliability_constraint}]
        
        result = minimize(
            objective,
            x0,
            method='SLSQP',
            bounds=self.bounds,
            constraints=constraints,
            options={'maxiter': 50}
        )
        
        x_opt = result.x
        opt_beta = self.form_analysis(x_opt)
        
        print(f"\n优化完成!")
        print(f"最优设计: 宽度={x_opt[0]:.2f}mm, 高度={x_opt[1]:.2f}mm")
        print(f"可靠性指数: {opt_beta:.3f}")
        
        # 蒙特卡洛验证
        opt_stats = self.beam.monte_carlo_simulation(x_opt, n_samples=5000)
        failure_count = np.sum(np.array(opt_stats['stress']['mean']) > 
                              self.beam.uncertain_params['sigma_allow'].nominal)
        actual_reliability = 1 - failure_count / 5000
        
        print(f"蒙特卡洛验证可靠性: {actual_reliability*100:.2f}%")
        
        return x_opt, {
            'x_opt': x_opt,
            'beta': opt_beta,
            'actual_reliability': actual_reliability
        }


def compare_robustness_measures():
    """比较不同鲁棒性度量方法"""
    print("=" * 70)
    print("鲁棒性度量方法对比")
    print("=" * 70)
    
    beam = CantileverBeam()
    
    results = {}
    
    for measure in [RobustnessMeasure.MEAN_VARIANCE, 
                    RobustnessMeasure.WORST_CASE,
                    RobustnessMeasure.TAGUCHI]:
        print(f"\n{'='*50}")
        print(f"方法: {measure.value}")
        print('='*50)
        
        optimizer = RobustOptimizer(beam, robustness_measure=measure)
        x_opt, info = optimizer.optimize()
        
        results[measure.name] = {
            'x_opt': x_opt,
            'stats': info['optimal_stats']
        }
    
    # 对比总结
    print("\n" + "=" * 70)
    print("对比总结")
    print("=" * 70)
    
    for name, result in results.items():
        stats = result['stats']
        print(f"\n{name}:")
        print(f"  设计: b={result['x_opt'][0]:.2f}mm, h={result['x_opt'][1]:.2f}mm")
        print(f"  重量: {stats['weight']['mean']:.3f} ± {stats['weight']['std']:.3f} kg")
        print(f"  应力: {stats['stress']['mean']:.2f} ± {stats['stress']['std']:.2f} MPa")
    
    # 可视化
    visualize_uncertainty_results(results)
    
    return results


def visualize_uncertainty_results(results: Dict):
    """可视化不确定性优化结果"""
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 重量对比
    ax = axes[0, 0]
    methods = list(results.keys())
    weights_mean = [results[m]['stats']['weight']['mean'] for m in methods]
    weights_std = [results[m]['stats']['weight']['std'] for m in methods]
    
    x_pos = np.arange(len(methods))
    bars = ax.bar(x_pos, weights_mean, yerr=weights_std, 
                  capsize=5, alpha=0.8, color='steelblue')
    ax.set_xticks(x_pos)
    ax.set_xticklabels(methods, rotation=15, ha='right')
    ax.set_ylabel('重量 (kg)')
    ax.set_title('重量对比（带标准差）')
    ax.grid(True, axis='y')
    
    # 2. 应力分布
    ax = axes[0, 1]
    
    stress_means = [results[m]['stats']['stress']['mean'] for m in methods]
    stress_stds = [results[m]['stats']['stress']['std'] for m in methods]
    
    ax.errorbar(methods, stress_means, yerr=stress_stds, 
               fmt='o-', capsize=5, linewidth=2, markersize=8)
    ax.axhline(y=250, color='r', linestyle='--', label='许用应力')
    ax.set_ylabel('应力 (MPa)')
    ax.set_title('应力对比')
    ax.legend()
    ax.grid(True)
    plt.setp(ax.xaxis.get_majorticklabels(), rotation=15, ha='right')
    
    # 3. 设计空间分布
    ax = axes[1, 0]
    
    for name, result in results.items():
        x_opt = result['x_opt']
        ax.scatter(x_opt[0], x_opt[1], s=200, label=name, alpha=0.7)
    
    ax.set_xlabel('宽度 b (mm)')
    ax.set_ylabel('高度 h (mm)')
    ax.set_title('最优设计分布')
    ax.legend()
    ax.grid(True)
    
    # 4. 安全因子分布
    ax = axes[1, 1]
    
    sf_data = []
    for name in methods:
        stats = results[name]['stats']
        sf_data.append([
            stats['safety_factor']['p5'],
            stats['safety_factor']['mean'],
            stats['safety_factor']['p95']
        ])
    
    sf_data = np.array(sf_data)
    
    x_pos = np.arange(len(methods))
    ax.bar(x_pos, sf_data[:, 1], alpha=0.7, color='green', label='均值')
    ax.errorbar(x_pos, sf_data[:, 1], 
               yerr=[sf_data[:, 1] - sf_data[:, 0], 
                     sf_data[:, 2] - sf_data[:, 1]],
               fmt='none', ecolor='black', capsize=5)
    
    ax.axhline(y=1.5, color='r', linestyle='--', label='最小要求')
    ax.set_xticks(x_pos)
    ax.set_xticklabels(methods, rotation=15, ha='right')
    ax.set_ylabel('安全因子')
    ax.set_title('安全因子分布')
    ax.legend()
    ax.grid(True, axis='y')
    
    plt.tight_layout()
    plt.savefig('uncertainty_optimization.png', dpi=150, bbox_inches='tight')
    print("\n图像已保存: uncertainty_optimization.png")
    plt.show()


def demonstrate_monte_carlo_convergence():
    """演示蒙特卡洛收敛性"""
    print("\n" + "=" * 70)
    print("蒙特卡洛收敛性分析")
    print("=" * 70)
    
    beam = CantileverBeam()
    design_vars = np.array([25, 50])
    
    sample_sizes = [10, 50, 100, 500, 1000, 5000]
    means = []
    stds = []
    
    for n in sample_sizes:
        stats = beam.monte_carlo_simulation(design_vars, n_samples=n)
        means.append(stats['stress']['mean'])
        stds.append(stats['stress']['std'])
        print(f"n={n:5d}: 应力均值={stats['stress']['mean']:.2f}, "
              f"标准差={stats['stress']['std']:.2f}")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    
    ax = axes[0]
    ax.plot(sample_sizes, means, 'o-', linewidth=2)
    ax.set_xlabel('样本数')
    ax.set_ylabel('应力均值 (MPa)')
    ax.set_title('蒙特卡洛收敛性 - 均值')
    ax.set_xscale('log')
    ax.grid(True)
    
    ax = axes[1]
    ax.plot(sample_sizes, stds, 'o-', linewidth=2, color='orange')
    ax.set_xlabel('样本数')
    ax.set_ylabel('应力标准差 (MPa)')
    ax.set_title('蒙特卡洛收敛性 - 标准差')
    ax.set_xscale('log')
    ax.grid(True)
    
    plt.tight_layout()
    plt.savefig('mc_convergence.png', dpi=150, bbox_inches='tight')
    print("\n图像已保存: mc_convergence.png")
    plt.show()


if __name__ == "__main__":
    # 主程序
    results = compare_robustness_measures()
    
    # 可靠性优化
    beam = CantileverBeam()
    reliability_optimizer = ReliabilityBasedOptimizer(beam, target_reliability=0.99)
    x_opt_rel, info_rel = reliability_optimizer.optimize()
    
    # 蒙特卡洛收敛性
    demonstrate_monte_carlo_convergence()