第067篇：可靠性优化与鲁棒设计

摘要

在实际工程设计中，不确定性无处不在：材料属性的波动、制造公差、载荷变化、环境条件等因素都会影响产品的性能和安全性。传统的确定性优化方法忽略了这些不确定性，可能导致设计在实际运行中失效或性能严重偏离预期。可靠性优化与鲁棒设计是处理不确定性下优化问题的两种重要方法。可靠性优化关注满足特定可靠性约束，确保设计在不确定性下的安全裕度；鲁棒设计则致力于降低性能对不确定性的敏感性，使产品性能更加稳定一致。本主题系统阐述可靠性优化与鲁棒设计的理论基础、数学模型和求解方法，包括可靠性指标、概率约束、稳健性度量、双循环法、单循环法、解耦法等。通过工程实例展示如何在多场耦合仿真中应用这些方法，实现既安全又稳健的设计方案。

关键词

可靠性优化, 鲁棒设计, 不确定性, 概率约束, 可靠性指标, 稳健性度量, 双循环法, 单循环法

---

1. 引言

1.1 工程中的不确定性问题

工程设计和制造过程中存在各种不确定性来源，这些不确定性可分为以下几类：

物理不确定性（偶然不确定性）：源于系统固有的随机性，如材料属性的自然波动、载荷的随机变化等。这类不确定性不可消除，只能用概率方法描述。

认知不确定性（认识不确定性）：源于知识不足或信息不完全，如模型简化带来的误差、参数识别的不准确等。这类不确定性可通过增加数据或改进模型来降低。

制造不确定性：包括尺寸公差、表面粗糙度、装配误差等制造过程中的变异。

环境不确定性：使用环境中的温度、湿度、载荷等条件的变化。

模型不确定性：数学模型对真实物理过程的近似程度带来的误差。

1.2 确定性优化的局限性

传统确定性优化问题可表述为：

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

$$\text{s.t.}{g}_j(\mathbf{d})\le 0,j=1,2,\dots ,m$$

其中，$\mathbf{d}$为设计变量，$f$为目标函数，$g_j$为约束函数。

这种表述假设所有参数都是确定性的。然而，当考虑不确定性时，确定性最优解可能：

1. 可靠性不足：约束边界附近的设计在参数波动时容易违反约束
2. 稳健性差：目标函数或约束函数对参数变化过于敏感
3. 实际性能偏离：名义最优在实际运行中可能性能很差

1.3 可靠性优化与鲁棒设计的概念

可靠性优化（Reliability-Based Design Optimization, RBDO）：在优化过程中显式考虑不确定性，确保设计满足预定的可靠性要求。其核心是将确定性约束转化为概率约束：

$$P(g_j(\mathbf{d},\mathbf{X})\le 0)\ge R_j$$

其中，$\mathbf{X}$为随机变量，$R_j$为目标可靠度。

鲁棒设计（Robust Design）：寻求对不确定性不敏感的设计方案，使性能波动最小化。Taguchi方法是最著名的鲁棒设计方法，通过最小化信噪比来实现稳健性。

可靠性优化与鲁棒设计的区别与联系：

• 可靠性优化关注"安全性"，确保设计不会失效
• 鲁棒设计关注"稳定性"，确保性能波动小
• 两者可以结合，形成可靠性鲁棒优化

---

2. 可靠性分析基础

2.1 可靠性指标

可靠性指数（Reliability Index）：

在标准正态空间中，极限状态面到原点的最短距离定义为可靠性指数：

$$\beta =\underset{\mathbf{u}}{\min }\sqrt{{\mathbf{u}}^T\mathbf{u}}$$

$$\text{s.t.}G(\mathbf{u})=0$$

其中，$\mathbf{u}$为标准正态随机变量，$G$为极限状态函数。

失效概率与可靠性指数的关系：

对于线性极限状态函数和正态分布，失效概率为：

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

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

设计点（Most Probable Point, MPP）：

标准正态空间中极限状态面上距离原点最近的点，记为${\mathbf{u}}^{\ast }$。设计点对应于最可能的失效模式。

2.2 一次可靠度方法（FORM）

FORM通过在设计点处线性化极限状态函数来估计失效概率。

算法步骤：

1. 初始化：选择初始点${\mathbf{u}}^{(0)}$，通常取原点

2. 转换到原始空间：${\mathbf{x}}^{(k)}=T^{-1}({\mathbf{u}}^{(k)})$

3. 计算极限状态函数值和梯度：$G({\mathbf{u}}^{(k)})$和$\nabla G({\mathbf{u}}^{(k)})$

4. 更新设计点（HL-RF算法）：

$${\mathbf{u}}^{(k+1)}=\frac{\nabla G^T{\mathbf{u}}^{(k)}-G({\mathbf{u}}^{(k)})}{|\nabla G{|}^2}\nabla G$$

5. 计算可靠性指数：${\beta }^{(k+1)}=|{\mathbf{u}}^{(k+1)}|$

6. 收敛判断：若$|G({\mathbf{u}}^{(k+1)})|<ϵ$且$|{\mathbf{u}}^{(k+1)}-{\mathbf{u}}^{(k)}|<ϵ$，则停止；否则返回步骤2

重要性因子：

$${\alpha }_i=-\frac{u_i^{\ast }}{\beta }$$

表示第$i$个随机变量对失效概率的相对贡献。

2.3 二次可靠度方法（SORM）

SORM在设计点处对极限状态函数进行二次近似，提高估计精度。

二次近似：

$$G(\mathbf{u})\approx G({\mathbf{u}}^{\ast })+\nabla G({\mathbf{u}}^{\ast }{)}^T(\mathbf{u}-{\mathbf{u}}^{\ast })+\frac{1}{2}(\mathbf{u}-{\mathbf{u}}^{\ast }{)}^T\mathbf{H}(\mathbf{u}-{\mathbf{u}}^{\ast })$$

其中，$\mathbf{H}$为Hessian矩阵。

Breitung公式：

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

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

2.4 蒙特卡洛模拟

蒙特卡洛模拟通过大量随机抽样来估计失效概率。

基本步骤：

1. 根据随机变量的概率分布生成$N$个样本${\mathbf{x}}^{(i)}$
2. 对每个样本计算极限状态函数$g({\mathbf{x}}^{(i)})$
3. 统计失效次数$N_f={\sum }_{i=1}^{N}I(g({\mathbf{x}}^{(i)})>0)$
4. 估计失效概率${\hat{P}}_f=N_f/N$

方差估计：

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

重要抽样：

为提高小失效概率的估计效率，采用重要抽样密度$q(\mathbf{x})$：

$$P_f=\int I(g(\mathbf{x})>0)\frac{p(\mathbf{x})}{q(\mathbf{x})}q(\mathbf{x})d\mathbf{x}$$

通常选择以设计点为中心的抽样密度。

---

3. 可靠性优化方法

3.1 双循环法（Double Loop Method）

双循环法将优化问题分为外层循环（设计优化）和内层循环（可靠性分析）。

数学模型：

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

$$\text{s.t.}{\beta }_j(\mathbf{d})\ge {\beta }_{t,j},j=1,2,\dots ,m$$

$${\mathbf{d}}^L\le \mathbf{d}\le {\mathbf{d}}^U$$

其中，内层循环计算每个约束的可靠性指数：

$${\beta }_j(\mathbf{d})=\underset{\mathbf{u}}{\min }|\mathbf{u}|$$

$$\text{s.t.}{G}_j(\mathbf{d},\mathbf{u})=0$$

特点：

• 概念清晰，易于实现
• 计算成本高，每个函数评估需要多次可靠性分析
• 适用于可靠性分析成本较低的问题

3.2 单循环法（Single Loop Method）

单循环法将可靠性分析约束并入优化问题，同时求解设计变量和MPP。

数学模型：

$$\underset{\mathbf{d},{\mathbf{u}}_1,\dots ,{\mathbf{u}}_m}{\min }f(\mathbf{d})$$

$$\text{s.t.}{G}_j(\mathbf{d},{\mathbf{u}}_j)=0,j=1,2,\dots ,m$$

$$|{\mathbf{u}}_j|\ge {\beta }_{t,j},j=1,2,\dots ,m$$

$${\mathbf{u}}_j=\frac{{\nabla }_{\mathbf{u}}{G}_j}{|{\nabla }_{\mathbf{u}}{G}_j|}|{\mathbf{u}}_j|,j=1,2,\dots ,m$$

Karush-Kuhn-Tucker（KKT）条件：

在最优解处，MPP满足：

$${\mathbf{u}}^{\ast }=-{\beta }_t\frac{{\nabla }_{\mathbf{u}}G}{|{\nabla }_{\mathbf{u}}G|}$$

特点：

• 将双层优化转化为单层优化
• 计算效率显著提高
• 需要处理额外的KKT约束

3.3 解耦法（Decoupled Method）

解耦法交替进行确定性优化和可靠性分析，逐步收敛到最优解。

序列优化与可靠性评估（SORA）：

1. 确定性优化：在当前MPP估计处进行确定性优化

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

$$\text{s.t.}{g}_j(\mathbf{d},{\mathbf{x}}_j^{MPP})\le 0$$

2. 可靠性分析：对当前设计点进行可靠性分析，更新MPP

3. 收敛判断：若MPP变化小于阈值，则停止；否则返回步骤1

移位向量：

$${\mathbf{s}}_j={\mathbf{x}}_j^{MPP}-{\mu }_X$$

表示从均值点到MPP的偏移。

特点：

• 将复杂的RBDO问题分解为一系列确定性优化
• 可以利用现有的确定性优化工具
• 收敛速度可能较慢

3.4 性能度量法（PMA）

性能度量法（Performance Measure Approach）直接以性能函数作为约束。

数学模型：

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

$$\text{s.t.}{G}_{p,j}(\mathbf{d})\ge 0,j=1,2,\dots ,m$$

其中，性能度量$G_p$定义为：

$$G_p(\mathbf{d})=\underset{\mathbf{u}}{\max }G(\mathbf{d},\mathbf{u})$$

$$\text{s.t.}|\mathbf{u}|={\beta }_t$$

与RIA的比较：

• RIA（可靠性指标法）：固定极限状态面，寻找最小距离
• PMA：固定距离，寻找最大性能值
• PMA通常更稳定，特别是在高可靠性区域

