多场耦合优化-主题036-鲁棒性优化设计
·
主题036：鲁棒性优化设计
一、引言
1.1 什么是鲁棒性优化
鲁棒性优化（Robust Optimization, RO）是一种在不确定性环境下寻求最优解的优化方法。与传统的确定性优化不同，鲁棒性优化考虑设计参数、环境条件或模型本身的不确定性，寻找对这些不确定性不敏感的最优设计。
核心思想：一个好的设计不仅要在名义条件下表现优异，还要在参数波动、制造误差、环境变化等不确定因素的影响下保持稳定性能。
1.2 鲁棒性优化的重要性
在实际工程应用中，不确定性无处不在：
- 制造误差：加工精度限制导致实际尺寸与设计值存在偏差
- 材料分散性：材料性能参数（弹性模量、强度等）存在批次差异
- 环境变化：工作温度、湿度、载荷等环境条件波动
- 模型不确定性：仿真模型简化带来的误差
- 老化退化：长期使用导致的性能衰减
传统优化的局限性：
- 确定性优化得到的"最优解"往往位于约束边界
- 微小的参数扰动可能导致约束违反或性能急剧下降
- 实际产品的合格率和可靠性难以保证
鲁棒性优化的优势：
- 提高产品质量的一致性
- 降低对制造精度的苛刻要求
- 增强产品适应环境变化的能力
- 提高整体可靠性和经济性
1.3 鲁棒性优化的发展历程
早期阶段（1970s-1980s）：
- Taguchi方法在工业界的广泛应用
- 强调减少产品质量的变异性
- 通过正交试验设计优化稳健性
理论发展阶段（1990s-2000s）：
- 鲁棒优化数学理论的建立
- Ben-Tal和Nemirovski的凸优化鲁棒方法
- 随机规划与鲁棒优化的结合
工程应用阶段（2000s至今）：
- 结构鲁棒性优化设计
- 多学科鲁棒优化
- 基于仿真的鲁棒性优化
- 智能优化算法在鲁棒优化中的应用
二、鲁棒性优化的数学基础
2.1 不确定性建模
2.1.1 不确定性的分类
按来源分类：
-
偶然不确定性（Aleatory Uncertainty）：
- 固有的随机性，不可消除
- 如材料性能的随机波动
- 通常用概率方法描述
-
认知不确定性（Epistemic Uncertainty）：
- 由于知识不足导致
- 如模型简化带来的误差
- 可用区间分析、模糊理论描述
按表现形式分类：
- 参数不确定性：设计参数或环境参数的波动
- 模型不确定性：数学模型与实际系统的差异
- 载荷不确定性：外部载荷的随机变化
2.1.2 不确定性描述方法
概率方法：
当不确定性具有统计规律时，用随机变量描述：
X ~ N(μ, σ²) # 正态分布
X ~ U(a, b) # 均匀分布
区间方法：
当只知道参数变化范围时：
x ∈ [x^L, x^U] # 区间不确定性
模糊方法：
当不确定性具有模糊性时，使用模糊集或可能性理论。
凸集方法：
用凸集描述不确定性集合：
U = {u : ||u|| ≤ ρ} # 范数有界不确定性
2.2 鲁棒性优化的数学模型
2.2.1 一般形式
考虑不确定性的优化问题可表示为：
min f(x, ξ)
s.t. g_j(x, ξ) ≤ 0, j = 1, ..., m
x ∈ X
其中：
- x：设计变量
- ξ：不确定性参数
- f：目标函数（可能依赖于不确定性）
- g_j：约束函数
2.2.2 鲁棒性优化的分类
1. 鲁棒可行性问题：
寻找对所有不确定性都可行的解：
min f(x)
s.t. g_j(x, ξ) ≤ 0, ∀ξ ∈ U, j = 1, ..., m
2. 鲁棒最优性问题：
在最坏情况下优化目标：
min max_{ξ ∈ U} f(x, ξ)
s.t. g_j(x, ξ) ≤ 0, ∀ξ ∈ U, j = 1, ..., m
3. 可调节鲁棒优化：
考虑决策可以分阶段进行：
min c^T x + max_{ξ ∈ U} min_{y(ξ)} d^T y(ξ)
s.t. Ax + By(ξ) ≤ b, ∀ξ ∈ U
2.3 鲁棒性度量指标
2.3.1 性能鲁棒性
衡量目标函数对不确定性的敏感程度：
均值性能：
μ_f = E[f(x, ξ)]
性能方差：
σ_f² = Var[f(x, ξ)]
信噪比（Taguchi）：
SN = 10 log(μ_f² / σ_f²)
2.3.2 可行性鲁棒性
衡量约束在不确定性下的满足程度：
约束满足概率：
P(g_j(x, ξ) ≤ 0) ≥ R_j
最坏情况约束裕度：
Δg_j = min_{ξ ∈ U} (-g_j(x, ξ))
2.3.3 灵敏度指标
一阶灵敏度：
S_i = ∂f/∂ξ_i
二阶灵敏度（曲率）：
H_ij = ∂²f/∂ξ_i∂ξ_j
三、鲁棒性优化方法
3.1 基于梯度的鲁棒性优化
3.1.1 一阶泰勒展开法
假设不确定性较小，对目标函数进行一阶泰勒展开：
f(x, ξ) ≈ f(x, μ_ξ) + ∇_ξ f(x, μ_ξ)^T (ξ - μ_ξ)
均值近似：
μ_f ≈ f(x, μ_ξ)
方差近似：
σ_f² ≈ ∇_ξ f(x, μ_ξ)^T Σ_ξ ∇_ξ f(x, μ_ξ)
其中 Σ_ξ 是 ξ 的协方差矩阵。
优化问题：
min f(x, μ_ξ) + w·σ_f
s.t. g_j(x, μ_ξ) + k·σ_{g_j} ≤ 0
其中 w 和 k 是权重系数。
3.1.2 二阶泰勒展开法
对于非线性较强的问题，使用二阶展开：
f(x, ξ) ≈ f(x, μ_ξ) + ∇_ξ f^T (ξ - μ_ξ) + 1/2 (ξ - μ_ξ)^T H (ξ - μ_ξ)
均值修正：
μ_f ≈ f(x, μ_ξ) + 1/2 tr(H Σ_ξ)
方差修正：
σ_f² ≈ ∇_ξ f^T Σ_ξ ∇_ξ f + 1/2 tr(H Σ_ξ H Σ_ξ)
3.2 基于抽样的鲁棒性优化
3.2.1 蒙特卡洛模拟法
通过大量随机抽样估计统计特性：
μ_f ≈ (1/N) Σ_{i=1}^N f(x, ξ_i)
σ_f² ≈ (1/(N-1)) Σ_{i=1}^N (f(x, ξ_i) - μ_f)²
优点：
- 适用于任意复杂的非线性问题
- 可以估计完整的概率分布
缺点：
- 计算成本高
- 需要大量样本保证精度
3.2.2 拉丁超立方抽样（LHS）
分层抽样方法，提高抽样效率：
1. 将每个变量的取值范围分成N个等概率区间
2. 在每个区间随机抽取一个样本
3. 组合各变量的样本，形成N个样本点
优点：
- 样本覆盖更均匀
- 方差估计更稳定
- 所需样本数少于纯随机抽样
3.3 最坏情况优化（Worst-Case Optimization）
3.3.1 基本思想
在最坏的不确定性情况下优化性能：
min max_{ξ ∈ U} f(x, ξ)
s.t. max_{ξ ∈ U} g_j(x, ξ) ≤ 0, j = 1, ..., m
3.3.2 区间不确定性处理
对于区间不确定性 ξ ∈ [ξ^L, ξ^U]，最坏情况分析：
线性函数：
max f(x, ξ) = f(x, ξ^U) (若系数为正)
非线性函数：
需要求解内部最大化问题：
ξ* = argmax_{ξ ∈ U} f(x, ξ)
3.3.3 保守性控制
最坏情况优化往往过于保守，可通过以下方式调节：
预算不确定性：
U(Γ) = {ξ : Σ_i |ξ_i - μ_i|/σ_i ≤ Γ}
其中 Γ 是不确定性预算，控制保守程度。
3.4 基于Taguchi方法的鲁棒性设计
3.4.1 Taguchi质量损失函数
质量损失与偏离目标值的平方成正比：
L(y) = k(y - m)²
其中：
- y：实际性能值
- m：目标值
- k：损失系数
3.4.2 三阶段设计
1. 系统设计：
确定基本设计方案和结构。
2. 参数设计：
优化设计参数，使性能对噪声不敏感。
3. 容差设计：
确定各参数的允许偏差范围。
3.4.3 正交试验设计
使用正交表安排试验，高效评估因素效应：
L_9(3^4)：9次试验，4个因素，每个因素3个水平
L_16(4^5)：16次试验，5个因素，每个因素4个水平
3.5 随机规划方法
3.5.1 两阶段随机规划
决策分为两个阶段：
第一阶段（在不确定性实现前）：
做出"这里-and-now"决策 x。
第二阶段（在不确定性实现后）：
根据实际实现 ξ，做出"等待-and-see"决策 y(ξ)。
数学模型：
min c^T x + E[Q(x, ξ)]
s.t. Ax = b, x ≥ 0
其中 Q(x, ξ) = min_y {q(ξ)^T y : Wy = h(ξ) - T(ξ)x, y ≥ 0}
3.5.2 机会约束规划
允许约束以一定概率不满足：
min f(x)
s.t. P(g_j(x, ξ) ≤ 0) ≥ 1 - α_j, j = 1, ..., m
其中 α_j 是允许的风险水平。
四、鲁棒性优化的求解策略
4.1 双层优化方法
鲁棒性优化常可表述为双层优化问题：
上层：优化设计变量 x
下层：在最坏情况下评估性能
min_x f(x, ξ*(x))
s.t. g_j(x, ξ*(x)) ≤ 0
其中 ξ*(x) = argmax_{ξ ∈ U} f(x, ξ)
求解策略：
- 交替求解上下层问题
- 使用KKT条件将下层问题转化为约束
- 近似方法（如泰勒展开）
4.2 序列优化方法
4.2.1 序列鲁棒优化（SRO）
迭代 k：
1. 在当前设计点 x_k 进行不确定性分析
2. 构建局部鲁棒性近似模型
3. 求解近似鲁棒优化问题，得到 x_{k+1}
4. 检查收敛性
4.2.2 信赖域方法
在信赖域内建立代理模型：
min f̃_k(x) + w·σ̃_f(x)
s.t. ||x - x_k|| ≤ Δ_k
g̃_j(x) + k·σ̃_{g_j}(x) ≤ 0
其中 f̃_k 是局部代理模型，Δ_k 是信赖域半径。
4.3 智能优化算法
4.3.1 遗传算法（GA）
适应度函数设计：
Fitness(x) = -[f(x, μ_ξ) + w·σ_f(x)] - P·max(0, g_max(x))
特殊算子：
- 不确定性抽样评估个体鲁棒性
- 多目标处理（均值和方差）
4.3.2 粒子群优化（PSO）
速度更新考虑不确定性：
v_i^{k+1} = w·v_i^k + c1·r1·(pbest_i - x_i^k) + c2·r2·(gbest - x_i^k) + c3·r3·∇_x σ_f
4.3.3 多目标进化算法
将鲁棒性优化转化为多目标问题：
min [f(x, μ_ξ), σ_f(x)]
s.t. μ_{g_j}(x) + k·σ_{g_j}(x) ≤ 0
使用NSGA-II、MOEA/D等算法求解Pareto前沿。
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
鲁棒性优化设计仿真程序
主题036：多场耦合优化仿真
本程序包含4个鲁棒性优化仿真案例：
1. 悬臂梁鲁棒性优化设计（泰勒展开法）
2. 多目标鲁棒性优化（Pareto前沿）
3. 最坏情况优化与保守性分析
4. 基于抽样的鲁棒性优化对比
作者：仿真专家
日期：2026-02-28
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle, Ellipse
from matplotlib.collections import PatchCollection
import matplotlib.animation as animation
from scipy.optimize import minimize, differential_evolution, NonlinearConstraint
from scipy.stats import norm, uniform
import os
import warnings
from PIL import Image
warnings.filterwarnings('ignore')
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150
# 创建输出目录
output_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'output')
os.makedirs(output_dir, exist_ok=True)
print("=" * 80)
print("鲁棒性优化设计仿真")
print("主题036：多场耦合优化仿真")
print("=" * 80)
# =============================================================================
# 案例1：悬臂梁鲁棒性优化设计（泰勒展开法）
# =============================================================================
print("\n" + "=" * 80)
print("案例1：悬臂梁鲁棒性优化设计")
print("=" * 80)
def case1_cantilever_robust():
"""悬臂梁鲁棒性优化 - 基于泰勒展开法"""
print("\n【问题描述】")
print("设计矩形截面悬臂梁，在考虑材料属性和载荷不确定性的情况下，")
print("最小化重量同时保证挠度约束的鲁棒性满足。")
# 参数设置
L = 2.0 # 梁长度 (m)
rho = 7850 # 钢材密度 (kg/m³)
# 不确定性参数
E_mean = 200e9 # 弹性模量均值 (Pa)
sigma_E = 20e9 # 弹性模量标准差 (Pa, COV=0.1)
P_mean = 10e3 # 载荷均值 (N)
sigma_P = 2e3 # 载荷标准差 (N, COV=0.2)
delta_allow = 0.01 # 许用挠度 (m)
# 设计变量: x = [b, h] (宽度和高度)
print(f"\n【参数设置】")
print(f" 梁长度 L = {L} m")
print(f" 密度 ρ = {rho} kg/m³")
print(f" 弹性模量 E ~ N({E_mean/1e9:.0f}, {sigma_E/1e9:.1f}²) GPa")
print(f" 载荷 P ~ N({P_mean/1e3:.0f}, {sigma_P/1e3:.0f}²) kN")
print(f" 许用挠度 δ_allow = {delta_allow*1000:.1f} mm")
# 确定性优化（对比基准）
print(f"\n【确定性优化】")
def objective_det(x):
b, h = x
return rho * L * b * h # 重量
def constraint_det(x):
b, h = x
I = b * h**3 / 12 # 惯性矩
delta = P_mean * L**3 / (3 * E_mean * I)
return delta_allow - delta
result_det = minimize(
objective_det,
[0.05, 0.1],
method='SLSQP',
bounds=[(0.01, 0.5), (0.02, 0.5)],
constraints={'type': 'ineq', 'fun': constraint_det}
)
b_det, h_det = result_det.x
weight_det = result_det.fun
I_det = b_det * h_det**3 / 12
delta_det = P_mean * L**3 / (3 * E_mean * I_det)
print(f" 最优设计: b={b_det*1000:.2f}mm, h={h_det*1000:.2f}mm")
print(f" 重量: {weight_det:.2f} kg")
print(f" 挠度: {delta_det*1000:.2f} mm")
# 鲁棒性优化 - 一阶泰勒展开法
print(f"\n【鲁棒性优化（一阶泰勒展开法）】")
def compute_robust_metrics(x):
"""计算鲁棒性指标（均值和标准差）"""
b, h = x
I = b * h**3 / 12
# 名义挠度
delta_nom = P_mean * L**3 / (3 * E_mean * I)
# 对E和P的偏导数
d_delta_dE = -P_mean * L**3 / (3 * E_mean**2 * I)
d_delta_dP = L**3 / (3 * E_mean * I)
# 挠度的标准差（一阶近似）
sigma_delta = np.sqrt((d_delta_dE * sigma_E)**2 + (d_delta_dP * sigma_P)**2)
return delta_nom, sigma_delta
def objective_robust(x, w_sigma=1000):
"""鲁棒目标函数：重量 + 权重×挠度标准差"""
b, h = x
weight = rho * L * b * h
_, sigma_delta = compute_robust_metrics(x)
return weight + w_sigma * sigma_delta
def constraint_robust(x, k=3):
"""鲁棒约束：均值 + k×标准差 ≤ 许用值"""
delta_mean, sigma_delta = compute_robust_metrics(x)
return delta_allow - (delta_mean + k * sigma_delta)
# 不同权重下的鲁棒优化
weights = [0, 500, 1000, 2000]
results_robust = []
for w in weights:
result = minimize(
lambda x: objective_robust(x, w),
[0.05, 0.1],
method='SLSQP',
bounds=[(0.01, 0.5), (0.02, 0.5)],
constraints={'type': 'ineq', 'fun': constraint_robust}
)
b, h = result.x
weight = rho * L * b * h
delta_mean, sigma_delta = compute_robust_metrics(result.x)
results_robust.append({
'w': w,
'b': b,
'h': h,
'weight': weight,
'delta_mean': delta_mean,
'sigma_delta': sigma_delta,
'success': result.success
})
print(f" 权重w={w:4d}: b={b*1000:.2f}mm, h={h*1000:.2f}mm, "
f"W={weight:.2f}kg, δ_mean={delta_mean*1000:.2f}mm, σ_δ={sigma_delta*1000:.3f}mm")
# 蒙特卡洛验证
print(f"\n【蒙特卡洛验证（n=10000）】")
np.random.seed(42)
n_mc = 10000
E_samples = np.random.normal(E_mean, sigma_E, n_mc)
P_samples = np.random.normal(P_mean, sigma_P, n_mc)
for i, res in enumerate([results_robust[0], results_robust[-1]]):
b, h = res['b'], res['h']
I = b * h**3 / 12
delta_samples = P_samples * L**3 / (3 * E_samples * I)
delta_mean_mc = np.mean(delta_samples)
delta_std_mc = np.std(delta_samples)
violation_rate = np.sum(delta_samples > delta_allow) / n_mc * 100
label = "确定性设计" if i == 0 else "鲁棒设计"
print(f" {label}:")
print(f" 蒙特卡洛均值: {delta_mean_mc*1000:.2f} mm")
print(f" 蒙特卡洛标准差: {delta_std_mc*1000:.3f} mm")
print(f" 约束违反率: {violation_rate:.2f}%")
# ========== 可视化 ==========
fig = plt.figure(figsize=(16, 12))
# 子图1: 悬臂梁示意图
ax1 = fig.add_subplot(2, 3, 1)
# 绘制梁（确定性设计）
y_beam = np.linspace(0, L, 50)
x_beam = np.zeros_like(y_beam)
ax1.fill_betweenx(y_beam, x_beam - b_det/2, x_beam + b_det/2,
alpha=0.3, color='blue', label='确定性设计')
# 绘制梁（鲁棒设计）
b_rob = results_robust[-1]['b']
h_rob = results_robust[-1]['h']
ax1.fill_betweenx(y_beam, x_beam - b_rob/2, x_beam + b_rob/2,
alpha=0.3, color='red', label='鲁棒设计')
# 固定端
ax1.plot([-0.05, 0.05], [0, 0], 'k-', linewidth=3)
ax1.plot([-0.03, -0.03], [-0.05, 0], 'k-', linewidth=2)
ax1.plot([0.03, 0.03], [-0.05, 0], 'k-', linewidth=2)
# 载荷箭头
ax1.arrow(0, L*0.9, 0.1, 0, head_width=0.08, head_length=0.03,
fc='green', ec='green', linewidth=2)
ax1.text(0.15, L*0.9, 'P', fontsize=12, color='green')
ax1.set_xlim(-0.15, 0.25)
ax1.set_ylim(-0.1, L*1.1)
ax1.set_aspect('equal')
ax1.set_xlabel('宽度方向 (m)', fontsize=11)
ax1.set_ylabel('长度方向 (m)', fontsize=11)
ax1.set_title('悬臂梁设计对比', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
# 子图2: 重量-鲁棒性权衡
ax2 = fig.add_subplot(2, 3, 2)
weights_plot = [r['w'] for r in results_robust]
weights_val = [r['weight'] for r in results_robust]
sigmas = [r['sigma_delta']*1000 for r in results_robust]
ax2_twin = ax2.twinx()
line1 = ax2.plot(weights_plot, weights_val, 'b-o', linewidth=2, markersize=8, label='重量')
line2 = ax2_twin.plot(weights_plot, sigmas, 'r-s', linewidth=2, markersize=8, label='挠度标准差')
ax2.set_xlabel('鲁棒性权重 w', fontsize=11)
ax2.set_ylabel('重量 (kg)', fontsize=11, color='blue')
ax2_twin.set_ylabel('挠度标准差 (mm)', fontsize=11, color='red')
ax2.tick_params(axis='y', labelcolor='blue')
ax2_twin.tick_params(axis='y', labelcolor='red')
ax2.set_title('重量-鲁棒性权衡曲线', fontsize=13, fontweight='bold')
ax2.grid(True, alpha=0.3)
lines = line1 + line2
labels = [l.get_label() for l in lines]
ax2.legend(lines, labels, fontsize=10)
# 子图3: 截面尺寸演化
ax3 = fig.add_subplot(2, 3, 3)
bs = [r['b']*1000 for r in results_robust]
hs = [r['h']*1000 for r in results_robust]
x_pos = np.arange(len(weights_plot))
width = 0.35
bars1 = ax3.bar(x_pos - width/2, bs, width, label='宽度 b', color='steelblue', edgecolor='navy')
bars2 = ax3.bar(x_pos + width/2, hs, width, label='高度 h', color='coral', edgecolor='darkred')
ax3.set_xlabel('鲁棒性权重', fontsize=11)
ax3.set_ylabel('尺寸 (mm)', fontsize=11)
ax3.set_title('截面尺寸随权重变化', fontsize=13, fontweight='bold')
ax3.set_xticks(x_pos)
ax3.set_xticklabels([f'{w}' for w in weights_plot])
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, axis='y')
# 子图4: 挠度分布对比（蒙特卡洛）
ax4 = fig.add_subplot(2, 3, 4)
# 确定性设计的挠度分布
I_det = b_det * h_det**3 / 12
delta_samples_det = P_samples * L**3 / (3 * E_samples * I_det)
# 鲁棒设计的挠度分布
I_rob = b_rob * h_rob**3 / 12
delta_samples_rob = P_samples * L**3 / (3 * E_samples * I_rob)
ax4.hist(delta_samples_det*1000, bins=50, alpha=0.5, color='blue',
label='确定性设计', density=True)
ax4.hist(delta_samples_rob*1000, bins=50, alpha=0.5, color='red',
label='鲁棒设计', density=True)
ax4.axvline(x=delta_allow*1000, color='k', linestyle='--', linewidth=2,
label=f'许用值={delta_allow*1000:.1f}mm')
ax4.set_xlabel('挠度 (mm)', fontsize=11)
ax4.set_ylabel('概率密度', fontsize=11)
ax4.set_title('挠度分布对比（蒙特卡洛）', fontsize=13, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
# 子图5: 灵敏度分析
ax5 = fig.add_subplot(2, 3, 5)
# 计算不同设计点对E和P的灵敏度
designs = ['确定性', '鲁棒(w=500)', '鲁棒(w=2000)']
indices = [0, 1, 3]
sens_E = []
sens_P = []
for idx in indices:
b, h = results_robust[idx]['b'], results_robust[idx]['h']
I = b * h**3 / 12
d_delta_dE = -P_mean * L**3 / (3 * E_mean**2 * I)
d_delta_dP = L**3 / (3 * E_mean * I)
sens_E.append(abs(d_delta_dE) * sigma_E)
sens_P.append(abs(d_delta_dP) * sigma_P)
x = np.arange(len(designs))
width = 0.35
bars1 = ax5.bar(x - width/2, sens_E, width, label='对E的灵敏度', color='skyblue', edgecolor='navy')
bars2 = ax5.bar(x + width/2, sens_P, width, label='对P的灵敏度', color='lightcoral', edgecolor='darkred')
ax5.set_ylabel('灵敏度系数', fontsize=11)
ax5.set_title('设计对不确定性的灵敏度', fontsize=13, fontweight='bold')
ax5.set_xticks(x)
ax5.set_xticklabels(designs, fontsize=10)
ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3, axis='y')
# 子图6: 约束满足概率
ax6 = fig.add_subplot(2, 3, 6)
violation_rates = []
for res in results_robust:
b, h = res['b'], res['h']
I = b * h**3 / 12
delta_samples = P_samples * L**3 / (3 * E_samples * I)
violation_rate = np.sum(delta_samples > delta_allow) / n_mc * 100
violation_rates.append(violation_rate)
colors = ['green' if v < 1 else 'orange' if v < 5 else 'red' for v in violation_rates]
bars = ax6.bar(range(len(weights_plot)), violation_rates, color=colors, edgecolor='black')
ax6.axhline(y=1, color='r', linestyle='--', linewidth=2, label='1%阈值')
ax6.set_xlabel('鲁棒性权重', fontsize=11)
ax6.set_ylabel('约束违反率 (%)', fontsize=11)
ax6.set_title('约束满足概率分析', fontsize=13, fontweight='bold')
ax6.set_xticks(range(len(weights_plot)))
ax6.set_xticklabels([f'{w}' for w in weights_plot])
ax6.legend(fontsize=10)
ax6.grid(True, alpha=0.3, axis='y')
# 添加数值标签
for i, (bar, rate) in enumerate(zip(bars, violation_rates)):
ax6.text(bar.get_x() + bar.get_width()/2., bar.get_height() + 0.1,
f'{rate:.2f}%', ha='center', fontsize=9)
plt.suptitle('案例1：悬臂梁鲁棒性优化设计', fontsize=14, fontweight='bold', y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.96])
output_file = os.path.join(output_dir, 'case1_cantilever_robust.png')
plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
print(f"\n✓ 图像已保存: {output_file}")
plt.close()
return results_robust
results_case1 = case1_cantilever_robust()
# =============================================================================
# 案例2：多目标鲁棒性优化（Pareto前沿）
# =============================================================================
print("\n" + "=" * 80)
print("案例2：多目标鲁棒性优化")
print("=" * 80)
def case2_multiobjective_robust():
"""多目标鲁棒性优化 - 生成Pareto前沿"""
print("\n【问题描述】")
print("通过多目标优化方法，同时优化性能和鲁棒性，")
print("生成Pareto前沿，展示性能与鲁棒性的权衡关系。")
# 问题设置：一个简单的数学优化问题
# 目标1: f1 = x1^2 + x2^2 (性能，越小越好)
# 目标2: f2 = 不确定性引起的方差 (鲁棒性，越小越好)
# 约束: g = x1 + x2 >= 1
# 不确定性：目标函数系数存在不确定性
# f_actual = (c1 + δc1)*x1^2 + (c2 + δc2)*x2^2
# 其中 δc1 ~ N(0, σc1^2), δc2 ~ N(0, σc2^2)
c1, c2 = 1.0, 1.0
sigma_c1, sigma_c2 = 0.2, 0.2
print(f"\n【参数设置】")
print(f" 名义系数: c1={c1}, c2={c2}")
print(f" 不确定性: σc1={sigma_c1}, σc2={sigma_c2}")
print(f" 约束: x1 + x2 >= 1")
# 使用加权和方法生成Pareto前沿
print(f"\n【生成Pareto前沿】")
n_points = 50
weights_f1 = np.linspace(0, 1, n_points)
pareto_front = []
for w1 in weights_f1:
w2 = 1 - w1
def objective(x):
x1, x2 = x
# 名义性能
f1 = c1 * x1**2 + c2 * x2**2
# 方差（鲁棒性度量）- 一阶近似
# Var(f) ≈ (∂f/∂c1)^2 * σc1^2 + (∂f/∂c2)^2 * σc2^2
# ∂f/∂c1 = x1^2, ∂f/∂c2 = x2^2
f2 = (x1**2)**2 * sigma_c1**2 + (x2**2)**2 * sigma_c2**2
return w1 * f1 + w2 * f2
def constraint(x):
return x[0] + x[1] - 1 # >= 0
result = minimize(
objective,
[0.5, 0.5],
method='SLSQP',
bounds=[(0, 2), (0, 2)],
constraints={'type': 'ineq', 'fun': constraint}
)
if result.success:
x1, x2 = result.x
f1 = c1 * x1**2 + c2 * x2**2
f2 = (x1**2)**2 * sigma_c1**2 + (x2**2)**2 * sigma_c2**2
pareto_front.append({'x': [x1, x2], 'f1': f1, 'f2': f2, 'w1': w1})
print(f" Pareto前沿点数: {len(pareto_front)}")
# 提取Pareto前沿数据
f1_values = [p['f1'] for p in pareto_front]
f2_values = [p['f2'] for p in pareto_front]
# 蒙特卡洛验证
print(f"\n【蒙特卡洛验证】")
np.random.seed(42)
n_mc = 5000
# 选择几个代表性的Pareto解进行验证
indices = [0, len(pareto_front)//4, len(pareto_front)//2, 3*len(pareto_front)//4, len(pareto_front)-1]
mc_results = []
for idx in indices:
p = pareto_front[idx]
x1, x2 = p['x']
# 生成不确定性样本
c1_samples = np.random.normal(c1, sigma_c1, n_mc)
c2_samples = np.random.normal(c2, sigma_c2, n_mc)
# 计算实际目标函数值
f_samples = c1_samples * x1**2 + c2_samples * x2**2
mc_results.append({
'idx': idx,
'x': [x1, x2],
'f1_nominal': p['f1'],
'f1_mean_mc': np.mean(f_samples),
'f1_std_mc': np.std(f_samples),
'f2_analytical': p['f2']
})
for r in mc_results:
print(f" 解 {r['idx']}: x=({r['x'][0]:.3f}, {r['x'][1]:.3f}), "
f"f1_nominal={r['f1_nominal']:.4f}, f1_mean_mc={r['f1_mean_mc']:.4f}, "
f"f1_std={r['f1_std_mc']:.4f}")
# ========== 可视化 ==========
fig = plt.figure(figsize=(16, 10))
# 子图1: Pareto前沿（目标空间）
ax1 = fig.add_subplot(2, 3, 1)
scatter = ax1.scatter(f1_values, f2_values, c=weights_f1, cmap='viridis',
s=50, edgecolors='black', linewidth=0.5)
plt.colorbar(scatter, ax=ax1, label='性能权重 w₁')
# 标记极端点
ax1.scatter(f1_values[0], f2_values[0], c='red', s=200, marker='*',
label='最鲁棒', zorder=5)
ax1.scatter(f1_values[-1], f2_values[-1], c='blue', s=200, marker='*',
label='最优性能', zorder=5)
ax1.set_xlabel('性能目标 f₁', fontsize=11)
ax1.set_ylabel('方差目标 f₂', fontsize=11)
ax1.set_title('Pareto前沿（目标空间）', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
# 子图2: 设计空间中的Pareto解
ax2 = fig.add_subplot(2, 3, 2)
x1_values = [p['x'][0] for p in pareto_front]
x2_values = [p['x'][1] for p in pareto_front]
scatter2 = ax2.scatter(x1_values, x2_values, c=weights_f1, cmap='viridis',
s=50, edgecolors='black', linewidth=0.5)
plt.colorbar(scatter2, ax=ax2, label='性能权重 w₁')
# 绘制约束边界
x1_boundary = np.linspace(0, 2, 100)
x2_boundary = 1 - x1_boundary
ax2.plot(x1_boundary, x2_boundary, 'r--', linewidth=2, label='约束边界')
ax2.fill_between(x1_boundary, 0, x2_boundary, alpha=0.2, color='red', label='不可行域')
ax2.set_xlabel('x₁', fontsize=11)
ax2.set_ylabel('x₂', fontsize=11)
ax2.set_title('Pareto解（设计空间）', fontsize=13, fontweight='bold')
ax2.set_xlim(0, 2)
ax2.set_ylim(0, 2)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
# 子图3: 权重-目标值关系
ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(weights_f1, f1_values, 'b-', linewidth=2, label='性能 f₁')
ax3_twin = ax3.twinx()
ax3_twin.plot(weights_f1, f2_values, 'r-', linewidth=2, label='方差 f₂')
ax3.set_xlabel('性能权重 w₁', fontsize=11)
ax3.set_ylabel('性能 f₁', fontsize=11, color='blue')
ax3_twin.set_ylabel('方差 f₂', fontsize=11, color='red')
ax3.tick_params(axis='y', labelcolor='blue')
ax3_twin.tick_params(axis='y', labelcolor='red')
ax3.set_title('权重与目标值关系', fontsize=13, fontweight='bold')
ax3.grid(True, alpha=0.3)
# 子图4: 蒙特卡洛验证对比
ax4 = fig.add_subplot(2, 3, 4)
indices_plot = [r['idx'] for r in mc_results]
f1_nominal_plot = [r['f1_nominal'] for r in mc_results]
f1_mean_mc_plot = [r['f1_mean_mc'] for r in mc_results]
x = np.arange(len(indices_plot))
width = 0.35
bars1 = ax4.bar(x - width/2, f1_nominal_plot, width, label='名义值',
color='steelblue', edgecolor='navy')
bars2 = ax4.bar(x + width/2, f1_mean_mc_plot, width, label='蒙特卡洛均值',
color='coral', edgecolor='darkred')
ax4.set_xlabel('Pareto解编号', fontsize=11)
ax4.set_ylabel('性能 f₁', fontsize=11)
ax4.set_title('名义值与蒙特卡洛验证对比', fontsize=13, fontweight='bold')
ax4.set_xticks(x)
ax4.set_xticklabels([f'{i}' for i in indices_plot])
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3, axis='y')
# 子图5: 不确定性影响可视化
ax5 = fig.add_subplot(2, 3, 5)
# 选择两个极端解进行比较
p_robust = pareto_front[0] # 最鲁棒
p_performance = pareto_front[-1] # 最优性能
# 生成大量样本
c1_samples = np.random.normal(c1, sigma_c1, 10000)
c2_samples = np.random.normal(c2, sigma_c2, 10000)
f_robust_samples = c1_samples * p_robust['x'][0]**2 + c2_samples * p_robust['x'][1]**2
f_perf_samples = c1_samples * p_performance['x'][0]**2 + c2_samples * p_performance['x'][1]**2
ax5.hist(f_robust_samples, bins=50, alpha=0.5, color='green',
label=f'鲁棒解 (σ={np.std(f_robust_samples):.3f})', density=True)
ax5.hist(f_perf_samples, bins=50, alpha=0.5, color='orange',
label=f'性能解 (σ={np.std(f_perf_samples):.3f})', density=True)
ax5.axvline(x=p_robust['f1'], color='green', linestyle='--', linewidth=2)
ax5.axvline(x=p_performance['f1'], color='orange', linestyle='--', linewidth=2)
ax5.set_xlabel('实际性能值', fontsize=11)
ax5.set_ylabel('概率密度', fontsize=11)
ax5.set_title('不确定性下的性能分布', fontsize=13, fontweight='bold')
ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3)
# 子图6: 鲁棒性指标对比
ax6 = fig.add_subplot(2, 3, 6)
f1_stds = [r['f1_std_mc'] for r in mc_results]
f2_analytical = [r['f2_analytical'] for r in mc_results]
ax6.scatter(f2_analytical, f1_stds, c='blue', s=100, edgecolors='black', linewidth=2)
# 绘制y=x参考线
max_val = max(max(f2_analytical), max(f1_stds))
ax6.plot([0, max_val], [0, max_val], 'r--', linewidth=2, label='y=x')
ax6.set_xlabel('解析方差 f₂', fontsize=11)
ax6.set_ylabel('蒙特卡洛标准差', fontsize=11)
ax6.set_title('鲁棒性指标验证', fontsize=13, fontweight='bold')
ax6.legend(fontsize=10)
ax6.grid(True, alpha=0.3)
plt.suptitle('案例2：多目标鲁棒性优化', fontsize=14, fontweight='bold', y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.96])
output_file = os.path.join(output_dir, 'case2_multiobjective_robust.png')
plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
print(f"\n✓ 图像已保存: {output_file}")
plt.close()
return pareto_front
pareto_front = case2_multiobjective_robust()
# =============================================================================
# 案例3：最坏情况优化与保守性分析
# =============================================================================
print("\n" + "=" * 80)
print("案例3：最坏情况优化与保守性分析")
print("=" * 80)
def case3_worst_case_optimization():
"""最坏情况优化 - 区间不确定性分析"""
print("\n【问题描述】")
print("考虑区间不确定性，使用最坏情况优化方法，")
print("分析不同保守性水平对优化结果的影响。")
# 问题设置：桁架结构优化
# 设计变量：杆件截面积 A1, A2
# 不确定性：载荷 P 在 [P_nom - ΔP, P_nom + ΔP] 范围内变化
L1, L2 = 1.0, 1.0 # 杆件长度 (m)
rho = 7850 # 密度 (kg/m³)
E = 200e9 # 弹性模量 (Pa)
sigma_allow = 250e6 # 许用应力 (Pa)
P_nom = 100e3 # 名义载荷 (N)
Delta_P = 20e3 # 载荷不确定范围 (N)
print(f"\n【参数设置】")
print(f" 杆件1长度 L1 = {L1} m")
print(f" 杆件2长度 L2 = {L2} m")
print(f" 密度 ρ = {rho} kg/m³")
print(f" 弹性模量 E = {E/1e9:.0f} GPa")
print(f" 许用应力 σ_allow = {sigma_allow/1e6:.0f} MPa")
print(f" 名义载荷 P = {P_nom/1e3:.0f} ± {Delta_P/1e3:.0f} kN")
# 两杆桁架结构
# 杆件1：与水平方向夹角 θ1 = 60°
# 杆件2：与水平方向夹角 θ2 = 45°
theta1 = np.radians(60)
theta2 = np.radians(45)
# 确定性优化
print(f"\n【确定性优化】")
def weight_det(A):
A1, A2 = A
return rho * (A1 * L1 + A2 * L2)
def stress_constraint_det(A):
A1, A2 = A
# 节点平衡方程求解内力
N2 = P_nom / (np.sin(theta1) * np.cos(theta2)/np.cos(theta1) + np.sin(theta2))
N1 = N2 * np.cos(theta2) / np.cos(theta1)
sigma1 = abs(N1) / A1
sigma2 = abs(N2) / A2
return min(sigma_allow - sigma1, sigma_allow - sigma2)
result_det = minimize(
weight_det,
[1e-4, 1e-4],
method='SLSQP',
bounds=[(1e-5, 1e-2), (1e-5, 1e-2)],
constraints={'type': 'ineq', 'fun': stress_constraint_det}
)
A1_det, A2_det = result_det.x
weight_det_val = result_det.fun
print(f" 最优设计: A1={A1_det*1e4:.2f} cm², A2={A2_det*1e4:.2f} cm²")
print(f" 重量: {weight_det_val:.2f} kg")
# 最坏情况优化（不同保守性水平）
print(f"\n【最坏情况优化】")
gammas = [0.0, 0.25, 0.5, 0.75, 1.0]
results_wc = []
for gamma in gammas:
# 有效不确定范围
Delta_P_eff = gamma * Delta_P
P_min = P_nom - Delta_P_eff
P_max = P_nom + Delta_P_eff
def worst_case_constraint(A):
A1, A2 = A
# 最坏情况：最大载荷
P_wc = P_max
N2 = P_wc / (np.sin(theta1) * np.cos(theta2)/np.cos(theta1) + np.sin(theta2))
N1 = N2 * np.cos(theta2) / np.cos(theta1)
sigma1 = abs(N1) / A1
sigma2 = abs(N2) / A2
return min(sigma_allow - sigma1, sigma_allow - sigma2)
result = minimize(
weight_det,
[1e-4, 1e-4],
method='SLSQP',
bounds=[(1e-5, 1e-2), (1e-5, 1e-2)],
constraints={'type': 'ineq', 'fun': worst_case_constraint}
)
if result.success:
A1, A2 = result.x
weight_val = result.fun
# 计算在最坏情况下的应力
N2_wc = P_max / (np.sin(theta1) * np.cos(theta2)/np.cos(theta1) + np.sin(theta2))
N1_wc = N2_wc * np.cos(theta2) / np.cos(theta1)
sigma1_wc = abs(N1_wc) / A1
sigma2_wc = abs(N2_wc) / A2
results_wc.append({
'gamma': gamma,
'A1': A1,
'A2': A2,
'weight': weight_val,
'sigma1_wc': sigma1_wc,
'sigma2_wc': sigma2_wc
})
print(f" γ={gamma:.2f}: A1={A1*1e4:.2f}cm², A2={A2*1e4:.2f}cm², "
f"W={weight_val:.2f}kg, σ_max={max(sigma1_wc, sigma2_wc)/1e6:.1f}MPa")
# ========== 可视化 ==========
fig = plt.figure(figsize=(16, 10))
# 子图1: 桁架结构示意图
ax1 = fig.add_subplot(2, 3, 1)
# 绘制桁架（确定性设计）
# 节点坐标
node0 = np.array([0, 0]) # 固定节点
node1 = np.array([L1 * np.cos(theta1), L1 * np.sin(theta1)]) # 连接节点
node2 = np.array([L2 * np.cos(theta2), -L2 * np.sin(theta2)]) # 另一个固定节点
# 绘制杆件
ax1.plot([node0[0], node1[0]], [node0[1], node1[1]], 'b-', linewidth=8, alpha=0.5, label='杆件1')
ax1.plot([node0[0], node2[0]], [node0[1], node2[1]], 'r-', linewidth=8, alpha=0.5, label='杆件2')
# 绘制节点
ax1.plot(*node0, 'ko', markersize=15)
ax1.plot(*node1, 'go', markersize=12)
ax1.plot(*node2, 'ko', markersize=15)
# 绘制支座
ax1.plot([node1[0]-0.1, node1[0]+0.1], [node1[1], node1[1]], 'k-', linewidth=3)
# 载荷箭头
ax1.arrow(node0[0], node0[1]+0.1, 0, -0.15, head_width=0.05, head_length=0.05,
fc='green', ec='green', linewidth=2)
ax1.text(node0[0]+0.1, node0[1]-0.1, 'P', fontsize=12, color='green')
ax1.set_xlim(-0.3, 1.2)
ax1.set_ylim(-0.5, 1.0)
ax1.set_aspect('equal')
ax1.set_xlabel('x (m)', fontsize=11)
ax1.set_ylabel('y (m)', fontsize=11)
ax1.set_title('两杆桁架结构', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
# 子图2: 保守性-重量关系
ax2 = fig.add_subplot(2, 3, 2)
gammas_plot = [r['gamma'] for r in results_wc]
weights_plot = [r['weight'] for r in results_wc]
ax2.plot(gammas_plot, weights_plot, 'b-o', linewidth=2, markersize=8)
ax2.fill_between(gammas_plot, weights_plot, alpha=0.3, color='blue')
ax2.set_xlabel('保守性参数 γ', fontsize=11)
ax2.set_ylabel('重量 (kg)', fontsize=11)
ax2.set_title('保守性-重量权衡', fontsize=13, fontweight='bold')
ax2.grid(True, alpha=0.3)
# 添加数值标签
for g, w in zip(gammas_plot, weights_plot):
ax2.text(g, w+0.5, f'{w:.2f}', ha='center', fontsize=9)
# 子图3: 截面积随保守性变化
ax3 = fig.add_subplot(2, 3, 3)
A1s = [r['A1']*1e4 for r in results_wc]
A2s = [r['A2']*1e4 for r in results_wc]
x = np.arange(len(gammas_plot))
width = 0.35
bars1 = ax3.bar(x - width/2, A1s, width, label='A1', color='steelblue', edgecolor='navy')
bars2 = ax3.bar(x + width/2, A2s, width, label='A2', color='coral', edgecolor='darkred')
ax3.set_xlabel('保守性参数 γ', fontsize=11)
ax3.set_ylabel('截面积 (cm²)', fontsize=11)
ax3.set_title('截面积随保守性变化', fontsize=13, fontweight='bold')
ax3.set_xticks(x)
ax3.set_xticklabels([f'{g:.2f}' for g in gammas_plot])
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, axis='y')
# 子图4: 最坏情况应力分析
ax4 = fig.add_subplot(2, 3, 4)
sigma1s = [r['sigma1_wc']/1e6 for r in results_wc]
sigma2s = [r['sigma2_wc']/1e6 for r in results_wc]
ax4.plot(gammas_plot, sigma1s, 'b-o', linewidth=2, markersize=8, label='杆件1应力')
ax4.plot(gammas_plot, sigma2s, 'r-s', linewidth=2, markersize=8, label='杆件2应力')
ax4.axhline(y=sigma_allow/1e6, color='k', linestyle='--', linewidth=2, label='许用应力')
ax4.set_xlabel('保守性参数 γ', fontsize=11)
ax4.set_ylabel('应力 (MPa)', fontsize=11)
ax4.set_title('最坏情况应力分析', fontsize=13, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
# 子图5: 载荷不确定性影响
ax5 = fig.add_subplot(2, 3, 5)
# 在不同保守性水平下，分析载荷变化对重量的影响
P_range = np.linspace(P_nom - Delta_P, P_nom + Delta_P, 100)
# 安全地选择要显示的结果
n_results = len(results_wc)
if n_results >= 3:
selected_indices = [0, n_results // 2, n_results - 1]
elif n_results >= 2:
selected_indices = [0, n_results - 1]
else:
selected_indices = [0]
colors = ['blue', 'green', 'red']
for i, idx in enumerate(selected_indices):
res = results_wc[idx]
A1, A2 = res['A1'], res['A2']
stresses = []
for P in P_range:
N2 = P / (np.sin(theta1) * np.cos(theta2)/np.cos(theta1) + np.sin(theta2))
N1 = N2 * np.cos(theta2) / np.cos(theta1)
sigma1 = abs(N1) / A1
sigma2 = abs(N2) / A2
stresses.append(max(sigma1, sigma2) / 1e6)
label = f"γ={res['gamma']:.1f}"
color = colors[i % len(colors)]
ax5.plot(P_range/1e3, stresses, color=color, linewidth=2, label=label)
ax5.axhline(y=sigma_allow/1e6, color='k', linestyle='--', linewidth=2, label='许用应力')
ax5.fill_betweenx([0, sigma_allow/1e6*1.2], P_nom/1e3 - Delta_P/1e3, P_nom/1e3 + Delta_P/1e3,
alpha=0.2, color='gray', label='不确定范围')
ax5.set_xlabel('载荷 (kN)', fontsize=11)
ax5.set_ylabel('最大应力 (MPa)', fontsize=11)
ax5.set_title('载荷-应力关系', fontsize=13, fontweight='bold')
ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3)
# 子图6: 安全裕度分析
ax6 = fig.add_subplot(2, 3, 6)
safety_margins = []
for res in results_wc:
sigma_max = max(res['sigma1_wc'], res['sigma2_wc'])
margin = (sigma_allow - sigma_max) / sigma_allow * 100
safety_margins.append(margin)
colors = ['red' if m < 0 else 'orange' if m < 5 else 'green' for m in safety_margins]
bars = ax6.bar(range(len(gammas_plot)), safety_margins, color=colors, edgecolor='black')
ax6.axhline(y=0, color='k', linestyle='-', linewidth=1)
ax6.set_xlabel('保守性参数 γ', fontsize=11)
ax6.set_ylabel('安全裕度 (%)', fontsize=11)
ax6.set_title('安全裕度分析', fontsize=13, fontweight='bold')
ax6.set_xticks(range(len(gammas_plot)))
ax6.set_xticklabels([f'{g:.2f}' for g in gammas_plot])
ax6.grid(True, alpha=0.3, axis='y')
# 添加数值标签
for i, (bar, margin) in enumerate(zip(bars, safety_margins)):
ax6.text(bar.get_x() + bar.get_width()/2., bar.get_height() + (1 if margin >= 0 else -3),
f'{margin:.1f}%', ha='center', fontsize=9)
plt.suptitle('案例3：最坏情况优化与保守性分析', fontsize=14, fontweight='bold', y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.96])
output_file = os.path.join(output_dir, 'case3_worst_case.png')
plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
print(f"\n✓ 图像已保存: {output_file}")
plt.close()
return results_wc
results_case3 = case3_worst_case_optimization()
# =============================================================================
# 案例4：基于抽样的鲁棒性优化对比
# =============================================================================
print("\n" + "=" * 80)
print("案例4：基于抽样的鲁棒性优化对比")
print("=" * 80)
def case4_sampling_based_robust():
"""基于抽样的鲁棒性优化 - 不同抽样方法对比"""
print("\n【问题描述】")
print("比较不同抽样方法（蒙特卡洛、LHS、Sobol序列）")
print("在鲁棒性优化中的效率和精度。")
# 问题设置：悬臂梁优化
L = 2.0
rho = 7850
sigma_allow = 250e6
# 不确定性参数
E_mean, sigma_E = 200e9, 20e9
P_mean, sigma_P = 10e3, 2e3
print(f"\n【参数设置】")
print(f" 梁长度 L = {L} m")
print(f" 弹性模量 E ~ N({E_mean/1e9:.0f}, {sigma_E/1e9:.1f}²) GPa")
print(f" 载荷 P ~ N({P_mean/1e3:.0f}, {sigma_P/1e3:.0f}²) kN")
print(f" 许用应力 σ_allow = {sigma_allow/1e6:.0f} MPa")
# 定义抽样方法
def lhs_sampling(n, E_mean, sigma_E, P_mean, sigma_P):
"""拉丁超立方抽样"""
# 生成LHS样本
u1 = np.random.rand(n)
u2 = np.random.rand(n)
# 分层
strata1 = np.floor(u1 * n).astype(int)
strata2 = np.floor(u2 * n).astype(int)
# 在层内随机化
v1 = (strata1 + np.random.rand(n)) / n
v2 = (strata2 + np.random.rand(n)) / n
# 转换为正态分布
E_samples = norm.ppf(v1, E_mean, sigma_E)
P_samples = norm.ppf(v2, P_mean, sigma_P)
return E_samples, P_samples
def sobol_sampling(n, E_mean, sigma_E, P_mean, sigma_P):
"""Sobol序列抽样（简化实现）"""
# 使用简单的低差异序列近似
i = np.arange(1, n+1)
# Van der Corput序列
def van_der_corput(n, base=2):
result = np.zeros(n)
for i in range(n):
n_temp = i + 1
f = 1.0
r = 0.0
while n_temp > 0:
f = f / base
r = r + f * (n_temp % base)
n_temp = n_temp // base
result[i] = r
return result
v1 = van_der_corput(n, 2)
v2 = van_der_corput(n, 3)
# 转换为正态分布
E_samples = norm.ppf(v1, E_mean, sigma_E)
P_samples = norm.ppf(v2, P_mean, sigma_P)
return E_samples, P_samples
# 计算应力的函数
def compute_stress(b, h, E, P):
I = b * h**3 / 12
M_max = P * L
sigma = M_max * (h/2) / I
return sigma
def weight(b, h):
return rho * L * b * h
# 不同抽样方法对比
print(f"\n【抽样方法对比】")
sampling_methods = {
'蒙特卡洛': lambda n: (np.random.normal(E_mean, sigma_E, n),
np.random.normal(P_mean, sigma_P, n)),
'LHS': lambda n: lhs_sampling(n, E_mean, sigma_E, P_mean, sigma_P),
'Sobol': lambda n: sobol_sampling(n, E_mean, sigma_E, P_mean, sigma_P)
}
# 不同样本量
sample_sizes = [50, 100, 200, 500, 1000]
results_comparison = {}
for method_name, sampler in sampling_methods.items():
print(f"\n【{method_name}方法】")
method_results = []
for n in sample_sizes:
# 生成样本
E_samples, P_samples = sampler(n)
# 定义鲁棒目标函数（使用样本估计）
def robust_objective(x):
b, h = x
stresses = [compute_stress(b, h, E, P) for E, P in zip(E_samples, P_samples)]
mean_stress = np.mean(stresses)
std_stress = np.std(stresses)
weight_val = weight(b, h)
# 鲁棒目标：重量 + 惩罚项
penalty = max(0, mean_stress + 3*std_stress - sigma_allow) * 1000
return weight_val + penalty
# 优化
result = minimize(
robust_objective,
[0.05, 0.1],
method='SLSQP',
bounds=[(0.01, 0.5), (0.02, 0.5)]
)
b, h = result.x
w_val = weight(b, h)
# 用大量样本验证
E_val = np.random.normal(E_mean, sigma_E, 10000)
P_val = np.random.normal(P_mean, sigma_P, 10000)
stresses_val = [compute_stress(b, h, E, P) for E, P in zip(E_val, P_val)]
mean_stress_val = np.mean(stresses_val)
std_stress_val = np.std(stresses_val)
violation_rate = np.sum(np.array(stresses_val) > sigma_allow) / 10000 * 100
method_results.append({
'n': n,
'b': b,
'h': h,
'weight': w_val,
'mean_stress': mean_stress_val,
'std_stress': std_stress_val,
'violation_rate': violation_rate
})
print(f" n={n:4d}: b={b*1000:.2f}mm, h={h*1000:.2f}mm, "
f"W={w_val:.2f}kg, σ_mean={mean_stress_val/1e6:.1f}MPa, "
f"违反率={violation_rate:.2f}%")
results_comparison[method_name] = method_results
# ========== 可视化 ==========
fig = plt.figure(figsize=(16, 12))
# 子图1: 重量收敛性对比
ax1 = fig.add_subplot(2, 3, 1)
for method_name, results in results_comparison.items():
ns = [r['n'] for r in results]
weights = [r['weight'] for r in results]
ax1.plot(ns, weights, 'o-', linewidth=2, markersize=8, label=method_name)
ax1.set_xlabel('样本量 n', fontsize=11)
ax1.set_ylabel('重量 (kg)', fontsize=11)
ax1.set_title('重量收敛性对比', fontsize=13, fontweight='bold')
ax1.set_xscale('log')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
# 子图2: 约束违反率对比
ax2 = fig.add_subplot(2, 3, 2)
for method_name, results in results_comparison.items():
ns = [r['n'] for r in results]
violations = [r['violation_rate'] for r in results]
ax2.plot(ns, violations, 's-', linewidth=2, markersize=8, label=method_name)
ax2.axhline(y=1, color='r', linestyle='--', linewidth=2, label='1%阈值')
ax2.set_xlabel('样本量 n', fontsize=11)
ax2.set_ylabel('约束违反率 (%)', fontsize=11)
ax2.set_title('约束违反率对比', fontsize=13, fontweight='bold')
ax2.set_xscale('log')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
# 子图3: 应力均值对比
ax3 = fig.add_subplot(2, 3, 3)
for method_name, results in results_comparison.items():
ns = [r['n'] for r in results]
mean_stresses = [r['mean_stress']/1e6 for r in results]
ax3.plot(ns, mean_stresses, '^-', linewidth=2, markersize=8, label=method_name)
ax3.axhline(y=sigma_allow/1e6, color='k', linestyle='--', linewidth=2, label='许用应力')
ax3.set_xlabel('样本量 n', fontsize=11)
ax3.set_ylabel('应力均值 (MPa)', fontsize=11)
ax3.set_title('应力均值收敛性', fontsize=13, fontweight='bold')
ax3.set_xscale('log')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
# 子图4: 截面尺寸对比（n=1000时）
ax4 = fig.add_subplot(2, 3, 4)
methods = list(results_comparison.keys())
bs = [results_comparison[m][-1]['b']*1000 for m in methods]
hs = [results_comparison[m][-1]['h']*1000 for m in methods]
x = np.arange(len(methods))
width = 0.35
bars1 = ax4.bar(x - width/2, bs, width, label='宽度 b', color='steelblue', edgecolor='navy')
bars2 = ax4.bar(x + width/2, hs, width, label='高度 h', color='coral', edgecolor='darkred')
ax4.set_ylabel('尺寸 (mm)', fontsize=11)
ax4.set_title('截面尺寸对比（n=1000）', fontsize=13, fontweight='bold')
ax4.set_xticks(x)
ax4.set_xticklabels(methods)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3, axis='y')
# 添加数值标签
for i, (b, h) in enumerate(zip(bs, hs)):
ax4.text(i - width/2, b + 1, f'{b:.1f}', ha='center', fontsize=9)
ax4.text(i + width/2, h + 1, f'{h:.1f}', ha='center', fontsize=9)
# 子图5: 样本分布可视化（2D散点图）
ax5 = fig.add_subplot(2, 3, 5)
np.random.seed(42)
n_plot = 100
# 蒙特卡洛
E_mc = np.random.normal(E_mean, sigma_E, n_plot)
P_mc = np.random.normal(P_mean, sigma_P, n_plot)
ax5.scatter(E_mc/1e9, P_mc/1e3, alpha=0.5, s=30, label='蒙特卡洛', color='blue')
# LHS
E_lhs, P_lhs = lhs_sampling(n_plot, E_mean, sigma_E, P_mean, sigma_P)
ax5.scatter(E_lhs/1e9, P_lhs/1e3, alpha=0.5, s=30, label='LHS', color='red')
ax5.set_xlabel('弹性模量 E (GPa)', fontsize=11)
ax5.set_ylabel('载荷 P (kN)', fontsize=11)
ax5.set_title('样本分布对比（n=100）', fontsize=13, fontweight='bold')
ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3)
# 子图6: 效率-精度权衡
ax6 = fig.add_subplot(2, 3, 6)
# 计算收敛速度（以n=1000为参考）
ref_weight = results_comparison['蒙特卡洛'][-1]['weight']
for method_name, results in results_comparison.items():
ns = [r['n'] for r in results]
errors = [abs(r['weight'] - ref_weight) / ref_weight * 100 for r in results]
ax6.plot(ns, errors, 'o-', linewidth=2, markersize=8, label=method_name)
ax6.set_xlabel('样本量 n', fontsize=11)
ax6.set_ylabel('相对误差 (%)', fontsize=11)
ax6.set_title('收敛速度对比', fontsize=13, fontweight='bold')
ax6.set_xscale('log')
ax6.set_yscale('log')
ax6.legend(fontsize=10)
ax6.grid(True, alpha=0.3)
plt.suptitle('案例4：基于抽样的鲁棒性优化对比', fontsize=14, fontweight='bold', y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.96])
output_file = os.path.join(output_dir, 'case4_sampling_comparison.png')
plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
print(f"\n✓ 图像已保存: {output_file}")
plt.close()
return results_comparison
results_case4 = case4_sampling_based_robust()
# =============================================================================
# 生成Pareto前沿演化GIF动画
# =============================================================================
print("\n" + "=" * 80)
print("生成Pareto前沿演化GIF动画")
print("=" * 80)
def create_pareto_evolution_gif():
"""创建Pareto前沿演化过程的GIF动画"""
print("\n【生成动画】")
# 重新计算Pareto前沿（用于动画）
c1, c2 = 1.0, 1.0
sigma_c1, sigma_c2 = 0.2, 0.2
n_frames = 30
weights_f1 = np.linspace(0, 1, n_frames)
# 存储每一帧的数据
frames_data = []
for frame, w1 in enumerate(weights_f1):
w2 = 1 - w1
def objective(x):
x1, x2 = x
f1 = c1 * x1**2 + c2 * x2**2
f2 = (x1**2)**2 * sigma_c1**2 + (x2**2)**2 * sigma_c2**2
return w1 * f1 + w2 * f2
def constraint(x):
return x[0] + x[1] - 1
result = minimize(
objective,
[0.5, 0.5],
method='SLSQP',
bounds=[(0, 2), (0, 2)],
constraints={'type': 'ineq', 'fun': constraint}
)
if result.success:
x1, x2 = result.x
f1 = c1 * x1**2 + c2 * x2**2
f2 = (x1**2)**2 * sigma_c1**2 + (x2**2)**2 * sigma_c2**2
frames_data.append({'frame': frame, 'w1': w1, 'f1': f1, 'f2': f2, 'x': [x1, x2]})
# 创建动画帧
temp_dir = os.path.join(output_dir, 'temp_frames')
os.makedirs(temp_dir, exist_ok=True)
frame_files = []
for i, data in enumerate(frames_data):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
# 左图：目标空间中的Pareto前沿演化
f1_vals = [d['f1'] for d in frames_data[:i+1]]
f2_vals = [d['f2'] for d in frames_data[:i+1]]
w1_vals = [d['w1'] for d in frames_data[:i+1]]
scatter = ax1.scatter(f1_vals, f2_vals, c=w1_vals, cmap='viridis',
s=100, edgecolors='black', linewidth=1.5)
# 高亮当前点
ax1.scatter(data['f1'], data['f2'], c='red', s=300, marker='*',
edgecolors='black', linewidth=2, zorder=5)
ax1.set_xlabel('性能目标 f₁', fontsize=12)
ax1.set_ylabel('方差目标 f₂', fontsize=12)
ax1.set_title(f'Pareto前沿演化 (w₁={data["w1"]:.2f})', fontsize=14, fontweight='bold')
ax1.set_xlim(0, 2.5)
ax1.set_ylim(0, 0.25)
ax1.grid(True, alpha=0.3)
plt.colorbar(scatter, ax=ax1, label='性能权重 w₁')
# 右图：设计空间中的解演化
x1_vals = [d['x'][0] for d in frames_data[:i+1]]
x2_vals = [d['x'][1] for d in frames_data[:i+1]]
scatter2 = ax2.scatter(x1_vals, x2_vals, c=w1_vals, cmap='viridis',
s=100, edgecolors='black', linewidth=1.5)
# 高亮当前点
ax2.scatter(data['x'][0], data['x'][1], c='red', s=300, marker='*',
edgecolors='black', linewidth=2, zorder=5)
# 绘制约束边界
x1_boundary = np.linspace(0, 2, 100)
x2_boundary = 1 - x1_boundary
ax2.plot(x1_boundary, x2_boundary, 'r--', linewidth=2, label='约束边界')
ax2.fill_between(x1_boundary, 0, x2_boundary, alpha=0.2, color='red')
ax2.set_xlabel('x₁', fontsize=12)
ax2.set_ylabel('x₂', fontsize=12)
ax2.set_title(f'设计空间演化 (w₁={data["w1"]:.2f})', fontsize=14, fontweight='bold')
ax2.set_xlim(0, 2)
ax2.set_ylim(0, 2)
ax2.grid(True, alpha=0.3)
ax2.legend(fontsize=10)
plt.colorbar(scatter2, ax=ax2, label='性能权重 w₁')
plt.suptitle(f'多目标鲁棒性优化演化过程 (帧 {i+1}/{len(frames_data)})',
fontsize=15, fontweight='bold')
plt.tight_layout(rect=[0, 0, 1, 0.95])
frame_file = os.path.join(temp_dir, f'frame_{i:03d}.png')
plt.savefig(frame_file, dpi=100, bbox_inches='tight', facecolor='white')
frame_files.append(frame_file)
plt.close()
if (i + 1) % 10 == 0:
print(f" 已生成 {i+1}/{len(frames_data)} 帧")
# 创建GIF
print("\n【合成GIF动画】")
images = []
for frame_file in frame_files:
img = Image.open(frame_file)
images.append(img)
gif_file = os.path.join(output_dir, 'pareto_evolution.gif')
images[0].save(
gif_file,
save_all=True,
append_images=images[1:],
duration=200,
loop=0
)
print(f"✓ GIF动画已保存: {gif_file}")
# 清理临时文件
for frame_file in frame_files:
os.remove(frame_file)
os.rmdir(temp_dir)
return gif_file
gif_file = create_pareto_evolution_gif()
# =============================================================================
# 程序结束
# =============================================================================
print("\n" + "=" * 80)
print("所有仿真案例已完成！")
print("=" * 80)
print(f"\n生成的文件:")
print(f" 1. {os.path.join(output_dir, 'case1_cantilever_robust.png')}")
print(f" 2. {os.path.join(output_dir, 'case2_multiobjective_robust.png')}")
print(f" 3. {os.path.join(output_dir, 'case3_worst_case.png')}")
print(f" 4. {os.path.join(output_dir, 'case4_sampling_comparison.png')}")
print(f" 5. {os.path.join(output_dir, 'pareto_evolution.gif')}")
print("\n" + "=" * 80)
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
2
0- 0
已为社区贡献84条内容
所有评论(0)