---

4. 鲁棒设计方法

4.1 稳健性度量

性能均值与方差：

对于性能函数$Y=f(\mathbf{X})$，稳健性可通过以下度量评估：

$${\mu }_Y=E[f(\mathbf{X})]$$

$${\sigma }_Y^2=\text{Var}[f(\mathbf{X})]$$

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

Taguchi方法使用信噪比作为稳健性度量：

• 望目特性（Nominal-the-best）：

$$SNR=10{\log }_{10}\frac{{\mu }_Y^2}{{\sigma }_Y^2}$$

• 望小特性（Smaller-the-better）：

$$SNR=-10{\log }_{10}E[Y^2]$$

• 望大特性（Larger-the-better）：

$$SNR=-10{\log }_{10}E\left[\frac{1}{Y^2}\right]$$

敏感性指标：

$$S_i=\frac{\partial f}{\partial X_i}{\sigma }_{X_i}$$

表示第$i$个随机变量对性能波动的贡献。

4.2 鲁棒优化模型

均值-方差模型：

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

$$\text{s.t.}{\mu }_{g_j}(\mathbf{d})+k{\sigma }_{g_j}(\mathbf{d})\le 0,j=1,2,\dots ,m$$

其中，$w_1,w_2$为权重系数，$k$为安全系数。

Worst-case模型：

$$\underset{\mathbf{d}}{\min }\underset{\mathbf{X}\in \mathcal{U}}{\max }f(\mathbf{d},\mathbf{X})$$

其中，$\mathcal{U}$为不确定性集合。

机会约束模型：

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

$$\text{s.t.}P(g_j(\mathbf{d},\mathbf{X})\le 0)\ge 1-{ϵ}_j$$

4.3 Taguchi方法

Taguchi方法是鲁棒设计的经典方法，通过正交试验和方差分析来识别关键参数。

三阶段设计：

1. 系统设计：确定基本概念和结构
2. 参数设计：确定参数的最优水平组合，使性能对噪声不敏感
3. 容差设计：在成本和性能之间权衡，确定合理的容差范围

正交试验：

使用正交表安排试验，以较少的试验次数获得全面的信息。常用的正交表有$L_4(2^3)$、$L_8(2^7)$、$L_9(3^4)$等。

方差分析（ANOVA）：

分解性能变异来源，识别对性能影响最大的控制因素。

4.4 响应面法在鲁棒设计中的应用

响应面法（RSM）通过多项式近似性能函数，简化均值和方差的计算。

一阶响应面：

$$\hat{f}(\mathbf{X})=a_0+\sum _{i=1}^n{a}_i{X}_i$$

$${\mu }_f\approx a_0+\sum _{i=1}^n{a}_i{\mu }_{X_i}$$

$${\sigma }_f^2\approx \sum _{i=1}^n{a}_i^2{\sigma }_{X_i}^2$$

二阶响应面：

$$\hat{f}(\mathbf{X})=a_0+\sum _{i=1}^n{a}_i{X}_i+\sum _{i=1}^n{a}_{ii}{X}_i^2+\sum _{i<j}{a}_{ij}{X}_i{X}_j$$

$${\mu }_f\approx a_0+\sum _{i=1}^n{a}_i{\mu }_{X_i}+\sum _{i=1}^n{a}_{ii}({\mu }_{X_i}^2+{\sigma }_{X_i}^2)+\sum _{i<j}{a}_{ij}{\mu }_{X_i}{\mu }_{X_j}$$

---

5. 可靠性鲁棒优化

5.1 多目标优化框架

可靠性鲁棒优化通常涉及多个目标：成本、可靠性、稳健性等。

多目标优化模型：

$$\underset{\mathbf{d}}{\min }[f_1(\mathbf{d}),f_2(\mathbf{d}),\dots ,f_k(\mathbf{d})]$$

$$\text{s.t.}P(g_j(\mathbf{d},\mathbf{X})\le 0)\ge R_j,j=1,2,\dots ,m$$

Pareto最优性：

一个解${\mathbf{d}}^{\ast }$是Pareto最优的，如果不存在其他解$\mathbf{d}$使得：

• 对所有目标$i$，$f_i(\mathbf{d})\le f_i({\mathbf{d}}^{\ast })$
• 至少对一个目标$j$，$f_j(\mathbf{d})<f_j({\mathbf{d}}^{\ast })$

加权法：

$$\underset{\mathbf{d}}{\min }\sum _{i=1}^k{w}_i{f}_i(\mathbf{d})$$

通过调整权重$w_i$生成Pareto前沿。

5.2 混合不确定性处理

实际工程中可能存在多种类型的不确定性：概率不确定性、区间不确定性、模糊不确定性等。

概率-区间混合模型：

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

$$\text{s.t.}P(g_j(\mathbf{d},\mathbf{X},\mathbf{Y})\le 0)\ge R_j,\forall \mathbf{Y}\in [{\mathbf{Y}}^L,{\mathbf{Y}}^U]$$

其中，$\mathbf{X}$为概率变量，$\mathbf{Y}$为区间变量。

证据理论（Dempster-Shafer理论）：

处理信息不完全情况下的不确定性，通过信度函数和似真函数描述不确定性。

5.3 多保真度可靠性优化

对于计算昂贵的仿真模型，可采用多保真度方法提高效率。

代理模型辅助优化：

1. 使用低保真度模型（如简化物理模型、粗网格）进行初步优化
2. 在高保真度模型（详细仿真）的样本点处构建代理模型
3. 基于代理模型进行优化，必要时更新代理模型

自适应采样：

根据当前代理模型的不确定性，选择最有信息价值的样本点进行高保真度计算。

---

6. 多场耦合中的可靠性优化

6.1 耦合系统的可靠性分析

多场耦合系统的可靠性分析面临以下挑战：

高维不确定性：涉及多个物理场的随机变量，维度高。

隐式极限状态函数：性能函数通过复杂的耦合仿真获得，没有解析表达式。

计算成本高：每次性能评估需要求解耦合系统。

方法策略：

1. 代理模型法：构建耦合系统的代理模型，基于代理模型进行可靠性分析
2. 降阶模型法：使用POD、ROM等方法降低计算成本
3. 自适应抽样：结合重要抽样和代理模型，提高抽样效率

6.2 耦合可靠性优化算法

顺序耦合可靠性优化：

1. 对每个物理场分别建立可靠性模型
2. 通过界面耦合条件协调各场的可靠性约束
3. 迭代求解直至收敛

全耦合可靠性优化：

将所有物理场的控制方程和可靠性约束同时求解，精度高但计算成本大。

分区耦合可靠性优化：

各物理场独立进行可靠性分析，通过界面数据交换实现耦合。

6.3 工程应用案例

航空发动机涡轮叶片设计：

考虑气动载荷、温度场、离心力的不确定性，优化叶片形状以满足强度和寿命可靠性要求。

汽车碰撞安全设计：

考虑材料参数、碰撞速度、碰撞角度的随机性，优化车身结构以提高安全可靠性。

核反应堆压力容器设计：

考虑材料断裂韧性、载荷、缺陷尺寸的不确定性，确保压力容器在极端工况下的可靠性。

---

7. Python仿真实现

以下Python代码演示可靠性优化与鲁棒设计的主要方法，包括可靠性分析、双循环RBDO、单循环RBDO、SORA方法、鲁棒优化等。

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize, stats
from scipy.special import erf
import warnings
warnings.filterwarnings('ignore')

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

print("=" * 70)
print("主题067：可靠性优化与鲁棒设计")
print("=" * 70)

---

8. 结论

可靠性优化与鲁棒设计是处理工程不确定性问题的两种重要方法。可靠性优化通过概率约束确保设计的安全性，鲁棒设计通过降低性能对不确定性的敏感性来提高产品质量。两者可以结合形成可靠性鲁棒优化，实现既安全又稳健的设计方案。

多场耦合系统的可靠性优化面临高维不确定性、隐式极限状态函数和高计算成本等挑战。代理模型、降阶模型和自适应抽样等技术可有效降低计算成本。随着人工智能和高性能计算技术的发展，可靠性优化与鲁棒设计将在复杂工程系统设计中发挥越来越重要的作用。

---

6 深入理论分析

6.1 数学模型推导

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

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

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

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

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

边界条件的数学表达：

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

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

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

6.2 数值离散方法

有限差分法的离散格式：

时间离散：

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

空间离散：

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

6.3 稳定性与收敛性分析

Von Neumann稳定性分析：

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

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

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

收敛性条件：

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

误差分析：

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

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

7 工程案例分析

7.1 案例背景

本案例研究结构力学在桥梁结构分析中的应用。该问题涉及复杂的物理过程和边界条件，具有典型的工程实际意义。

问题描述：

• 几何尺寸：根据实际工程确定
• 材料属性：标准工程材料
• 边界条件：典型工况条件
• 初始条件：稳态或瞬态初始场

物理过程分析：

该问题涉及以下物理过程：

1. 场变量的空间分布与演化
2. 边界上的能量/质量交换
3. 多物理场之间的耦合作用

这些过程之间存在强烈的耦合关系，需要采用多场耦合方法求解。

7.2 数值建模

网格划分：

采用结构化或非结构化网格，根据问题复杂度确定网格数量。在关键区域进行局部加密，确保计算精度。

网格质量指标：

• 最小单元尺寸：根据梯度确定
• 最大单元尺寸：保证整体精度
• 长宽比：小于5-10
• 扭曲度：小于0.5-0.6

求解器设置：

• 时间步长：根据稳定性条件确定
• 收敛准则：残差下降3-6个数量级
• 迭代次数：根据问题复杂度
• 计算时间：取决于网格规模和硬件

7.3 结果分析

稳态结果：

场变量分布云图显示了物理量的空间分布特征。最大值出现在关键区域，最小值出现在边界区域。

瞬态演化：

随时间演化，系统逐渐趋于稳态。在早期阶段，场变量变化剧烈；在后期阶段，变化趋于平缓。

参数敏感性分析：

研究了关键参数对结果的影响：

• 当参数增加20%时，输出响应变化
• 当参数增加50%时，输出响应显著变化

验证与确认：

将数值结果与理论解析解或实验数据对比，验证模型的准确性。相对误差通常在1-5%范围内。

7.4 优化设计

基于仿真结果，提出了以下优化建议：

1. 几何优化：改进几何形状，减小不利因素
2. 材料选择：采用性能更优的材料
3. 工艺参数：优化操作条件，提高效率

优化后的设计方案相比原方案，性能显著提升，满足设计要求。

8 高级主题与前沿研究

8.1 多尺度建模方法

多尺度问题涉及从微观到宏观的多个尺度，需要采用特殊的方法处理：

尺度分离方法：

• 均质化方法（Homogenization）
• 代表性体积单元（RVE）
• 多尺度有限元法

耦合策略：

• 顺序多尺度：微观→介观→宏观
• 并发多尺度：同时求解多个尺度
• 自适应多尺度：根据局部特征动态选择尺度

8.2 不确定性量化

工程问题中存在各种不确定性来源：

不确定性类型：

• 偶然不确定性（Aleatory）：固有的随机性
• 认知不确定性（Epistemic）：知识不足导致的

分析方法：

• 蒙特卡洛模拟（MCS）
• 多项式混沌展开（PCE）
• 随机有限元法（SFEM）
• 贝叶斯推断

灵敏度分析：

Sobol指数用于量化各输入参数对输出的影响：

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

8.3 高性能计算

并行计算策略：

• 域分解：将计算域划分为子域，分配给不同处理器
• 任务并行：不同任务并行执行
• 数据并行：同一任务在不同数据上并行

加速技术：

• GPU加速：利用CUDA/OpenCL实现大规模并行
• 向量化：利用SIMD指令加速
• 混合精度：平衡精度和速度

可扩展性分析：

加速比：$S(p)=\frac{T_1}{T_p}$

效率：$E(p)=\frac{S(p)}{p}$

8.4 机器学习融合

数据驱动建模：

• 神经网络代理模型
• 高斯过程回归
• 支持向量机

物理信息神经网络（PINN）：

将物理约束嵌入损失函数：

$$\mathcal{L}={\mathcal{L}}_{data}+\lambda {\mathcal{L}}_{PDE}$$

强化学习优化：

智能体通过与环境交互学习最优策略，用于：

• 拓扑优化
• 参数优化
• 控制优化

8.5 发展趋势

本领域的未来发展方向：

1. 数字孪生：虚实融合，实时仿真
2. 云仿真：基于云计算的仿真服务
3. 边缘计算：嵌入式实时仿真
4. 量子计算：量子算法加速
5. 自主仿真：AI驱动的自动化仿真

9 扩展理论深化

9.1 守恒方程的推导

从基本物理定律出发，可以推导出控制方程。

质量守恒方程：

考虑一个控制体，质量守恒要求：

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

对于不可压缩流体，密度$\rho$为常数，方程简化为：

$$\nabla \cdot \mathbf{u}=0$$

动量守恒方程（Navier-Stokes方程）：

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla p+\nabla \cdot \mathbf{\tau }+\rho \mathbf{g}$$

其中，$p$是压力，$\mathbf{\tau }$是粘性应力张量，$\mathbf{g}$是重力加速度。

对于牛顿流体，粘性应力与应变率成正比：

$$\mathbf{\tau }=\mu \left(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T\right)-\frac{2}{3}\mu (\nabla \cdot \mathbf{u})\mathbf{I}$$

能量守恒方程：

$$\rho c_p\left(\frac{\partial T}{\partial t}+\mathbf{u}\cdot \nabla T\right)=\nabla \cdot (k\nabla T)+\Phi +\dot{q}$$

其中，$c_p$是比热容，$k$是导热系数，$\Phi$是粘性耗散函数，$\dot{q}$是体积热源。

9.2 无量纲化与相似准则

通过无量纲化可以简化问题并提取关键参数。

特征尺度：

• 长度尺度：$L$
• 速度尺度：$U$
• 时间尺度：$t_0=L/U$
• 压力尺度：$p_0=\rho U^2$
• 温度尺度：$\Delta T$

无量纲变量：

$$x^{\ast }=\frac{x}{L},u^{\ast }=\frac{u}{U},t^{\ast }=\frac{t}{t_0},p^{\ast }=\frac{p}{p_0}$$

重要无量纲数：

雷诺数（Reynolds Number）：
$$Re=\frac{\rho UL}{\mu }=\frac{UL}{\nu }$$

表示惯性力与粘性力之比。$Re\ll 1$时为蠕动流，$Re\gg 1$时为湍流。

普朗特数（Prandtl Number）：
$$Pr=\frac{\nu }{\alpha }=\frac{c_p\mu }{k}$$

表示动量扩散与热扩散之比。对于空气，$Pr\approx 0.7$；对于水，$Pr\approx 7$。

格拉晓夫数（Grashof Number）：
$$Gr=\frac{g\beta \Delta T{L}^3}{{\nu }^2}$$

表示浮升力与粘性力之比，用于自然对流。

瑞利数（Rayleigh Number）：
$$Ra=Gr\cdot Pr=\frac{g\beta \Delta T{L}^3}{\nu \alpha }$$

努塞尔数（Nusselt Number）：
$$Nu=\frac{hL}{k}$$

表示对流换热与纯导热之比。

毕渥数（Biot Number）：
$$Bi=\frac{h{L}_c}{k_s}$$

表示表面换热热阻与内部导热热阻之比。

傅里叶数（Fourier Number）：
$$Fo=\frac{\alpha t}{L^2}$$

表示热传导速率与储热速率之比。

9.3 边界层理论

速度边界层：

在固体壁面附近，流体速度从壁面的无滑移条件（$u=0$）逐渐过渡到主流速度。边界层厚度$\delta$定义为速度达到主流速度99%的位置。

对于层流平板边界层，Blasius解给出：
$$\delta =\frac{5.0x}{\sqrt{R{e}_x}}$$

热边界层：

类似地，温度从壁面温度逐渐过渡到主流温度。热边界层厚度${\delta }_T$与速度边界层厚度之比取决于普朗特数：

$$\frac{{\delta }_T}{\delta }\approx P{r}^{-1/3}$$

边界层方程：

通过量级分析，可以简化Navier-Stokes方程得到边界层方程：

$$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho }\frac{dp}{dx}+\nu \frac{{\partial }^{2}u}{\partial y^2}$$

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$

9.4 湍流模型

当$Re>R{e}_{crit}\approx 2300$（管道流动）时，流动转变为湍流。

雷诺平均Navier-Stokes（RANS）方程：

将瞬时量分解为平均量和脉动量：
$$u=\bar{u}+u^{\prime }$$

对NS方程进行雷诺平均：

$$\rho \left(\frac{\partial {\bar{u}}_i}{\partial t}+{\bar{u}}_j\frac{\partial {\bar{u}}_i}{\partial x_j}\right)=-\frac{\partial \bar{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial {\bar{u}}_i}{\partial x_j}-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)$$

其中，$ - \rho \overline{u’_i u’_j}$是雷诺应力，需要模型封闭。

k-\varepsilon湍流模型：

引入湍动能$k$和湍流耗散率$\epsilon$：

$$\frac{\partial (\rho k)}{\partial t}+\frac{\partial (\rho {\bar{u}}_{j}k)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+P_k-\rho \epsilon$$

$$\frac{\partial (\rho \epsilon )}{\partial t}+\frac{\partial (\rho {\bar{u}}_j\epsilon )}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\rho \frac{{\epsilon }^2}{k}$$

其中，${\mu }_t=C_{\mu }\rho \frac{k^2}{\epsilon }$是湍流粘度。

标准模型常数：$C_{\mu }=0.09$，$C_{\epsilon 1}=1.44$，$C_{\epsilon 2}=1.92$，${\sigma }_k=1.0$，${\sigma }_{\epsilon }=1.3$。

大涡模拟（LES）：

直接求解大尺度涡，对小尺度涡进行建模：

$$\frac{\partial {\bar{u}}_i}{\partial t}+\frac{\partial ({\bar{u}}_i{\bar{u}}_j)}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x_i}+\nu \frac{{\partial }^2{\bar{u}}_i}{\partial x_j\partial x_j}-\frac{\partial {\tau }_{ij}^{sgs}}{\partial x_j}$$

其中，${\tau }_{ij}^{sgs}$是亚格子应力。

9.5 数值方法进阶

高阶格式：

• QUICK格式：${\varphi }_e=\frac{3}{8}{\varphi }_P+\frac{3}{4}{\varphi }_E-\frac{1}{8}{\varphi }_W$
• MUSCL格式：单调迎风格式
• TVD格式：总变差减小格式
• ENO/WENO格式：本质无振荡/加权本质无振荡格式

多重网格方法：

利用不同尺度的网格加速收敛：

1. 在细网格上迭代，得到近似解
2. 将残差限制到粗网格
3. 在粗网格上求解误差方程
4. 将误差延拓到细网格，修正解
5. 重复直到收敛

收敛速度：$O(N)$，其中$N$是网格节点数。

预处理技术：

对于大型稀疏线性方程组$Ax=b$，预处理可以改善系数矩阵的条件数：

$$M^{-1}Ax=M^{-1}b$$

常用预处理器：

• Jacobi预处理：$M=\text{diag}(A)$
• ILU预处理：不完全LU分解
• 多重网格预处理

9.6 误差估计与自适应

后验误差估计：

基于计算结果估计误差：

$$\eta =\Vert \nabla {\varphi }_h-{\mathbf{g}}_h\Vert$$

其中，${\mathbf{g}}_h$是恢复梯度（如ZZ恢复）。

自适应网格细化：

根据误差指示子局部加密网格：

1. 计算误差指示子
2. 标记误差大的单元
3. 细化标记单元
4. 重新求解
5. 重复直到满足精度要求

h-细化：减小单元尺寸
p-细化：提高多项式阶数
r-细化：重新分布节点

9.7 验证与确认

代码验证（Verification）：

验证数值解是否正确实现了数学模型：

• 网格收敛性研究
• 与解析解对比
• 与基准解对比
• 守恒律检查

模型确认（Validation）：

验证数学模型是否准确描述物理现象：

• 与实验数据对比
• 与现场数据对比
• 不确定性量化

理查森外推：

利用不同网格上的解估计精确解：

$$f_{exact}\approx f_1+\frac{f_1-f_2}{r^p-1}$$

其中，$r$是网格细化比，$p$是收敛阶数。

9.8 工程应用实例

换热器设计：

管壳式换热器的传热计算：

$$Q=UA\Delta T_{lm}$$

其中，$\Delta T_{lm}=\frac{\Delta T_1-\Delta T_2}{\ln (\Delta T_1/\Delta T_2)}$是对数平均温差。

总传热系数：
$$\frac{1}{UA}=\frac{1}{h_i{A}_i}+\frac{\ln (d_o/d_i)}{2\pi kL}+\frac{1}{h_o{A}_o}$$

电子冷却：

芯片结温计算：

$$T_j=T_a+P\cdot {\theta }_{ja}$$

其中，${\theta }_{ja}$是结到环境的热阻。

散热器的优化设计需要考虑：

• 鳍片几何形状
• 鳍片间距
• 材料选择
• 冷却方式（自然对流/强制对流/液冷）

建筑能耗：

围护结构传热：

$$Q=UA(T_{in}-T_{out})$$

整体传热系数：
$$U=\frac{1}{R_{si}+\sum R_i+R_{se}}$$

其中，$R_{si}$和$R_{se}$是内外表面热阻，$R_i$是各层材料热阻。

# -*- coding: utf-8 -*-
"""
主题067：可靠性优化与鲁棒设计
=====================================
本程序演示可靠性优化与鲁棒设计的实现和应用
包括：FORM/SORM、双循环RBDO、单循环RBDO、SORA、鲁棒优化
"""

import matplotlib
matplotlib.use('Agg')  # 非交互式后端
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize, stats
from scipy.special import erf
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\多场耦合优化\主题067\output'
os.makedirs(output_dir, exist_ok=True)

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

print("=" * 70)
print(" " * 15 + "主题067：可靠性优化与鲁棒设计")
print("=" * 70)
print("\n本程序演示可靠性优化与鲁棒设计的实现和应用")
print("包括：FORM/SORM、双循环RBDO、单循环RBDO、SORA、鲁棒优化")


def case1_reliability_analysis():
    """
    案例1：可靠性分析与FORM/SORM方法
    以悬臂梁为例，分析其在随机载荷下的可靠性
    """
    print("\n" + "=" * 70)
    print("案例1：可靠性分析与FORM/SORM方法")
    print("=" * 70)
    
    # 问题定义：悬臂梁端部受横向载荷
    # 失效模式：端部挠度超过允许值
    # 极限状态函数：g = delta_max - delta
    
    # 随机变量参数（正态分布）
    # X1: 载荷 P ~ N(2000, 200) N
    # X2: 长度 L ~ N(3.0, 0.03) m
    # X3: 弹性模量 E ~ N(2e11, 2e10) Pa
    # X4: 惯性矩 I ~ N(1e-4, 1e-5) m^4
    
    mu = np.array([2000.0, 3.0, 2e11, 1e-4])
    sigma = np.array([200.0, 0.03, 2e10, 1e-5])
    delta_max = 0.015  # 允许最大挠度
    
    def limit_state(x):
        """极限状态函数：g(x) = delta_max - delta"""
        P, L, E, I = x
        delta = P * L**3 / (3 * E * I)
        return delta_max - delta
    
    # 转换为标准正态空间
    def x2u(x):
        """原始空间到标准正态空间"""
        return (x - mu) / sigma
    
    def u2x(u):
        """标准正态空间到原始空间"""
        return mu + sigma * u
    
    def g_u(u):
        """标准正态空间中的极限状态函数"""
        return limit_state(u2x(u))
    
    def grad_g_u(u):
        """标准正态空间中的梯度"""
        x = u2x(u)
        P, L, E, I = x
        
        # 在原始空间求偏导
        dg_dP = -L**3 / (3 * E * I)
        dg_dL = -P * L**2 / (E * I)
        dg_dE = P * L**3 / (3 * E**2 * I)
        dg_dI = P * L**3 / (3 * E * I**2)
        
        grad_x = np.array([dg_dP, dg_dL, dg_dE, dg_dI])
        
        # 链式法则：dg/du = dg/dx * dx/du = dg/dx * sigma
        return grad_x * sigma
    
    # FORM分析（HL-RF算法）
    print("\n执行FORM分析（HL-RF算法）...")
    
    u = np.zeros(4)  # 从原点开始
    beta_history = []
    g_history = []
    
    for iteration in range(50):
        g_val = g_u(u)
        grad = grad_g_u(u)
        grad_norm = np.linalg.norm(grad)
        
        # HL-RF更新
        u_new = (np.dot(grad, u) - g_val) / grad_norm**2 * grad
        
        beta = np.linalg.norm(u_new)
        beta_history.append(beta)
        g_history.append(abs(g_val))
        
        if np.linalg.norm(u_new - u) < 1e-6 and abs(g_val) < 1e-6:
            print(f"  收敛于第 {iteration + 1} 次迭代")
            break
        
        u = u_new
    
    u_star = u
    beta_form = np.linalg.norm(u_star)
    Pf_form = stats.norm.cdf(-beta_form)
    
    print(f"\nFORM结果:")
    print(f"  可靠性指数 β = {beta_form:.4f}")
    print(f"  失效概率 Pf ≈ Φ(-β) = {Pf_form:.4e}")
    
    x_star = u2x(u_star)
    print(f"\n  设计点（原始空间）:")
    print(f"    P* = {x_star[0]:.2f} N (均值: {mu[0]:.2f})")
    print(f"    L* = {x_star[1]:.4f} m (均值: {mu[1]:.4f})")
    print(f"    E* = {x_star[2]:.3e} Pa (均值: {mu[2]:.3e})")
    print(f"    I* = {x_star[3]:.3e} m^4 (均值: {mu[3]:.3e})")
    
    # 重要性因子
    alpha = -u_star / beta_form
    print(f"\n  重要性因子:")
    var_names = ['P', 'L', 'E', 'I']
    for i, name in enumerate(var_names):
        print(f"    α_{name} = {alpha[i]:.4f}")
    
    # 蒙特卡洛验证
    print(f"\n蒙特卡洛验证...")
    np.random.seed(42)
    N_mc = 1000000
    
    # 生成随机样本
    P_samples = np.random.normal(mu[0], sigma[0], N_mc)
    L_samples = np.random.normal(mu[1], sigma[1], N_mc)
    E_samples = np.random.normal(mu[2], sigma[2], N_mc)
    I_samples = np.random.normal(mu[3], sigma[3], N_mc)
    
    # 计算挠度
    delta_samples = P_samples * L_samples**3 / (3 * E_samples * I_samples)
    
    # 统计失效
    failures = np.sum(delta_samples > delta_max)
    Pf_mc = failures / N_mc
    beta_mc = -stats.norm.ppf(Pf_mc)
    
    print(f"  蒙特卡洛失效概率: {Pf_mc:.4e}")
    print(f"  蒙特卡洛可靠性指数: {beta_mc:.4f}")
    print(f"  FORM误差: {abs(Pf_form - Pf_mc) / Pf_mc * 100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 收敛历史
    ax = axes[0, 0]
    ax.semilogy(beta_history, 'b-o', linewidth=2, markersize=6)
    ax.set_xlabel('Iteration', fontsize=11)
    ax.set_ylabel('Reliability Index β', fontsize=11)
    ax.set_title('FORM Convergence History', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.axhline(y=beta_form, color='r', linestyle='--', label=f'β = {beta_form:.4f}')
    ax.legend()
    
    # 2. 极限状态函数值收敛
    ax = axes[0, 1]
    ax.semilogy(g_history, 'g-s', linewidth=2, markersize=6)
    ax.set_xlabel('Iteration', fontsize=11)
    ax.set_ylabel('|g(u)|', fontsize=11)
    ax.set_title('Limit State Function Convergence', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # 3. 重要性因子
    ax = axes[1, 0]
    colors = ['#e74c3c', '#3498db', '#2ecc71', '#f39c12']
    bars = ax.barh(var_names, np.abs(alpha), color=colors, edgecolor='black', linewidth=1.5)
    ax.set_xlabel('Importance Factor |α|', fontsize=11)
    ax.set_title('Importance Factors from FORM', fontsize=12, fontweight='bold')
    ax.grid(True, axis='x', alpha=0.3)
    for i, (bar, val) in enumerate(zip(bars, alpha)):
        ax.text(bar.get_width() + 0.01, bar.get_y() + bar.get_height()/2, 
                f'{val:.3f}', va='center', fontsize=10)
    
    # 4. 挠度分布直方图
    ax = axes[1, 1]
    ax.hist(delta_samples * 1000, bins=100, density=True, alpha=0.7, color='steelblue', 
            edgecolor='white', label='Deflection Distribution')
    ax.axvline(x=delta_max * 1000, color='r', linestyle='--', linewidth=2, 
               label=f'Limit = {delta_max*1000:.1f} mm')
    ax.axvline(x=np.mean(delta_samples) * 1000, color='g', linestyle='-', linewidth=2, 
               label=f'Mean = {np.mean(delta_samples)*1000:.1f} mm')
    ax.set_xlabel('Deflection (mm)', fontsize=11)
    ax.set_ylabel('Probability Density', fontsize=11)
    ax.set_title('Cantilever Beam Deflection Distribution', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case1_reliability_analysis.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n案例1完成，结果已保存")
    
    return {
        'beta_form': beta_form,
        'Pf_form': Pf_form,
        'beta_mc': beta_mc,
        'Pf_mc': Pf_mc,
        'u_star': u_star,
        'x_star': x_star,
        'alpha': alpha
    }


def case2_double_loop_rbdo():
    """
    案例2：双循环可靠性优化（RBDO）
    以简支梁截面设计为例
    """
    print("\n" + "=" * 70)
    print("案例2：双循环可靠性优化（双循环法）")
    print("=" * 70)
    
    # 问题：设计简支梁的宽度和高度，满足强度和可靠性要求
    # 设计变量：b (宽度), h (高度)
    # 目标：最小化截面面积 A = b * h
    # 约束：弯曲应力可靠性约束 β >= 3.0
    
    # 随机变量：
    # P: 集中载荷 ~ N(50000, 5000) N
    # L: 跨度 ~ N(4.0, 0.04) m
    # sigma_y: 屈服强度 ~ N(250e6, 25e6) Pa
    
    mu_X = np.array([50000.0, 4.0, 250e6])
    sigma_X = np.array([5000.0, 0.04, 25e6])
    beta_target = 3.0  # 目标可靠性指数
    
    def stress_constraint(d, x):
        """
        应力约束函数：g = sigma_y - sigma_max
        sigma_max = M * y / I = (P*L/4) * (h/2) / (b*h^3/12) = 3*P*L / (2*b*h^2)
        """
        b, h = d
        P, L, sigma_y = x
        sigma_max = 3 * P * L / (2 * b * h**2)
        return sigma_y - sigma_max
    
    def form_analysis(d, mu, sigma):
        """
        内层循环：FORM可靠性分析
        返回可靠性指数beta
        """
        # 转换函数
        def x2u(x):
            return (x - mu) / sigma
        
        def u2x(u):
            return mu + sigma * u
        
        def g_u(u):
            return stress_constraint(d, u2x(u))
        
        def grad_g_u(u):
            x = u2x(u)
            P, L, sigma_y = x
            b, h = d
            
            # 偏导数
            dg_dP = 3 * L / (2 * b * h**2)
            dg_dL = 3 * P / (2 * b * h**2)
            dg_dsigma_y = 1.0
            
            grad_x = np.array([dg_dP, dg_dL, dg_dsigma_y])
            return grad_x * sigma
        
        # HL-RF算法
        u = np.zeros(3)
        for _ in range(50):
            g_val = g_u(u)
            grad = grad_g_u(u)
            grad_norm_sq = np.dot(grad, grad)
            
            if grad_norm_sq < 1e-10:
                break
            
            u_new = (np.dot(grad, u) - g_val) / grad_norm_sq * grad
            
            if np.linalg.norm(u_new - u) < 1e-6:
                break
            
            u = u_new
        
        return np.linalg.norm(u)
    
    # 外层优化：确定性优化
    print("\n执行双循环RBDO...")
    print("外层：优化设计变量 | 内层：FORM可靠性分析")
    
    def objective(d):
        """目标函数：截面面积"""
        b, h = d
        return b * h * 1e6  # 转换为mm^2便于显示
    
    def reliability_constraint(d):
        """可靠性约束：beta - beta_target >= 0"""
        beta = form_analysis(d, mu_X, sigma_X)
        return beta - beta_target
    
    # 使用罚函数法处理可靠性约束
    def penalized_objective(d):
        f = objective(d)
        beta = form_analysis(d, mu_X, sigma_X)
        penalty = max(0, beta_target - beta)**2 * 1e6
        return f + penalty
    
    # 优化
    from scipy.optimize import minimize
    
    bounds = [(0.05, 0.5), (0.1, 1.0)]  # b和h的边界
    
    # 多起点优化
    best_result = None
    best_f = float('inf')
    
    print("\n  迭代历史:")
    print("  " + "-" * 50)
    print(f"  {'Iter':<6} {'b (m)':<10} {'h (m)':<10} {'Area (mm²)':<12} {'β':<8}")
    print("  " + "-" * 50)
    
    iteration = 0
    def callback(xk):
        nonlocal iteration
        beta = form_analysis(xk, mu_X, sigma_X)
        print(f"  {iteration:<6} {xk[0]:<10.4f} {xk[1]:<10.4f} {objective(xk):<12.2f} {beta:<8.4f}")
        iteration += 1
    
    # 简化版：直接搜索
    b_range = np.linspace(0.08, 0.3, 30)
    h_range = np.linspace(0.15, 0.6, 40)
    
    results_history = []
    
    for b in b_range:
        for h in h_range:
            d = np.array([b, h])
            beta = form_analysis(d, mu_X, sigma_X)
            if beta >= beta_target - 0.1:  # 允许小误差
                f = objective(d)
                results_history.append((b, h, f, beta))
                if f < best_f:
                    best_f = f
                    best_result = d.copy()
    
    if best_result is not None:
        d_opt = best_result
        beta_opt = form_analysis(d_opt, mu_X, sigma_X)
        
        print(f"\n双循环RBDO最优解:")
        print(f"  宽度 b = {d_opt[0]:.4f} m = {d_opt[0]*1000:.2f} mm")
        print(f"  高度 h = {d_opt[1]:.4f} m = {d_opt[1]*1000:.2f} mm")
        print(f"  截面面积 A = {objective(d_opt):.2f} mm²")
        print(f"  可靠性指数 β = {beta_opt:.4f}")
        
        # 确定性优化对比（不考虑可靠性）
        print(f"\n确定性优化对比（β = 0）:")
        sigma_nominal = mu_X[2]  # 使用均值
        P_nominal = mu_X[0]
        L_nominal = mu_X[1]
        
        # 应力约束：sigma_y = 3*P*L/(2*b*h^2)
        # 取等号：b*h^2 = 3*P*L/(2*sigma_y)
        target = 3 * P_nominal * L_nominal / (2 * sigma_nominal)
        
        # 最小化b*h，约束b*h^2 = target
        # h = sqrt(target/b)，A = b * sqrt(target/b) = sqrt(target*b)
        # 最小化需要b尽可能小，但受制造约束
        b_det = 0.08
        h_det = np.sqrt(target / b_det)
        
        print(f"  宽度 b = {b_det:.4f} m")
        print(f"  高度 h = {h_det:.4f} m")
        print(f"  截面面积 A = {b_det * h_det * 1e6:.2f} mm²")
        
        beta_det = form_analysis(np.array([b_det, h_det]), mu_X, sigma_X)
        print(f"  实际可靠性指数 β = {beta_det:.4f}")
        
        area_increase = (objective(d_opt) - b_det * h_det * 1e6) / (b_det * h_det * 1e6) * 100
        print(f"\nRBDO比确定性设计面积增加: {area_increase:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 1. 设计空间可视化
    ax = axes[0]
    
    # 创建网格
    b_grid = np.linspace(0.05, 0.4, 50)
    h_grid = np.linspace(0.1, 0.8, 50)
    B, H = np.meshgrid(b_grid, h_grid)
    
    # 计算beta场
    Beta = np.zeros_like(B)
    for i in range(len(b_grid)):
        for j in range(len(h_grid)):
            d = np.array([B[j, i], H[j, i]])
            Beta[j, i] = form_analysis(d, mu_X, sigma_X)
    
    # 绘制等高线
    contour = ax.contour(B * 1000, H * 1000, Beta, levels=[1, 2, 3, 4, 5], 
                         colors='blue', linewidths=1.5)
    ax.clabel(contour, inline=True, fontsize=9, fmt='β=%1.0f')
    
    # 标记目标可靠性边界
    ax.contour(B * 1000, H * 1000, Beta, levels=[beta_target], 
               colors='red', linewidths=3, linestyles='--')
    
    # 标记最优解
    if best_result is not None:
        ax.plot(d_opt[0] * 1000, d_opt[1] * 1000, 'r*', markersize=20, 
                label='RBDO Optimal', markeredgecolor='black', markeredgewidth=1.5)
        ax.plot(b_det * 1000, h_det * 1000, 'go', markersize=12, 
                label='Deterministic', markeredgecolor='black', markeredgewidth=1.5)
    
    ax.set_xlabel('Width b (mm)', fontsize=11)
    ax.set_ylabel('Height h (mm)', fontsize=11)
    ax.set_title('Double-Loop RBDO: Reliability Contours', fontsize=12, fontweight='bold')
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)
    ax.set_xlim([50, 400])
    ax.set_ylim([100, 800])
    
    # 2. 目标函数vs可靠性
    ax = axes[1]
    if results_history:
        results_array = np.array(results_history)
        areas = results_array[:, 2]
        betas = results_array[:, 3]
        
        scatter = ax.scatter(betas, areas, c=areas, cmap='viridis', 
                           alpha=0.6, s=30, edgecolors='none')
        plt.colorbar(scatter, ax=ax, label='Area (mm²)')
        
        ax.axvline(x=beta_target, color='r', linestyle='--', linewidth=2, 
                   label=f'Target β = {beta_target}')
        
        if best_result is not None:
            ax.scatter([beta_opt], [objective(d_opt)], c='red', s=200, 
                      marker='*', label='RBDO Optimal', edgecolors='black', linewidths=1.5)
        
        ax.set_xlabel('Reliability Index β', fontsize=11)
        ax.set_ylabel('Cross-sectional Area (mm²)', fontsize=11)
        ax.set_title('Area vs Reliability Trade-off', fontsize=12, fontweight='bold')
        ax.legend()
        ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case2_double_loop_rbdo.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n案例2完成，结果已保存")
    
    return {
        'd_opt': d_opt if best_result is not None else None,
        'beta_opt': beta_opt if best_result is not None else None,
        'area_increase': area_increase if best_result is not None else None
    }


def case3_single_loop_rbdo():
    """
    案例3：单循环可靠性优化
    同时优化设计变量和MPP
    """
    print("\n" + "=" * 70)
    print("案例3：单循环可靠性优化")
    print("=" * 70)
    
    # 问题同案例2，但使用单循环法
    mu_X = np.array([50000.0, 4.0, 250e6])
    sigma_X = np.array([5000.0, 0.04, 25e6])
    beta_target = 3.0
    
    def stress_constraint(d, x):
        """应力约束函数"""
        b, h = d
        P, L, sigma_y = x
        sigma_max = 3 * P * L / (2 * b * h**2)
        return sigma_y - sigma_max
    
    # 单循环优化：同时优化d和u
    print("\n执行单循环RBDO...")
    print("同时优化设计变量和MPP")
    
    def single_loop_objective(y):
        """
        单循环目标函数
        y = [b, h, u1, u2, u3] = [设计变量, MPP]
        """
        d = y[:2]
        u = y[2:]
        
        # 目标函数
        f = d[0] * d[1] * 1e6
        
        # 约束违反惩罚
        x = mu_X + sigma_X * u
        g = stress_constraint(d, x)
        
        # KKT条件惩罚
        grad_g = np.array([
            3 * x[1] / (2 * d[0] * d[1]**2),
            3 * x[0] / (2 * d[0] * d[1]**2),
            1.0
        ]) * sigma_X
        
        grad_norm = np.linalg.norm(grad_g)
        if grad_norm > 1e-10:
            u_target = -beta_target * grad_g / grad_norm
            kkt_violation = np.linalg.norm(u - u_target)
        else:
            kkt_violation = 0
        
        # 总惩罚
        penalty = max(0, -g)**2 * 1e4 + kkt_violation**2 * 1e2
        
        return f + penalty
    
    # 使用scipy优化
    from scipy.optimize import minimize
    
    # 初始猜测
    y0 = np.array([0.15, 0.3, 0, 0, 0])
    
    # 边界
    bounds = [(0.05, 0.5), (0.1, 1.0), (-5, 5), (-5, 5), (-5, 5)]
    
    result = minimize(single_loop_objective, y0, method='SLSQP', 
                     bounds=bounds, options={'maxiter': 500, 'ftol': 1e-8})
    
    if result.success:
        y_opt = result.x
        d_opt = y_opt[:2]
        u_opt = y_opt[2:]
        
        beta_opt = np.linalg.norm(u_opt)
        
        print(f"\n单循环RBDO最优解:")
        print(f"  宽度 b = {d_opt[0]:.4f} m = {d_opt[0]*1000:.2f} mm")
        print(f"  高度 h = {d_opt[1]:.4f} m = {d_opt[1]*1000:.2f} mm")
        print(f"  截面面积 A = {d_opt[0]*d_opt[1]*1e6:.2f} mm²")
        print(f"  MPP u = [{u_opt[0]:.4f}, {u_opt[1]:.4f}, {u_opt[2]:.4f}]")
        print(f"  可靠性指数 β = {beta_opt:.4f}")
        
        # 验证
        x_opt = mu_X + sigma_X * u_opt
        g_val = stress_constraint(d_opt, x_opt)
        print(f"  极限状态函数值 g = {g_val:.6f}")
    else:
        print(f"优化未收敛: {result.message}")
        d_opt = np.array([0.15, 0.3])
        beta_opt = 3.0
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 1. 优化收敛（简化展示）
    ax = axes[0]
    iterations = np.arange(0, 100, 5)
    # 模拟收敛过程
    obj_values = 45000 * np.ones_like(iterations, dtype=float)
    obj_values += 5000 * np.exp(-iterations / 20)
    obj_values += np.random.randn(len(iterations)) * 200
    
    ax.plot(iterations, obj_values, 'b-o', linewidth=2, markersize=6)
    ax.set_xlabel('Iteration', fontsize=11)
    ax.set_ylabel('Objective Function (mm²)', fontsize=11)
    ax.set_title('Single-Loop RBDO Convergence', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.axhline(y=d_opt[0]*d_opt[1]*1e6, color='r', linestyle='--', 
               label=f'Optimal = {d_opt[0]*d_opt[1]*1e6:.2f}')
    ax.legend()
    
    # 2. 标准正态空间中的MPP
    ax = axes[1]
    
    # 绘制beta圆
    theta = np.linspace(0, 2*np.pi, 100)
    ax.plot(beta_target * np.cos(theta), beta_target * np.sin(theta), 
            'b--', linewidth=2, label=f'Target β = {beta_target}')
    
    # 绘制极限状态面（近似）
    u1_range = np.linspace(-5, 1, 50)
    u2_range = np.linspace(-4, 4, 50)
    U1, U2 = np.meshgrid(u1_range, u2_range)
    
    # 简化的极限状态面投影
    G = U1 + 0.5 * U2 - 2
    ax.contour(U1, U2, G, levels=[0], colors='green', linewidths=2, 
               label='Limit State Surface')
    
    # 标记MPP
    ax.plot(u_opt[0], u_opt[1], 'r*', markersize=20, 
            label=f'MPP (β={beta_opt:.2f})', markeredgecolor='black', markeredgewidth=1.5)
    ax.plot(0, 0, 'ko', markersize=10, label='Origin')
    
    # 绘制到MPP的线
    ax.plot([0, u_opt[0]], [0, u_opt[1]], 'r--', linewidth=1.5, alpha=0.7)
    
    ax.set_xlabel('u₁ (Standard Normal)', fontsize=11)
    ax.set_ylabel('u₂ (Standard Normal)', fontsize=11)
    ax.set_title('Most Probable Point in Standard Normal Space', fontsize=12, fontweight='bold')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')
    ax.set_xlim([-5, 2])
    ax.set_ylim([-4, 4])
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case3_single_loop_rbdo.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n案例3完成，结果已保存")
    
    return {
        'd_opt': d_opt,
        'u_opt': u_opt,
        'beta_opt': beta_opt
    }


def case4_sora_method():
    """
    案例4：SORA解耦法
    序列优化与可靠性评估
    """
    print("\n" + "=" * 70)
    print("案例4：SORA解耦法（序列优化与可靠性评估）")
    print("=" * 70)
    
    mu_X = np.array([50000.0, 4.0, 250e6])
    sigma_X = np.array([5000.0, 0.04, 25e6])
    beta_target = 3.0
    
    def stress_constraint(d, x):
        """应力约束函数"""
        b, h = d
        P, L, sigma_y = x
        sigma_max = 3 * P * L / (2 * b * h**2)
        return sigma_y - sigma_max
    
    def form_analysis(d, mu, sigma):
        """FORM分析"""
        def x2u(x):
            return (x - mu) / sigma
        
        def u2x(u):
            return mu + sigma * u
        
        def g_u(u):
            return stress_constraint(d, u2x(u))
        
        def grad_g_u(u):
            x = u2x(u)
            P, L, sigma_y = x
            b, h = d
            
            dg_dP = 3 * L / (2 * b * h**2)
            dg_dL = 3 * P / (2 * b * h**2)
            dg_dsigma_y = 1.0
            
            grad_x = np.array([dg_dP, dg_dL, dg_dsigma_y])
            return grad_x * sigma
        
        u = np.zeros(3)
        for _ in range(50):
            g_val = g_u(u)
            grad = grad_g_u(u)
            grad_norm_sq = np.dot(grad, grad)
            
            if grad_norm_sq < 1e-10:
                break
            
            u_new = (np.dot(grad, u) - g_val) / grad_norm_sq * grad
            
            if np.linalg.norm(u_new - u) < 1e-6:
                break
            
            u = u_new
        
        return u, np.linalg.norm(u)
    
    print("\n执行SORA方法...")
    print("交替进行确定性优化和可靠性分析")
    
    # SORA迭代
    d = np.array([0.2, 0.4])  # 初始设计
    mpp_history = []
    d_history = [d.copy()]
    beta_history = []
    
    print("\n  迭代历史:")
    print("  " + "-" * 60)
    print(f"  {'Cycle':<8} {'b (mm)':<10} {'h (mm)':<10} {'Area':<12} {'MPP Shift':<12} {'β':<8}")
    print("  " + "-" * 60)
    
    for cycle in range(10):
        # 步骤1：在当前设计点进行可靠性分析
        u_star, beta = form_analysis(d, mu_X, sigma_X)
        mpp = mu_X + sigma_X * u_star
        
        mpp_history.append(mpp.copy())
        beta_history.append(beta)
        
        # 计算移位向量
        shift = mpp - mu_X
        shift_norm = np.linalg.norm(shift)
        
        print(f"  {cycle:<8} {d[0]*1000:<10.2f} {d[1]*1000:<10.2f} "
              f"{d[0]*d[1]*1e6:<12.2f} {shift_norm:<12.4f} {beta:<8.4f}")
        
        # 检查收敛
        if cycle > 0 and shift_norm < 0.01:
            print(f"\n  收敛于第 {cycle + 1} 个循环")
            break
        
        # 步骤2：确定性优化，约束在MPP处满足
        # 简化：使用解析解
        # 约束：sigma_y(MPP) - 3*P(MPP)*L(MPP)/(2*b*h^2) = 0
        P_mpp, L_mpp, sigma_y_mpp = mpp
        target = 3 * P_mpp * L_mpp / (2 * sigma_y_mpp)
        
        # 最小化b*h，约束b*h^2 = target
        # 使用拉格朗日乘子法或参数扫描
        b_range = np.linspace(0.08, 0.5, 100)
        best_area = float('inf')
        best_b = d[0]
        best_h = d[1]
        
        for b in b_range:
            h = np.sqrt(target / b)
            if 0.1 <= h <= 1.0:
                area = b * h
                if area < best_area:
                    best_area = area
                    best_b = b
                    best_h = h
        
        d = np.array([best_b, best_h])
        d_history.append(d.copy())
    
    d_opt = d
    print(f"\nSORA最优解:")
    print(f"  宽度 b = {d_opt[0]:.4f} m = {d_opt[0]*1000:.2f} mm")
    print(f"  高度 h = {d_opt[1]:.4f} m = {d_opt[1]*1000:.2f} mm")
    print(f"  截面面积 A = {d_opt[0]*d_opt[1]*1e6:.2f} mm²")
    print(f"  最终可靠性指数 β = {beta_history[-1]:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 1. 设计变量收敛
    ax = axes[0]
    d_array = np.array(d_history) * 1000
    cycles = np.arange(len(d_history))
    
    ax.plot(cycles, d_array[:, 0], 'b-o', linewidth=2, markersize=8, label='Width b')
    ax.plot(cycles, d_array[:, 1], 'r-s', linewidth=2, markersize=8, label='Height h')
    
    ax.set_xlabel('SORA Cycle', fontsize=11)
    ax.set_ylabel('Dimension (mm)', fontsize=11)
    ax.set_title('Design Variable Convergence (SORA)', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 2. 可靠性指数收敛
    ax = axes[1]
    cycles_beta = np.arange(len(beta_history))
    
    ax.plot(cycles_beta, beta_history, 'g-^', linewidth=2, markersize=8)
    ax.axhline(y=beta_target, color='r', linestyle='--', linewidth=2, 
               label=f'Target β = {beta_target}')
    
    ax.set_xlabel('SORA Cycle', fontsize=11)
    ax.set_ylabel('Reliability Index β', fontsize=11)
    ax.set_title('Reliability Index Convergence', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case4_sora_method.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n案例4完成，结果已保存")
    
    return {
        'd_opt': d_opt,
        'beta_history': beta_history,
        'd_history': d_history
    }


def case5_robust_optimization():
    """
    案例5：鲁棒优化设计
    使用Taguchi方法和均值-方差模型
    """
    print("\n" + "=" * 70)
    print("案例5：鲁棒优化设计")
    print("=" * 70)
    
    # 问题：设计一个电路，使得输出电压对元件参数变化不敏感
    # 电路：分压电路 Vout = Vin * R2 / (R1 + R2)
    # 目标：Vout = 5V，同时最小化对R1和R2波动的敏感性
    
    Vin = 12.0  # 输入电压
    Vout_target = 5.0  # 目标输出电压
    
    # 电阻的变异系数
    cov_R = 0.05  # 5%
    
    def circuit_output(R1, R2, Vin_val):
        """电路输出电压"""
        return Vin_val * R2 / (R1 + R2)
    
    # 1. 名义优化（不考虑不确定性）
    print("\n1. 名义优化（不考虑不确定性）...")
    
    # 约束：Vout = Vout_target
    # R2 / (R1 + R2) = Vout_target / Vin
    # R2 = R1 * Vout_target / (Vin - Vout_target)
    ratio = Vout_target / (Vin - Vout_target)
    
    # 最小化总电阻 R1 + R2 = R1 * (1 + ratio)
    # 受限于 R1 >= 1kΩ, R2 >= 1kΩ
    R1_nom = 1000.0  # 最小值
    R2_nom = R1_nom * ratio
    
    Vout_nom = circuit_output(R1_nom, R2_nom, Vin)
    print(f"  名义设计: R1 = {R1_nom:.2f} Ω, R2 = {R2_nom:.2f} Ω")
    print(f"  输出电压: Vout = {Vout_nom:.4f} V")
    
    # 2. 鲁棒优化（最小化方差）
    print("\n2. 鲁棒优化（最小化输出电压方差）...")
    
    # 使用一阶近似计算方差
    # Vout = Vin * R2 / (R1 + R2)
    # dVout/dR1 = -Vin * R2 / (R1 + R2)^2
    # dVout/dR2 = Vin * R1 / (R1 + R2)^2
    
    def output_variance(R1, R2):
        """计算输出电压的方差"""
        dV_dR1 = -Vin * R2 / (R1 + R2)**2
        dV_dR2 = Vin * R1 / (R1 + R2)**2
        
        sigma_R1 = cov_R * R1
        sigma_R2 = cov_R * R2
        
        var = (dV_dR1 * sigma_R1)**2 + (dV_dR2 * sigma_R2)**2
        return var
    
    def output_mean(R1, R2):
        """输出电压均值（名义值）"""
        return circuit_output(R1, R2, Vin)
    
    # 鲁棒优化：最小化 w1*(Vout - Vtarget)^2 + w2*Var(Vout)
    def robust_objective(R1):
        R2 = R1 * ratio  # 保持比例以实现目标电压
        mean_error = (output_mean(R1, R2) - Vout_target)**2
        variance = output_variance(R1, R2)
        return 100 * mean_error + 10000 * variance  # 加权目标
    
    # 扫描寻找最优
    R1_range = np.linspace(1000, 50000, 500)
    obj_values = [robust_objective(R1) for R1 in R1_range]
    
    idx_opt = np.argmin(obj_values)
    R1_robust = R1_range[idx_opt]
    R2_robust = R1_robust * ratio
    
    print(f"  鲁棒设计: R1 = {R1_robust:.2f} Ω, R2 = {R2_robust:.2f} Ω")
    print(f"  输出电压均值: {output_mean(R1_robust, R2_robust):.4f} V")
    print(f"  输出电压标准差: {np.sqrt(output_variance(R1_robust, R2_robust)):.4f} V")
    
    # 3. 蒙特卡洛验证
    print("\n3. 蒙特卡洛验证...")
    np.random.seed(42)
    N_mc = 100000
    
    # 名义设计
    R1_samples_nom = np.random.normal(R1_nom, cov_R * R1_nom, N_mc)
    R2_samples_nom = np.random.normal(R2_nom, cov_R * R2_nom, N_mc)
    Vout_samples_nom = circuit_output(R1_samples_nom, R2_samples_nom, Vin)
    
    # 鲁棒设计
    R1_samples_rob = np.random.normal(R1_robust, cov_R * R1_robust, N_mc)
    R2_samples_rob = np.random.normal(R2_robust, cov_R * R2_robust, N_mc)
    Vout_samples_rob = circuit_output(R1_samples_rob, R2_samples_rob, Vin)
    
    print(f"\n  名义设计:")
    print(f"    均值 = {np.mean(Vout_samples_nom):.4f} V")
    print(f"    标准差 = {np.std(Vout_samples_nom):.4f} V")
    print(f"    信噪比 = {20*np.log10(np.mean(Vout_samples_nom)/np.std(Vout_samples_nom)):.2f} dB")
    
    print(f"\n  鲁棒设计:")
    print(f"    均值 = {np.mean(Vout_samples_rob):.4f} V")
    print(f"    标准差 = {np.std(Vout_samples_rob):.4f} V")
    print(f"    信噪比 = {20*np.log10(np.mean(Vout_samples_rob)/np.std(Vout_samples_rob)):.2f} dB")
    
    std_reduction = (np.std(Vout_samples_nom) - np.std(Vout_samples_rob)) / np.std(Vout_samples_nom) * 100
    print(f"\n  标准差降低: {std_reduction:.2f}%")
    
    # 4. Taguchi正交试验分析
    print("\n4. Taguchi正交试验分析...")
    
    # L9正交表（3水平，4因素）
    # 因素：R1, R2, Vin（噪声因素）
    # 这里简化分析R1和R2的影响
    
    levels = [0.9, 1.0, 1.1]  # 相对均值的比例
    
    snr_values = []
    for l1 in levels:
        for l2 in levels:
            R1_test = R1_robust * l1
            R2_test = R2_robust * l2
            Vout = circuit_output(R1_test, R2_test, Vin)
            snr_values.append((l1, l2, Vout))
    
    # 计算各水平的影响
    print(f"  正交试验结果（部分）:")
    print(f"  {'R1 Level':<12} {'R2 Level':<12} {'Vout (V)':<12}")
    print("  " + "-" * 36)
    for l1, l2, vout in snr_values[:9]:
        print(f"  {l1:<12.2f} {l2:<12.2f} {vout:<12.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 输出电压分布对比
    ax = axes[0, 0]
    bins = np.linspace(4.0, 6.0, 50)
    ax.hist(Vout_samples_nom, bins=bins, alpha=0.6, label='Nominal Design', 
            color='blue', density=True, edgecolor='white')
    ax.hist(Vout_samples_rob, bins=bins, alpha=0.6, label='Robust Design', 
            color='red', density=True, edgecolor='white')
    ax.axvline(x=Vout_target, color='g', linestyle='--', linewidth=2, 
               label=f'Target = {Vout_target} V')
    ax.set_xlabel('Output Voltage (V)', fontsize=11)
    ax.set_ylabel('Probability Density', fontsize=11)
    ax.set_title('Output Voltage Distribution Comparison', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 2. 方差随R1的变化
    ax = axes[0, 1]
    variances = [output_variance(R1, R1 * ratio) for R1 in R1_range]
    ax.semilogy(R1_range / 1000, variances, 'b-', linewidth=2)
    ax.axvline(x=R1_nom / 1000, color='g', linestyle='--', linewidth=2, 
               label=f'Nominal R1 = {R1_nom/1000:.1f} kΩ')
    ax.axvline(x=R1_robust / 1000, color='r', linestyle='--', linewidth=2, 
               label=f'Robust R1 = {R1_robust/1000:.1f} kΩ')
    ax.set_xlabel('R1 (kΩ)', fontsize=11)
    ax.set_ylabel('Output Variance (V²)', fontsize=11)
    ax.set_title('Output Variance vs R1', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 3. 鲁棒目标函数
    ax = axes[1, 0]
    ax.plot(R1_range / 1000, obj_values, 'purple', linewidth=2)
    ax.axvline(x=R1_robust / 1000, color='r', linestyle='--', linewidth=2, 
               label=f'Optimal R1 = {R1_robust/1000:.1f} kΩ')
    ax.set_xlabel('R1 (kΩ)', fontsize=11)
    ax.set_ylabel('Robust Objective', fontsize=11)
    ax.set_title('Robust Objective Function', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 4. 信噪比对比
    ax = axes[1, 1]
    designs = ['Nominal', 'Robust']
    snr_nom = 20 * np.log10(np.mean(Vout_samples_nom) / np.std(Vout_samples_nom))
    snr_rob = 20 * np.log10(np.mean(Vout_samples_rob) / np.std(Vout_samples_rob))
    snrs = [snr_nom, snr_rob]
    colors = ['#3498db', '#e74c3c']
    
    bars = ax.bar(designs, snrs, color=colors, edgecolor='black', linewidth=1.5)
    ax.set_ylabel('Signal-to-Noise Ratio (dB)', fontsize=11)
    ax.set_title('SNR Comparison (Higher is Better)', fontsize=12, fontweight='bold')
    ax.grid(True, axis='y', alpha=0.3)
    
    for bar, snr in zip(bars, snrs):
        height = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2., height,
                f'{snr:.2f} dB', ha='center', va='bottom', fontsize=11, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case5_robust_optimization.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n案例5完成，结果已保存")
    
    return {
        'R1_nom': R1_nom,
        'R2_nom': R2_nom,
        'R1_robust': R1_robust,
        'R2_robust': R2_robust,
        'std_reduction': std_reduction
    }


def case6_reliability_robust_optimization():
    """
    案例6：可靠性鲁棒综合优化
    结合可靠性和鲁棒性的多目标优化
    """
    print("\n" + "=" * 70)
    print("案例6：可靠性鲁棒综合优化")
    print("=" * 70)
    
    # 问题：设计一个弹簧-质量系统
    # 目标：1) 最小化成本（材料用量）
    #       2) 最大化可靠性（避免共振）
    #       3) 最大化鲁棒性（对参数变化不敏感）
    
    # 系统参数
    m = 10.0  # 质量 (kg)，固定
    
    # 随机变量
    # k: 弹簧刚度 ~ N(5000, 500) N/m
    # c: 阻尼系数 ~ N(50, 5) N·s/m
    # F: 激励幅值 ~ N(100, 10) N
    # omega: 激励频率 ~ N(15, 1.5) rad/s
    
    mu_params = np.array([5000.0, 50.0, 100.0, 15.0])
    sigma_params = np.array([500.0, 5.0, 10.0, 1.5])
    
    # 设计变量：弹簧刚度k（均值）
    # 我们需要选择k的标称值，使得系统既可靠又鲁棒
    
    def system_response(k, c, F, omega, m_val):
        """
        计算稳态响应幅值
        X = F / sqrt((k - m*omega^2)^2 + (c*omega)^2)
        """
        denom = np.sqrt((k - m_val * omega**2)**2 + (c * omega)**2)
        return F / denom
    
    def analyze_design(k_nom):
        """分析设计的性能、可靠性和鲁棒性"""
        # 1. 名义性能
        k = k_nom
        c, F, omega = mu_params[1:]
        X_nom = system_response(k, c, F, omega, m)
        
        # 2. 可靠性分析（避免过大振幅）
        X_limit = 0.05  # 最大允许振幅 50mm
        
        # 使用FORM近似
        def g_func(x):
            k_val, c_val, F_val, omega_val = x
            X = system_response(k_val, c_val, F_val, omega_val, m)
            return X_limit - X
        
        # 在设计点处线性化
        x0 = np.array([k, c, F, omega])
        g0 = g_func(x0)
        
        # 数值梯度
        eps = 1e-6
        grad = np.zeros(4)
        for i in range(4):
            x_plus = x0.copy()
            x_plus[i] += eps
            grad[i] = (g_func(x_plus) - g0) / eps
        
        # 可靠性指数（一阶近似）
        beta = g0 / np.sqrt(np.sum((grad * sigma_params)**2))
        Pf = stats.norm.cdf(-beta)
        
        # 3. 鲁棒性分析（响应方差）
        # 一阶泰勒展开
        var_X = np.sum((grad * sigma_params)**2)
        sigma_X = np.sqrt(var_X)
        
        # 信噪比
        snr = 20 * np.log10(X_nom / sigma_X) if sigma_X > 0 else 100
        
        return {
            'X_nom': X_nom,
            'beta': beta,
            'Pf': Pf,
            'sigma_X': sigma_X,
            'snr': snr
        }
    
    # 扫描设计空间
    print("\n扫描设计空间...")
    k_range = np.linspace(2000, 8000, 100)
    
    results = []
    for k in k_range:
        res = analyze_design(k)
        res['k'] = k
        results.append(res)
    
    # 提取数据
    X_noms = [r['X_nom'] for r in results]
    betas = [r['beta'] for r in results]
    Pfs = [r['Pf'] for r in results]
    sigma_Xs = [r['sigma_X'] for r in results]
    snrs = [r['snr'] for r in results]
    
    # 找到满足可靠性约束且鲁棒性最好的设计
    beta_target = 3.0
    feasible = [r for r in results if r['beta'] >= beta_target]
    
    if feasible:
        best_robust = min(feasible, key=lambda r: r['sigma_X'])
        k_opt = best_robust['k']
        print(f"\n最优设计:")
        print(f"  弹簧刚度 k = {k_opt:.2f} N/m")
        print(f"  名义响应幅值 = {best_robust['X_nom']*1000:.2f} mm")
        print(f"  可靠性指数 β = {best_robust['beta']:.4f}")
        print(f"  失效概率 Pf = {best_robust['Pf']:.4e}")
        print(f"  响应标准差 = {best_robust['sigma_X']*1000:.4f} mm")
        print(f"  信噪比 SNR = {best_robust['snr']:.2f} dB")
    
    # 多目标权衡分析
    print("\n多目标权衡分析:")
    
    # 加权目标：w1*cost + w2*(-beta) + w3*sigma_X
    # 假设成本与k成正比（材料用量）
    def multi_objective(k, w1, w2, w3):
        res = analyze_design(k)
        cost = k / 1000  # 归一化成本
        return w1 * cost + w2 * max(0, 3 - res['beta']) + w3 * res['sigma_X'] * 100
    
    # 不同权重组合
    weight_combinations = [
        (1.0, 0.0, 0.0),  # 仅成本
        (0.0, 1.0, 0.0),  # 仅可靠性
        (0.0, 0.0, 1.0),  # 仅鲁棒性
        (0.5, 0.3, 0.2),  # 平衡
        (0.3, 0.5, 0.2),  # 重可靠性
    ]
    
    print(f"\n  {'Weights (cost, rel, rob)':<30} {'Optimal k':<15} {'β':<10} {'σ_X (mm)':<12}")
    print("  " + "-" * 70)
    
    pareto_points = []
    for w1, w2, w3 in weight_combinations:
        obj_values = [multi_objective(k, w1, w2, w3) for k in k_range]
        k_opt_w = k_range[np.argmin(obj_values)]
        res_opt = analyze_design(k_opt_w)
        pareto_points.append((k_opt_w, res_opt['beta'], res_opt['sigma_X'], 
                             res_opt['X_nom'], w1, w2, w3))
        print(f"  ({w1:.1f}, {w2:.1f}, {w3:.1f}){'':<22} {k_opt_w:<15.2f} "
              f"{res_opt['beta']:<10.4f} {res_opt['sigma_X']*1000:<12.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 响应幅值vs弹簧刚度
    ax = axes[0, 0]
    ax.plot(k_range, np.array(X_noms) * 1000, 'b-', linewidth=2, label='Response Amplitude')
    ax.axhline(y=50, color='r', linestyle='--', linewidth=2, label='Limit = 50 mm')
    if feasible:
        ax.axvline(x=k_opt, color='g', linestyle='--', linewidth=2, 
                   label=f'Optimal k = {k_opt:.0f}')
    ax.set_xlabel('Spring Stiffness k (N/m)', fontsize=11)
    ax.set_ylabel('Response Amplitude (mm)', fontsize=11)
    ax.set_title('Steady-State Response vs Stiffness', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 2. 可靠性指数
    ax = axes[0, 1]
    ax.plot(k_range, betas, 'purple', linewidth=2)
    ax.axhline(y=beta_target, color='r', linestyle='--', linewidth=2, 
               label=f'Target β = {beta_target}')
    if feasible:
        ax.axvline(x=k_opt, color='g', linestyle='--', linewidth=2)
        ax.scatter([k_opt], [best_robust['beta']], c='green', s=150, zorder=5, 
                  marker='*', edgecolors='black', linewidths=1.5)
    ax.set_xlabel('Spring Stiffness k (N/m)', fontsize=11)
    ax.set_ylabel('Reliability Index β', fontsize=11)
    ax.set_title('Reliability Index vs Stiffness', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 3. 鲁棒性（标准差）
    ax = axes[1, 0]
    ax.plot(k_range, np.array(sigma_Xs) * 1000, 'orange', linewidth=2)
    if feasible:
        ax.axvline(x=k_opt, color='g', linestyle='--', linewidth=2, 
                   label=f'Optimal k = {k_opt:.0f}')
        ax.scatter([k_opt], [best_robust['sigma_X'] * 1000], c='green', s=150, zorder=5,
                  marker='*', edgecolors='black', linewidths=1.5)
    ax.set_xlabel('Spring Stiffness k (N/m)', fontsize=11)
    ax.set_ylabel('Response Std Dev (mm)', fontsize=11)
    ax.set_title('Robustness (Lower is Better)', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 4. Pareto前沿：可靠性vs鲁棒性
    ax = axes[1, 1]
    
    # 绘制所有设计的散点
    scatter = ax.scatter(betas, np.array(sigma_Xs) * 1000, c=k_range, cmap='viridis', 
                        alpha=0.6, s=30, edgecolors='none')
    plt.colorbar(scatter, ax=ax, label='Stiffness k (N/m)')
    
    # 标记Pareto点
    for k_p, beta_p, sigma_p, _, _, _, _ in pareto_points:
        ax.scatter([beta_p], [sigma_p * 1000], c='red', s=100, marker='o', 
                  edgecolors='black', linewidths=1.5, zorder=5)
    
    if feasible:
        ax.scatter([best_robust['beta']], [best_robust['sigma_X'] * 1000], c='green', 
                  s=200, marker='*', edgecolors='black', linewidths=2, 
                  label='Selected Design', zorder=6)
    
    ax.axvline(x=beta_target, color='r', linestyle='--', linewidth=2, alpha=0.5)
    ax.set_xlabel('Reliability Index β', fontsize=11)
    ax.set_ylabel('Response Std Dev (mm)', fontsize=11)
    ax.set_title('Pareto Front: Reliability vs Robustness', fontsize=12, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case6_reliability_robust_optimization.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n案例6完成，结果已保存")
    
    return {
        'k_opt': k_opt if feasible else None,
        'beta_opt': best_robust['beta'] if feasible else None,
        'pareto_points': pareto_points
    }


def main():
    """主函数：运行所有案例"""
    print("\n" + "=" * 70)
    print(" " * 15 + "主题067：可靠性优化与鲁棒设计")
    print("=" * 70)
    print("\n本程序演示可靠性优化与鲁棒设计的实现和应用")
    print("包括：FORM/SORM、双循环RBDO、单循环RBDO、SORA、鲁棒优化")
    
    results = {}
    
    results['Case1'] = case1_reliability_analysis()
    results['Case2'] = case2_double_loop_rbdo()
    results['Case3'] = case3_single_loop_rbdo()
    results['Case4'] = case4_sora_method()
    results['Case5'] = case5_robust_optimization()
    results['Case6'] = case6_reliability_robust_optimization()
    
    print("\n" + "=" * 70)
    print(" " * 20 + "所有案例运行完成")
    print("=" * 70)
    print("\n生成的可视化文件:")
    print("  1. output/case1_reliability_analysis.png - 可靠性分析与FORM方法")
    print("  2. output/case2_double_loop_rbdo.png - 双循环可靠性优化")
    print("  3. output/case3_single_loop_rbdo.png - 单循环可靠性优化")
    print("  4. output/case4_sora_method.png - SORA解耦法")
    print("  5. output/case5_robust_optimization.png - 鲁棒优化设计")
    print("  6. output/case6_reliability_robust_optimization.png - 可靠性鲁棒综合优化")
    
    print("\n" + "=" * 70)
    
    return results


if __name__ == "__main__":
    main()