第七十七篇：电磁场多物理场耦合优化

摘要

电磁场多物理场耦合优化是现代工程设计中解决复杂电磁系统综合性能优化的关键技术。本主题系统介绍电磁-热-结构多物理场耦合的基本理论、耦合机制、数值求解方法以及协同优化策略。通过建立电磁场、温度场、应力场之间的耦合模型，实现电磁器件在电性能、热管理和结构可靠性等多方面的综合优化。本主题涵盖多物理场耦合理论、弱耦合与强耦合求解方法、多学科设计优化（MDO）框架、代理模型技术、灵敏度分析等内容，并通过Python实现微波器件、天线系统、功率电子等典型应用的多物理场耦合优化设计，为复杂电磁系统的综合性能优化提供理论指导和工程实践方法。

关键词

多物理场耦合，电磁-热耦合，电磁-结构耦合，多学科优化，代理模型，协同优化，灵敏度分析，Pareto优化

---

1. 多物理场耦合基础理论

1.1 耦合物理场概述

电磁场（Electromagnetic Field）：

电磁场由麦克斯韦方程组描述：

$$\nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t}$$
$$\nabla \times \mathbf{H}=\mathbf{J}+\frac{\partial \mathbf{D}}{\partial t}$$
$$\nabla \cdot \mathbf{D}=\rho$$
$$\nabla \cdot \mathbf{B}=0$$

其中 $\mathbf{E}$ 为电场强度，$\mathbf{H}$ 为磁场强度，$\mathbf{D}$ 为电位移矢量，$\mathbf{B}$ 为磁感应强度，$\mathbf{J}$ 为电流密度，$\rho$ 为电荷密度。

温度场（Thermal Field）：

温度场由热传导方程描述：

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+Q_{gen}$$

其中 $T$ 为温度，$\rho$ 为密度，$c_p$ 为比热容，$k$ 为热导率，$Q_{gen}$ 为热源项。

结构场（Structural Field）：

结构响应由弹性力学方程描述：

$$\nabla \cdot \mathbf{\sigma }+\mathbf{f}=\rho \frac{{\partial }^2\mathbf{u}}{\partial t^2}$$

其中 $\mathbf{\sigma }$ 为应力张量，$\mathbf{f}$ 为体积力，$\mathbf{u}$ 为位移矢量。

1.2 多物理场耦合机制

电磁-热耦合：

电磁损耗产生热量，形成热源：

$$Q_{EM}=\frac{1}{2}\omega {\epsilon }_0{\epsilon }^{\prime \prime }|\mathbf{E}{|}^2+\frac{1}{2}\omega {\mu }_0{\mu }^{\prime \prime }|\mathbf{H}{|}^2+\frac{1}{2}\sigma |\mathbf{E}{|}^2$$

其中 ${\epsilon }^{\prime \prime }$ 为介电损耗因子，${\mu }^{\prime \prime }$ 为磁损耗因子，$\sigma$ 为电导率。

热-结构耦合：

热膨胀产生热应力：

$${\mathbf{\epsilon }}_{thermal}=\alpha (T-T_0)\mathbf{I}$$

总应变为机械应变与热应变之和：

$${\mathbf{\epsilon }}_{total}={\mathbf{\epsilon }}_{mech}+{\mathbf{\epsilon }}_{thermal}$$

电磁-结构耦合：

电磁力（洛伦兹力）作用在结构上：

$${\mathbf{F}}_{EM}={\rho }_e\mathbf{E}+\mathbf{J}\times \mathbf{B}$$

压电效应（机电耦合）：

$$\mathbf{D}=\mathbf{d}:\mathbf{\sigma }+{\mathbf{\epsilon }}^T\cdot \mathbf{E}$$

其中 $\mathbf{d}$ 为压电常数张量。

1.3 耦合分类

单向耦合（One-way Coupling）：

物理场A影响物理场B，但B对A的影响可忽略。例如：电磁损耗产生热量，但温度变化对电磁参数影响较小。

双向耦合（Two-way Coupling）：

物理场A和B相互影响。例如：温度改变材料电磁参数，进而改变电磁场分布和损耗。

强耦合（Strong Coupling）：

物理场之间高度相互依赖，需要同时求解。例如：压电材料中的机电耦合。

弱耦合（Weak Coupling）：

物理场之间影响较弱，可以顺序求解。例如：稳态电磁-热耦合分析。

---

2. 多物理场数值求解方法

2.1 弱耦合求解策略

顺序求解法（Sequential Solution）：

1. 求解电磁场，计算损耗分布
2. 将损耗作为热源，求解温度场
3. 根据温度分布，求解结构响应

迭代求解法（Iterative Solution）：

当物理场之间存在双向耦合时，采用迭代方法：

$$\begin{array}{rl}\text{Step 1:}& {\mathbf{E}}^{(k)}\to Q^{(k)}\\ \text{Step 2:}& Q^{(k)}\to T^{(k)}\\ \text{Step 3:}& T^{(k)}\to \epsilon (T^{(k)}),\mu (T^{(k)})\\ \text{Step 4:}& \text{Check convergence: }\Vert {\mathbf{E}}^{(k)}-{\mathbf{E}}^{(k-1)}\Vert <ϵ\end{array}$$

2.2 强耦合求解策略

全耦合求解（Fully Coupled Solution）：

将所有物理场的控制方程组合成大型方程组同时求解：

$$\left[\begin{array}{ccc}{\mathbf{K}}_{EE}& {\mathbf{K}}_{ET}& {\mathbf{K}}_{ES}\\ {\mathbf{K}}_{TE}& {\mathbf{K}}_{TT}& {\mathbf{K}}_{TS}\\ {\mathbf{K}}_{SE}& {\mathbf{K}}_{ST}& {\mathbf{K}}_{SS}\end{array}\right]\left[\begin{array}{c}\mathbf{E}\\ \mathbf{T}\\ \mathbf{u}\end{array}\right]=\left[\begin{array}{c}{\mathbf{f}}_E\\ {\mathbf{f}}_T\\ {\mathbf{f}}_S\end{array}\right]$$

其中 ${\mathbf{K}}_{ij}$ 表示物理场 $i$ 和 $j$ 之间的耦合矩阵。

分区求解法（Partitioned Solution）：

保持各物理场求解器的独立性，通过界面传递数据：

$${\mathbf{K}}_{ii}{\mathbf{x}}_i={\mathbf{f}}_i-\sum _{j\ne i}{\mathbf{K}}_{ij}{\mathbf{x}}_j$$

2.3 时间耦合策略

同步时间步进（Synchronous Time Stepping）：

所有物理场使用相同的时间步长：

$$\Delta t=\min (\Delta t_{EM},\Delta t_{thermal},\Delta t_{structural})$$

异步时间步进（Asynchronous Time Stepping）：

各物理场根据自身特性选择时间步长，通过插值交换数据。

---

3. 多学科设计优化（MDO）

3.1 MDO架构

多学科可行法（MDF: Multidisciplinary Feasible）：

在每次优化迭代中，执行完整的多物理场分析，确保各学科一致性：

$$\underset{\mathbf{x}}{\min }f(\mathbf{x},\mathbf{y}(\mathbf{x}))\text{s.t.}\mathbf{g}(\mathbf{x},\mathbf{y}(\mathbf{x}))\le 0$$

其中 $\mathbf{y}$ 为状态变量（电磁场、温度、位移等）。

个体学科可行法（IDF: Individual Discipline Feasible）：

引入一致性约束，允许各学科独立分析：

$$\underset{\mathbf{x},\mathbf{y}}{\min }f(\mathbf{x},\mathbf{y})\text{s.t.}\mathbf{g}(\mathbf{x},\mathbf{y})\le 0,\mathbf{r}(\mathbf{x},\mathbf{y})=0$$

其中 $\mathbf{r}$ 为残差方程，保证各学科一致性。

协同优化（CO: Collaborative Optimization）：

将优化问题分解为系统级和学科级：

系统级：
$$\underset{{\mathbf{x}}_s,{\mathbf{y}}_s}{\min }f({\mathbf{x}}_s,{\mathbf{y}}_s)\text{s.t.}\sum _i\Vert {\mathbf{x}}_s-{\mathbf{x}}_i^{\ast }{\Vert }^2+\Vert {\mathbf{y}}_s-{\mathbf{y}}_i^{\ast }{\Vert }^2=0$$

学科级：
$$\underset{{\mathbf{x}}_i,{\mathbf{y}}_i}{\min }\Vert {\mathbf{x}}_s-{\mathbf{x}}_i{\Vert }^2+\Vert {\mathbf{y}}_s-{\mathbf{y}}_i{\Vert }^2\text{s.t.}{\mathbf{g}}_i({\mathbf{x}}_i,{\mathbf{y}}_i)\le 0$$

3.2 代理模型技术

响应面模型（RSM: Response Surface Methodology）：

使用多项式近似复杂的多物理场响应：

$$\hat{y}(\mathbf{x})={\beta }_0+\sum _i{\beta }_i{x}_i+\sum _i{\beta }_{ii}{x}_i^2+\sum _{i<j}{\beta }_{ij}{x}_i{x}_j$$

Kriging模型：

$$\hat{y}(\mathbf{x})=\mathbf{f}(\mathbf{x}{)}^T\mathbf{\beta }+Z(\mathbf{x})$$

其中 $Z(\mathbf{x})$ 为具有高斯分布的随机过程。

径向基函数（RBF）：

$$\hat{y}(\mathbf{x})=\sum _{i=1}^N{w}_i\varphi (\Vert \mathbf{x}-{\mathbf{x}}_i\Vert )$$

常用基函数包括高斯函数、多二次函数等。

3.3 灵敏度分析

直接法（Direct Method）：

$$\frac{d\mathbf{y}}{d\mathbf{x}}=-{\mathbf{K}}^{-1}\frac{\partial \mathbf{K}}{\partial \mathbf{x}}\mathbf{y}$$

伴随法（Adjoint Method）：

对于目标函数 $f(\mathbf{y}(\mathbf{x}))$：

$$\frac{df}{d\mathbf{x}}=\frac{\partial f}{\partial \mathbf{x}}+{\mathbf{\lambda }}^T\frac{\partial \mathbf{R}}{\partial \mathbf{x}}$$

其中伴随变量 $\mathbf{\lambda }$ 满足：

$${\mathbf{K}}^T\mathbf{\lambda }=-\frac{\partial f}{\partial \mathbf{y}}$$

---

4. 多物理场优化问题建模

4.1 优化问题一般形式

多目标优化：

min⁡x{f1(x),f2(x),...,fm(x)}\min_{\mathbf{x}} \{f_1(\mathbf{x}), f_2(\mathbf{x}), ..., f_m(\mathbf{x})\}xmin​{f1​(x),f2​(x),...,fm​(x)}

$$\text{s.t.}\mathbf{g}(\mathbf{x})\le 0,\mathbf{h}(\mathbf{x})=0$$

其中 $f_i$ 可包括：

• $f_{EM}$: 电磁性能（S参数、增益、效率等）
• $f_{thermal}$: 热性能（最高温度、温度均匀性等）
• $f_{struct}$: 结构性能（应力、变形、固有频率等）

约束条件：

• 几何约束：尺寸限制、制造公差
• 物理约束：材料属性范围、工作条件
• 性能约束：电磁指标、温度上限、应力限制

4.2 设计变量

几何参数：

• 天线尺寸（长度、宽度、厚度）
• 散热器翅片数量和尺寸
• 结构支撑位置

材料参数：

• 介电常数、磁导率
• 热导率、比热容
• 弹性模量、热膨胀系数

工作参数：

• 输入功率、频率
• 环境温度、冷却条件

4.3 Pareto优化

Pareto支配关系：

设计 ${\mathbf{x}}_A$ 支配设计 ${\mathbf{x}}_B$（记作 ${\mathbf{x}}_A\prec {\mathbf{x}}_B$），当且仅当：

fi(xA)≤fi(xB),∀i∈{1,...,m}f_i(\mathbf{x}_A) \leq f_i(\mathbf{x}_B), \quad \forall i \in \{1,...,m\}fi​(xA​)≤fi​(xB​),∀i∈{1,...,m}
fj(xA)<fj(xB),∃j∈{1,...,m}f_j(\mathbf{x}_A) < f_j(\mathbf{x}_B), \quad \exists j \in \{1,...,m\}fj​(xA​)<fj​(xB​),∃j∈{1,...,m}

Pareto前沿：

所有非支配解构成的集合，代表性能与鲁棒性的最优权衡。

---

5. 典型应用案例

5.1 微波功率器件优化

优化目标：

• 最大化功率传输效率
• 最小化器件最高温度
• 最小化热应力

设计变量：

• 散热基板尺寸和材料
• 匹配网络参数
• 偏置电路设计

5.2 相控阵天线热设计

优化目标：

• 保持波束指向精度
• 控制单元温度分布
• 最小化结构变形

设计变量：

• 单元间距和排列
• 冷却通道布局
• 支撑结构刚度

5.3 高频变压器优化

优化目标：

• 最大化耦合系数
• 最小化铁损和铜损
• 控制温升和绝缘应力

设计变量：

• 磁芯几何形状
• 绕组匝数和线径
• 绝缘层厚度

---

6. Python实现：多物理场耦合优化

6.1 多物理场求解器框架

class MultiphysicsSolver:
    """多物理场求解器基类"""
    
    def __init__(self, physics_models):
        self.models = physics_models
        self.coupling_matrix = None
    
    def solve_weak_coupling(self, design_vars, tol=1e-6, max_iter=100):
        """弱耦合顺序求解"""
        for iteration in range(max_iter):
            # 电磁场求解
            E_field = self.models['EM'].solve(design_vars)
            losses = self.models['EM'].compute_losses(E_field)
            
            # 温度场求解
            T_field = self.models['thermal'].solve(losses)
            
            # 结构场求解
            stress, displacement = self.models['structural'].solve(T_field)
            
            # 检查收敛
            if self.check_convergence(E_field, T_field, stress):
                break
            
            # 更新材料参数（双向耦合）
            self.update_material_properties(T_field)
        
        return E_field, T_field, stress, displacement

6.2 多学科优化框架

class MDOFramework:
    """多学科设计优化框架"""
    
    def __init__(self, disciplines, architecture='MDF'):
        self.disciplines = disciplines
        self.architecture = architecture
    
    def evaluate_objectives(self, x):
        """评估目标函数"""
        if self.architecture == 'MDF':
            # 多学科可行法
            y = self.solve_multidisciplinary(x)
            f = [discipline.objective(x, y) for discipline in self.disciplines]
        
        return f
    
    def optimize(self, x0, bounds):
        """执行优化"""
        result = optimize.minimize(
            self.evaluate_objectives, x0, 
            method='SLSQP', bounds=bounds
        )
        return result

---

7. 案例研究

7.1 案例1：微波功放模块多物理场优化

考虑电磁性能、热管理和结构应力的综合优化，通过Pareto分析找到最优设计。

7.2 案例2：天线-散热器一体化设计

优化天线辐射性能和散热效率的权衡，实现紧凑高效的热管理设计。

7.3 案例3：高频连接器热-结构优化

在高频信号完整性约束下，优化连接器的温度分布和机械可靠性。

7.4 案例4：功率分配器多目标优化

同时优化功率分配器的电磁性能、热特性和制造成本。

7.5 案例5：多物理场灵敏度分析

分析设计参数对电磁、热、结构性能的灵敏度，指导设计优化方向。

7.6 案例6：代理模型辅助优化

使用Kriging代理模型加速多物理场优化，提高计算效率。

---

8. 结果分析与讨论

8.1 优化结果验证

通过全波仿真和实验测试验证优化结果：

• 电磁性能：S参数、辐射方向图
• 热性能：红外热成像、温度传感器
• 结构性能：应变片测量、模态分析

8.2 计算效率分析

对比不同优化策略的计算成本：

• 直接优化 vs 代理模型辅助优化
• 单学科优化 vs 多学科协同优化
• 不同MDO架构的效率对比

8.3 设计空间探索

分析设计变量对多物理场性能的影响：

• 主效应分析
• 交互效应分析
• 设计空间可视化

---

9. 工程应用与扩展

9.1 5G基站天线优化

多频段、大规模MIMO天线的电磁-热-结构协同优化：

• 多频段电磁性能平衡
• 高功率密度热管理
• 风载荷结构可靠性

9.2 电动汽车无线充电系统

充电线圈的电磁-热-结构优化：

• 耦合效率最大化
• 温升控制
• 机械防护设计

9.3 卫星通信载荷

空间环境下电磁器件的多物理场优化：

• 真空热环境适应性
• 发射振动载荷
• 在轨热循环

---

10. 总结与展望

本主题系统介绍了电磁场多物理场耦合优化的理论基础和方法体系：

核心理论：

• 电磁-热-结构多物理场耦合机制
• 弱耦合与强耦合求解策略
• 多学科设计优化（MDO）框架

方法体系：

• 顺序求解与迭代求解方法
• 代理模型技术（RSM、Kriging、RBF）
• 灵敏度分析方法（直接法、伴随法）
• Pareto多目标优化

工程应用：

• 微波功率器件、相控阵天线、高频变压器优化
• 5G基站、电动汽车、卫星通信等应用场景

未来发展方向：

1. 高保真多物理场建模：结合CFD、FEA的高精度耦合分析
2. 机器学习加速：深度神经网络替代传统代理模型
3. 实时优化：数字孪生驱动的在线优化调整
4. 不确定性量化：考虑制造公差和环境变化的多物理场鲁棒优化

电磁场多物理场耦合优化为复杂电磁系统的综合性能提升提供了科学方法，在高频电子、通信系统、功率电子等领域具有重要应用价值。

附录：Python代码清单

# -*- coding: utf-8 -*-
"""
主题077：电磁场多物理场耦合优化仿真
Multiphysics Coupled Optimization for Electromagnetic Fields
"""

import matplotlib
matplotlib.use('Agg')  # 非交互式后端
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from scipy import optimize, interpolate
from scipy.spatial.distance import cdist
import warnings
warnings.filterwarnings('ignore')
import os

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

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

print("=" * 70)
print("主题077：电磁场多物理场耦合优化仿真")
print("=" * 70)
print("\n开始运行电磁场多物理场耦合优化仿真案例...")
print("=" * 70)

# =============================================================================
# 案例1：微波功放模块多物理场优化（含GIF动画）
# =============================================================================
def case1_power_amplifier_multiphysics():
    """案例1：微波功放模块多物理场优化 - Pareto优化"""
    print("\n案例1：微波功放模块多物理场优化")
    print("-" * 60)
    
    # 设计变量：散热片高度、散热片数量、基板厚度
    # x = [fin_height, num_fins, substrate_thickness]
    
    def electromagnetic_performance(x):
        """电磁性能：S11参数（越小越好）"""
        fin_height, num_fins, substrate_thickness = x
        # 散热片影响匹配
        mismatch = 0.1 * (fin_height - 5)**2 / 25 + 0.05 * (num_fins - 10)**2 / 100
        # 基板厚度影响阻抗
        impedance_deviation = 0.05 * (substrate_thickness - 1.0)**2
        return -(20 - mismatch - impedance_deviation)  # 负值因为最小化
    
    def thermal_performance(x):
        """热性能：最高温度（越小越好）"""
        fin_height, num_fins, substrate_thickness = x
        # 散热面积
        heat_dissipation = fin_height * num_fins * 0.1
        # 基板热阻
        thermal_resistance = substrate_thickness / 200
        T_max = 85 - heat_dissipation + thermal_resistance * 100
        return T_max
    
    def structural_performance(x):
        """结构性能：最大应力（越小越好）"""
        fin_height, num_fins, substrate_thickness = x
        # 高散热片增加应力集中
        stress_concentration = 0.02 * fin_height**2
        # 厚基板降低应力
        stress_reduction = substrate_thickness * 5
        sigma_max = 50 + stress_concentration - stress_reduction
        return sigma_max
    
    # 生成设计空间样本
    print("  生成设计空间样本...")
    np.random.seed(42)
    n_samples = 500
    
    fin_height_range = [2, 10]  # mm
    num_fins_range = [5, 20]
    substrate_thickness_range = [0.5, 2.0]  # mm
    
    samples = np.random.rand(n_samples, 3)
    samples[:, 0] = samples[:, 0] * (fin_height_range[1] - fin_height_range[0]) + fin_height_range[0]
    samples[:, 1] = samples[:, 1] * (num_fins_range[1] - num_fins_range[0]) + num_fins_range[0]
    samples[:, 2] = samples[:, 2] * (substrate_thickness_range[1] - substrate_thickness_range[0]) + substrate_thickness_range[0]
    
    # 评估所有样本
    f_EM = np.array([electromagnetic_performance(x) for x in samples])
    f_thermal = np.array([thermal_performance(x) for x in samples])
    f_struct = np.array([structural_performance(x) for x in samples])
    
    # 寻找Pareto前沿
    print("  计算Pareto前沿...")
    def is_pareto_dominated(point, objectives):
        """检查点是否被支配"""
        for other in objectives:
            if np.all(other <= point) and np.any(other < point):
                return True
        return False
    
    objectives = np.column_stack([f_EM, f_thermal, f_struct])
    pareto_mask = np.array([not is_pareto_dominated(obj, objectives) for obj in objectives])
    pareto_samples = samples[pareto_mask]
    pareto_objectives = objectives[pareto_mask]
    
    print(f"  找到 {len(pareto_samples)} 个Pareto最优解")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 3D Pareto前沿
    ax1 = fig.add_subplot(2, 3, 1, projection='3d')
    ax1.scatter(f_EM, f_thermal, f_struct, c='lightgray', alpha=0.3, s=20, label='All Designs')
    ax1.scatter(pareto_objectives[:, 0], pareto_objectives[:, 1], pareto_objectives[:, 2], 
               c='red', s=50, label='Pareto Front', edgecolors='black')
    ax1.set_xlabel('EM Performance (S11, dB)', fontsize=10)
    ax1.set_ylabel('Max Temperature (°C)', fontsize=10)
    ax1.set_zlabel('Max Stress (MPa)', fontsize=10)
    ax1.set_title('3D Pareto Front', fontsize=12, fontweight='bold')
    ax1.legend()
    
    # 电磁-热权衡
    ax2 = fig.add_subplot(2, 3, 2)
    ax2.scatter(f_EM, f_thermal, c='lightgray', alpha=0.3, s=20)
    ax2.scatter(pareto_objectives[:, 0], pareto_objectives[:, 1], c='red', s=50, edgecolors='black')
    ax2.set_xlabel('EM Performance (S11, dB)', fontsize=10)
    ax2.set_ylabel('Max Temperature (°C)', fontsize=10)
    ax2.set_title('EM-Thermal Trade-off', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 热-结构权衡
    ax3 = fig.add_subplot(2, 3, 3)
    ax3.scatter(f_thermal, f_struct, c='lightgray', alpha=0.3, s=20)
    ax3.scatter(pareto_objectives[:, 1], pareto_objectives[:, 2], c='red', s=50, edgecolors='black')
    ax3.set_xlabel('Max Temperature (°C)', fontsize=10)
    ax3.set_ylabel('Max Stress (MPa)', fontsize=10)
    ax3.set_title('Thermal-Structural Trade-off', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 设计变量分布
    ax4 = fig.add_subplot(2, 3, 4)
    ax4.scatter(pareto_samples[:, 0], pareto_samples[:, 1], c=pareto_objectives[:, 0], 
               cmap='viridis', s=50, edgecolors='black')
    ax4.set_xlabel('Fin Height (mm)', fontsize=10)
    ax4.set_ylabel('Number of Fins', fontsize=10)
    ax4.set_title('Pareto Designs: Fin Height vs Count', fontsize=12, fontweight='bold')
    cbar = plt.colorbar(ax4.collections[0], ax=ax4)
    cbar.set_label('EM Performance', fontsize=9)
    
    # 基板厚度分布
    ax5 = fig.add_subplot(2, 3, 5)
    ax5.hist(pareto_samples[:, 2], bins=15, color='steelblue', edgecolor='black', alpha=0.7)
    ax5.set_xlabel('Substrate Thickness (mm)', fontsize=10)
    ax5.set_ylabel('Frequency', fontsize=10)
    ax5.set_title('Distribution of Substrate Thickness', fontsize=12, fontweight='bold')
    ax5.grid(True, alpha=0.3)
    
    # 目标函数分布
    ax6 = fig.add_subplot(2, 3, 6)
    ax6.boxplot([f_EM, f_thermal, f_struct], labels=['EM', 'Thermal', 'Structural'])
    ax6.set_ylabel('Objective Value', fontsize=10)
    ax6.set_title('Objective Function Distributions', fontsize=12, fontweight='bold')
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case1_power_amplifier_pareto.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    # 创建GIF动画：Pareto前沿演化
    print("  正在生成Pareto前沿演化动画...")
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
    
    n_frames = 50
    
    def init():
        ax1.clear()
        ax2.clear()
        ax1.set_xlabel('EM Performance (S11, dB)', fontsize=11)
        ax1.set_ylabel('Max Temperature (°C)', fontsize=11)
        ax1.set_title('Pareto Front Evolution', fontsize=13, fontweight='bold')
        ax1.grid(True, alpha=0.3)
        ax1.set_xlim(f_EM.min()-1, f_EM.max()+1)
        ax1.set_ylim(f_thermal.min()-5, f_thermal.max()+5)
        
        ax2.set_xlabel('Iteration', fontsize=11)
        ax2.set_ylabel('Pareto Solutions Count', fontsize=11)
        ax2.set_title('Pareto Set Growth', fontsize=13, fontweight='bold')
        ax2.grid(True, alpha=0.3)
        ax2.set_xlim(0, n_frames)
        ax2.set_ylim(0, len(pareto_samples)+10)
        
        return []
    
    pareto_indices = np.where(pareto_mask)[0]
    
    def update(frame):
        ax1.clear()
        ax2.clear()
        
        # 当前帧显示的Pareto点数量
        n_show = min((frame + 1) * len(pareto_samples) // n_frames + 1, len(pareto_samples))
        current_indices = pareto_indices[:n_show]
        
        # 显示所有样本
        ax1.scatter(f_EM, f_thermal, c='lightgray', alpha=0.2, s=15)
        # 显示当前Pareto点
        ax1.scatter(f_EM[current_indices], f_thermal[current_indices], 
                   c='red', s=60, edgecolors='black', alpha=0.8, label='Pareto Front')
        
        ax1.set_xlabel('EM Performance (S11, dB)', fontsize=11)
        ax1.set_ylabel('Max Temperature (°C)', fontsize=11)
        ax1.set_title(f'Pareto Front Evolution (Frame {frame+1}/{n_frames})', fontsize=13, fontweight='bold')
        ax1.grid(True, alpha=0.3)
        ax1.legend()
        ax1.set_xlim(f_EM.min()-1, f_EM.max()+1)
        ax1.set_ylim(f_thermal.min()-5, f_thermal.max()+5)
        
        # Pareto解数量增长
        ax2.plot(range(frame+1), [(i+1) * len(pareto_samples) // n_frames + 1 for i in range(frame+1)], 
                'b-', linewidth=2)
        ax2.set_xlabel('Iteration', fontsize=11)
        ax2.set_ylabel('Pareto Solutions Count', fontsize=11)
        ax2.set_title('Pareto Set Growth', fontsize=13, fontweight='bold')
        ax2.grid(True, alpha=0.3)
        ax2.set_xlim(0, n_frames)
        ax2.set_ylim(0, len(pareto_samples)+10)
        
        return []
    
    anim = FuncAnimation(fig, update, init_func=init, frames=n_frames, 
                        interval=100, blit=False, repeat=True)
    anim.save(os.path.join(output_dir, 'case1_pareto_evolution.gif'), 
             writer=PillowWriter(fps=10), dpi=100)
    plt.close()
    
    print("  ✓ 案例1完成：微波功放模块多物理场优化")
    return True

# =============================================================================
# 案例2：天线-散热器一体化设计（含GIF动画）
# =============================================================================
def case2_antenna_heatsink_design():
    """案例2：天线-散热器一体化设计优化"""
    print("\n案例2：天线-散热器一体化设计")
    print("-" * 60)
    
    # 设计参数
    # x = [antenna_length, heatsink_size, fin_density, base_thickness]
    
    def antenna_gain(x):
        """天线增益（dBi）"""
        antenna_length, heatsink_size, fin_density, base_thickness = x
        # 天线长度影响增益
        gain = 2.15 + 10 * np.log10(antenna_length / 50)  # 相对于半波偶极子
        # 散热器对辐射的遮挡
        blockage = 0.5 * (heatsink_size / 100) ** 2
        return gain - blockage
    
    def thermal_efficiency(x):
        """散热效率（%）"""
        antenna_length, heatsink_size, fin_density, base_thickness = x
        # 散热面积
        surface_area = heatsink_size**2 * (1 + fin_density * 0.1)
        # 基板热阻
        thermal_resistance = base_thickness / 5
        efficiency = min(95, surface_area / 50 - thermal_resistance)
        return max(10, efficiency)
    
    def compactness(x):
        """紧凑度（体积倒数，越大越好）"""
        antenna_length, heatsink_size, fin_density, base_thickness = x
        volume = antenna_length * heatsink_size * base_thickness / 1000  # cm³
        return 100 / volume
    
    # 设计空间探索
    print("  探索设计空间...")
    n_points = 30
    antenna_lengths = np.linspace(30, 100, n_points)  # mm
    heatsink_sizes = np.linspace(20, 80, n_points)    # mm
    
    AL, HS = np.meshgrid(antenna_lengths, heatsink_sizes)
    
    # 固定其他参数
    fin_density = 0.5
    base_thickness = 2.0
    
    Gain = np.zeros_like(AL)
    Thermal = np.zeros_like(AL)
    Compact = np.zeros_like(AL)
    
    for i in range(n_points):
        for j in range(n_points):
            x = [AL[i,j], HS[i,j], fin_density, base_thickness]
            Gain[i,j] = antenna_gain(x)
            Thermal[i,j] = thermal_efficiency(x)
            Compact[i,j] = compactness(x)
    
    # 创建可视化
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 增益分布
    im1 = axes[0,0].contourf(AL, HS, Gain, levels=20, cmap='viridis')
    axes[0,0].set_xlabel('Antenna Length (mm)', fontsize=10)
    axes[0,0].set_ylabel('Heatsink Size (mm)', fontsize=10)
    axes[0,0].set_title('Antenna Gain (dBi)', fontsize=12, fontweight='bold')
    plt.colorbar(im1, ax=axes[0,0])
    
    # 散热效率分布
    im2 = axes[0,1].contourf(AL, HS, Thermal, levels=20, cmap='hot')
    axes[0,1].set_xlabel('Antenna Length (mm)', fontsize=10)
    axes[0,1].set_ylabel('Heatsink Size (mm)', fontsize=10)
    axes[0,1].set_title('Thermal Efficiency (%)', fontsize=12, fontweight='bold')
    plt.colorbar(im2, ax=axes[0,1])
    
    # 紧凑度分布
    im3 = axes[0,2].contourf(AL, HS, Compact, levels=20, cmap='coolwarm')
    axes[0,2].set_xlabel('Antenna Length (mm)', fontsize=10)
    axes[0,2].set_ylabel('Heatsink Size (mm)', fontsize=10)
    axes[0,2].set_title('Compactness Index', fontsize=12, fontweight='bold')
    plt.colorbar(im3, ax=axes[0,2])
    
    # 多目标权衡曲面
    ax4 = fig.add_subplot(2, 3, 4, projection='3d')
    surf = ax4.plot_surface(AL, HS, Gain, alpha=0.7, cmap='viridis')
    ax4.set_xlabel('Antenna Length (mm)', fontsize=9)
    ax4.set_ylabel('Heatsink Size (mm)', fontsize=9)
    ax4.set_zlabel('Gain (dBi)', fontsize=9)
    ax4.set_title('Gain Surface', fontsize=12, fontweight='bold')
    
    # Pareto前沿投影
    ax5 = axes[1,1]
    # 归一化目标函数
    Gain_norm = (Gain - Gain.min()) / (Gain.max() - Gain.min())
    Thermal_norm = (Thermal - Thermal.min()) / (Thermal.max() - Thermal.min())
    Compact_norm = (Compact - Compact.min()) / (Compact.max() - Compact.min())
    
    # 综合性能指标
    Overall = 0.4 * Gain_norm + 0.4 * Thermal_norm + 0.2 * Compact_norm
    im5 = ax5.contourf(AL, HS, Overall, levels=20, cmap='RdYlGn')
    ax5.set_xlabel('Antenna Length (mm)', fontsize=10)
    ax5.set_ylabel('Heatsink Size (mm)', fontsize=10)
    ax5.set_title('Overall Performance Index', fontsize=12, fontweight='bold')
    plt.colorbar(im5, ax=ax5)
    
    # 最优设计点
    max_idx = np.unravel_index(np.argmax(Overall), Overall.shape)
    ax5.plot(AL[max_idx], HS[max_idx], 'r*', markersize=20, markeredgecolor='black', markeredgewidth=2)
    
    # 设计参数敏感性
    ax6 = axes[1,2]
    param_names = ['Length', 'Size', 'Fin Density', 'Thickness']
    param_ranges = [[30, 100], [20, 80], [0.1, 1.0], [1.0, 5.0]]
    sensitivities = []
    
    base_design = [65, 50, 0.5, 2.0]
    base_performance = 0.4 * antenna_gain(base_design) + 0.4 * thermal_efficiency(base_design) + 0.2 * compactness(base_design)
    
    for i, (name, range_val) in enumerate(zip(param_names, param_ranges)):
        perturbed = base_design.copy()
        perturbed[i] = range_val[1]
        pert_perf = 0.4 * antenna_gain(perturbed) + 0.4 * thermal_efficiency(perturbed) + 0.2 * compactness(perturbed)
        sensitivities.append(abs(pert_perf - base_performance))
    
    ax6.barh(param_names, sensitivities, color='steelblue', edgecolor='black')
    ax6.set_xlabel('Sensitivity', fontsize=10)
    ax6.set_title('Design Parameter Sensitivity', fontsize=12, fontweight='bold')
    ax6.grid(True, alpha=0.3, axis='x')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case2_antenna_heatsink.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    # 创建GIF动画：设计空间探索
    print("  正在生成设计空间探索动画...")
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
    
    n_frames = 40
    
    def init():
        ax1.clear()
        ax2.clear()
        return []
    
    def update(frame):
        ax1.clear()
        ax2.clear()
        
        # 当前探索进度
        progress = (frame + 1) / n_frames
        n_explored = int(n_points * n_points * progress)
        
        # 绘制已探索区域
        explored = np.zeros_like(Overall)
        for idx in range(min(n_explored, n_points * n_points)):
            i, j = divmod(idx, n_points)
            explored[i, j] = Overall[i, j]
        
        im1 = ax1.contourf(AL, HS, explored, levels=20, cmap='RdYlGn', vmin=Overall.min(), vmax=Overall.max())
        ax1.set_xlabel('Antenna Length (mm)', fontsize=11)
        ax1.set_ylabel('Heatsink Size (mm)', fontsize=11)
        ax1.set_title(f'Design Space Exploration ({int(progress*100)}%)', fontsize=13, fontweight='bold')
        
        # 当前最优值
        if n_explored > 0:
            current_max = np.max(explored[explored > 0])
            ax2.bar(['Current Best'], [current_max], color='green', alpha=0.7, edgecolor='black')
            ax2.axhline(y=Overall.max(), color='red', linestyle='--', linewidth=2, label='Global Optimum')
        
        ax2.set_ylabel('Performance Index', fontsize=11)
        ax2.set_title('Optimization Progress', fontsize=13, fontweight='bold')
        ax2.set_ylim(0, Overall.max() * 1.1)
        ax2.legend()
        ax2.grid(True, alpha=0.3, axis='y')
        
        return []
    
    anim = FuncAnimation(fig, update, init_func=init, frames=n_frames, 
                        interval=100, blit=False, repeat=True)
    anim.save(os.path.join(output_dir, 'case2_design_exploration.gif'), 
             writer=PillowWriter(fps=10), dpi=100)
    plt.close()
    
    print(f"  最优设计：天线长度={AL[max_idx]:.1f}mm，散热器尺寸={HS[max_idx]:.1f}mm")
    print("  ✓ 案例2完成：天线-散热器一体化设计")
    return True

# =============================================================================
# 案例3：高频连接器热-结构优化
# =============================================================================
def case3_connector_thermal_structural():
    """案例3：高频连接器热-结构优化"""
    print("\n案例3：高频连接器热-结构优化")
    print("-" * 60)
    
    # 设计变量：接触面积、绝缘层厚度、外壳厚度、散热槽数量
    
    def signal_integrity(x):
        """信号完整性：阻抗匹配度（%）"""
        contact_area, insulation_thickness, shell_thickness, num_slots = x
        # 接触面积影响阻抗
        impedance_match = 100 - 5 * abs(contact_area - 2.0)
        # 绝缘层厚度影响
        insulation_effect = 100 - 10 * abs(insulation_thickness - 1.5)
        return min(100, max(50, (impedance_match + insulation_effect) / 2))
    
    def thermal_resistance(x):
        """热阻（K/W）- 越小越好"""
        contact_area, insulation_thickness, shell_thickness, num_slots = x
        # 接触面积增大降低热阻
        contact_resistance = 5.0 / contact_area
        # 外壳厚度增加热阻
        shell_resistance = 0.5 * shell_thickness
        # 散热槽降低热阻
        slot_effect = -0.3 * num_slots
        return max(0.5, contact_resistance + shell_resistance + slot_effect)
    
    def mechanical_strength(x):
        """机械强度（N）"""
        contact_area, insulation_thickness, shell_thickness, num_slots = x
        # 外壳厚度增加强度
        shell_strength = 200 * shell_thickness
        # 散热槽削弱强度
        slot_penalty = 5 * num_slots
        return max(50, shell_strength - slot_penalty)
    
    # 多目标优化
    print("  执行多目标优化...")
    
    def objective(x):
        """综合目标：最大化信号完整性，最小化热阻，最大化机械强度"""
        sig_int = signal_integrity(x)
        therm_res = thermal_resistance(x)
        mech_str = mechanical_strength(x)
        
        # 归一化并组合（最小化问题）
        f1 = -sig_int / 100  # 信号完整性（负值因为最小化）
        f2 = therm_res / 10  # 热阻
        f3 = -mech_str / 500  # 机械强度（负值）
        
        return [f1, f2, f3]
    
    # 使用NSGA-II风格的多目标优化（简化版）
    np.random.seed(42)
    n_population = 100
    n_generations = 50
    
    # 初始化种群
    bounds = [(0.5, 4.0), (0.5, 3.0), (0.5, 3.0), (0, 10)]
    population = np.random.rand(n_population, 4)
    for i, (lb, ub) in enumerate(bounds):
        population[:, i] = population[:, i] * (ub - lb) + lb
    
    # 评估初始种群
    objectives = np.array([objective(ind) for ind in population])
    
    # 简化进化过程
    for gen in range(n_generations):
        # 选择、交叉、变异（简化）
        new_population = population.copy()
        for i in range(n_population):
            if np.random.rand() < 0.1:  # 变异概率
                idx = np.random.randint(4)
                lb, ub = bounds[idx]
                new_population[i, idx] = np.random.rand() * (ub - lb) + lb
        
        # 评估新种群
        new_objectives = np.array([objective(ind) for ind in new_population])
        
        # 合并并选择
        combined_pop = np.vstack([population, new_population])
        combined_obj = np.vstack([objectives, new_objectives])
        
        # 非支配排序选择
        def non_dominated_sort(obj_values):
            """非支配排序"""
            n = len(obj_values)
            domination_count = np.zeros(n)
            dominated_solutions = [[] for _ in range(n)]
            ranks = np.zeros(n)
            
            for i in range(n):
                for j in range(i+1, n):
                    if np.all(obj_values[i] <= obj_values[j]) and np.any(obj_values[i] < obj_values[j]):
                        dominated_solutions[i].append(j)
                        domination_count[j] += 1
                    elif np.all(obj_values[j] <= obj_values[i]) and np.any(obj_values[j] < obj_values[i]):
                        dominated_solutions[j].append(i)
                        domination_count[i] += 1
            
            current_rank = 0
            current_front = np.where(domination_count == 0)[0]
            
            while len(current_front) > 0:
                ranks[current_front] = current_rank
                next_front = []
                for i in current_front:
                    for j in dominated_solutions[i]:
                        domination_count[j] -= 1
                        if domination_count[j] == 0:
                            next_front.append(j)
                current_front = np.array(next_front)
                current_rank += 1
            
            return ranks
        
        ranks = non_dominated_sort(combined_obj)
        
        # 选择前n_population个
        selected_indices = np.argsort(ranks)[:n_population]
        population = combined_pop[selected_indices]
        objectives = combined_obj[selected_indices]
    
    # 最终Pareto前沿
    final_ranks = non_dominated_sort(objectives)
    pareto_indices = np.where(final_ranks == 0)[0]
    pareto_population = population[pareto_indices]
    pareto_objectives = objectives[pareto_indices]
    
    print(f"  最终Pareto解数量：{len(pareto_population)}")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 10))
    
    # 3D Pareto前沿
    ax1 = fig.add_subplot(2, 3, 1, projection='3d')
    ax1.scatter(objectives[:, 0], objectives[:, 1], objectives[:, 2], 
               c='lightgray', alpha=0.3, s=20, label='All Solutions')
    ax1.scatter(pareto_objectives[:, 0], pareto_objectives[:, 1], pareto_objectives[:, 2], 
               c='red', s=60, label='Pareto Front', edgecolors='black')
    ax1.set_xlabel('Signal Integrity (norm)', fontsize=9)
    ax1.set_ylabel('Thermal Resistance (norm)', fontsize=9)
    ax1.set_zlabel('Mechanical Strength (norm)', fontsize=9)
    ax1.set_title('3D Pareto Front', fontsize=12, fontweight='bold')
    ax1.legend()
    
    # 信号完整性 vs 热阻
    ax2 = fig.add_subplot(2, 3, 2)
    ax2.scatter(objectives[:, 0], objectives[:, 1], c='lightgray', alpha=0.3, s=20)
    ax2.scatter(pareto_objectives[:, 0], pareto_objectives[:, 1], c='red', s=60, edgecolors='black')
    ax2.set_xlabel('Signal Integrity (normalized)', fontsize=10)
    ax2.set_ylabel('Thermal Resistance (normalized)', fontsize=10)
    ax2.set_title('Signal vs Thermal Trade-off', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 热阻 vs 机械强度
    ax3 = fig.add_subplot(2, 3, 3)
    ax3.scatter(objectives[:, 1], objectives[:, 2], c='lightgray', alpha=0.3, s=20)
    ax3.scatter(pareto_objectives[:, 1], pareto_objectives[:, 2], c='red', s=60, edgecolors='black')
    ax3.set_xlabel('Thermal Resistance (normalized)', fontsize=10)
    ax3.set_ylabel('Mechanical Strength (normalized)', fontsize=10)
    ax3.set_title('Thermal vs Mechanical Trade-off', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 设计变量分布
    ax4 = fig.add_subplot(2, 3, 4)
    ax4.scatter(pareto_population[:, 0], pareto_population[:, 1], 
               c=pareto_objectives[:, 0], cmap='viridis', s=50, edgecolors='black')
    ax4.set_xlabel('Contact Area (mm²)', fontsize=10)
    ax4.set_ylabel('Insulation Thickness (mm)', fontsize=10)
    ax4.set_title('Pareto Designs: Contact vs Insulation', fontsize=12, fontweight='bold')
    plt.colorbar(ax4.collections[0], ax=ax4, label='Signal Integrity')
    
    # 外壳设计空间
    ax5 = fig.add_subplot(2, 3, 5)
    ax5.scatter(pareto_population[:, 2], pareto_population[:, 3], 
               c=pareto_objectives[:, 1], cmap='hot', s=50, edgecolors='black')
    ax5.set_xlabel('Shell Thickness (mm)', fontsize=10)
    ax5.set_ylabel('Number of Slots', fontsize=10)
    ax5.set_title('Pareto Designs: Shell vs Slots', fontsize=12, fontweight='bold')
    plt.colorbar(ax5.collections[0], ax=ax5, label='Thermal Resistance')
    
    # 收敛曲线
    ax6 = fig.add_subplot(2, 3, 6)
    # 模拟收敛过程
    generations = np.arange(n_generations)
    pareto_size_history = 5 + 15 * (1 - np.exp(-generations / 15))
    ax6.plot(generations, pareto_size_history, 'b-', linewidth=2, label='Pareto Set Size')
    ax6.set_xlabel('Generation', fontsize=10)
    ax6.set_ylabel('Pareto Solutions Count', fontsize=10)
    ax6.set_title('Optimization Convergence', fontsize=12, fontweight='bold')
    ax6.grid(True, alpha=0.3)
    ax6.legend()
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case3_connector_optimization.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print("  ✓ 案例3完成：高频连接器热-结构优化")
    return True

# =============================================================================
# 案例4：功率分配器多目标优化（含GIF动画）
# =============================================================================
def case4_power_divider_optimization():
    """案例4：功率分配器多目标优化"""
    print("\n案例4：功率分配器多目标优化")
    print("-" * 60)
    
    # 设计变量：传输线宽度、隔离电阻值、基板厚度、工作频率
    
    def insertion_loss(x):
        """插入损耗（dB）- 越小越好"""
        line_width, resistor_value, substrate_thickness, frequency = x
        # 传输线损耗
        conductor_loss = 0.01 * frequency / line_width
        # 基板损耗
        dielectric_loss = 0.001 * frequency * substrate_thickness
        return conductor_loss + dielectric_loss
    
    def isolation(x):
        """隔离度（dB）- 越大越好"""
        line_width, resistor_value, substrate_thickness, frequency = x
        # 隔离电阻影响
        isolation_db = 20 * np.log10(resistor_value / 50 + 1e-10)
        # 频率影响
        freq_effect = -0.1 * abs(frequency - 2.4)
        return max(10, isolation_db + freq_effect)
    
    def amplitude_balance(x):
        """幅度平衡（dB）- 越小越好"""
        line_width, resistor_value, substrate_thickness, frequency = x
        # 制造公差影响
        manufacturing_tolerance = 0.05 / line_width
        return manufacturing_tolerance
    
    def cost(x):
        """制造成本（相对值）- 越小越好"""
        line_width, resistor_value, substrate_thickness, frequency = x
        # 材料成本
        material_cost = substrate_thickness * 10
        # 精度成本
        precision_cost = 5 / line_width
        return material_cost + precision_cost
    
    # 多目标优化
    print("  执行四目标优化...")
    
    np.random.seed(42)
    n_samples = 800
    
    # 设计空间
    bounds = [(0.5, 3.0), (80, 150), (0.5, 2.0), (1.0, 5.0)]
    samples = np.random.rand(n_samples, 4)
    for i, (lb, ub) in enumerate(bounds):
        samples[:, i] = samples[:, i] * (ub - lb) + lb
    
    # 评估所有样本
    f1 = np.array([insertion_loss(x) for x in samples])
    f2 = np.array([isolation(x) for x in samples])
    f3 = np.array([amplitude_balance(x) for x in samples])
    f4 = np.array([cost(x) for x in samples])
    
    # 归一化（处理可能的除零情况）
    def safe_normalize(x, invert=False):
        x_range = x.max() - x.min()
        if x_range < 1e-10:
            return np.zeros_like(x)
        normalized = (x - x.min()) / x_range
        if invert:
            normalized = 1 - normalized
        return normalized
    
    f1_norm = safe_normalize(f1)
    f2_norm = safe_normalize(f2, invert=True)  # 反转（越大越好）
    f3_norm = safe_normalize(f3)
    f4_norm = safe_normalize(f4)
    
    objectives_norm = np.column_stack([f1_norm, f2_norm, f3_norm, f4_norm])
    
    # 寻找Pareto前沿
    def is_pareto_dominated(point, objectives):
        for other in objectives:
            if np.all(other <= point) and np.any(other < point):
                return True
        return False
    
    pareto_mask = np.array([not is_pareto_dominated(obj, objectives_norm) for obj in objectives_norm])
    pareto_samples = samples[pareto_mask]
    pareto_objectives = objectives_norm[pareto_mask]
    
    print(f"  找到 {len(pareto_samples)} 个Pareto最优解")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 平行坐标图
    ax1 = fig.add_subplot(2, 3, 1)
    obj_names = ['Insertion Loss', 'Isolation', 'Amplitude Balance', 'Cost']
    for i in range(min(50, len(pareto_objectives))):
        ax1.plot(range(4), pareto_objectives[i], alpha=0.5, color='steelblue')
    ax1.set_xticks(range(4))
    ax1.set_xticklabels(obj_names, rotation=15, ha='right')
    ax1.set_ylabel('Normalized Objective', fontsize=10)
    ax1.set_title('Parallel Coordinates Plot', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 雷达图示例
    ax2 = fig.add_subplot(2, 3, 2, projection='polar')
    angles = np.linspace(0, 2*np.pi, 4, endpoint=False).tolist()
    angles += angles[:1]
    
    # 绘制几个典型解
    for i in [0, len(pareto_objectives)//2, -1]:
        values = pareto_objectives[i].tolist()
        values += values[:1]
        ax2.plot(angles, values, 'o-', linewidth=2, label=f'Design {i}')
        ax2.fill(angles, values, alpha=0.1)
    ax2.set_xticks(angles[:-1])
    ax2.set_xticklabels(obj_names)
    ax2.set_title('Radar Chart of Designs', fontsize=12, fontweight='bold')
    ax2.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0))
    
    # 目标函数相关性
    ax3 = fig.add_subplot(2, 3, 3)
    corr_matrix = np.corrcoef(objectives_norm.T)
    im = ax3.imshow(corr_matrix, cmap='coolwarm', vmin=-1, vmax=1)
    ax3.set_xticks(range(4))
    ax3.set_yticks(range(4))
    ax3.set_xticklabels(obj_names, rotation=45, ha='right')
    ax3.set_yticklabels(obj_names)
    ax3.set_title('Objective Correlation Matrix', fontsize=12, fontweight='bold')
    plt.colorbar(im, ax=ax3)
    
    # 插入损耗 vs 隔离度
    ax4 = fig.add_subplot(2, 3, 4)
    scatter = ax4.scatter(f1, f2, c=f4, cmap='viridis', alpha=0.6, s=30)
    ax4.scatter(f1[pareto_mask], f2[pareto_mask], c='red', s=60, edgecolors='black', label='Pareto')
    ax4.set_xlabel('Insertion Loss (dB)', fontsize=10)
    ax4.set_ylabel('Isolation (dB)', fontsize=10)
    ax4.set_title('Insertion Loss vs Isolation', fontsize=12, fontweight='bold')
    ax4.legend()
    plt.colorbar(scatter, ax=ax4, label='Cost')
    
    # 设计变量分布
    ax5 = fig.add_subplot(2, 3, 5)
    ax5.scatter(samples[:, 0], samples[:, 1], c='lightgray', alpha=0.3, s=20)
    ax5.scatter(pareto_samples[:, 0], pareto_samples[:, 1], c='red', s=50, edgecolors='black')
    ax5.set_xlabel('Line Width (mm)', fontsize=10)
    ax5.set_ylabel('Resistor Value (Ω)', fontsize=10)
    ax5.set_title('Design Space: Line Width vs Resistor', fontsize=12, fontweight='bold')
    ax5.grid(True, alpha=0.3)
    
    # 综合性能分布
    ax6 = fig.add_subplot(2, 3, 6)
    overall_performance = 1 - np.mean(objectives_norm, axis=1)
    # 过滤NaN值
    overall_performance = overall_performance[np.isfinite(overall_performance)]
    if len(overall_performance) > 0:
        ax6.hist(overall_performance, bins=30, color='steelblue', edgecolor='black', alpha=0.7)
        pareto_overall = 1 - np.mean(objectives_norm[pareto_mask], axis=1)
        pareto_overall = pareto_overall[np.isfinite(pareto_overall)]
        if len(pareto_overall) > 0:
            ax6.axvline(np.mean(pareto_overall), color='red', linestyle='--', 
                       linewidth=2, label='Pareto Mean')
    ax6.set_xlabel('Overall Performance', fontsize=10)
    ax6.set_ylabel('Frequency', fontsize=10)
    ax6.set_title('Overall Performance Distribution', fontsize=12, fontweight='bold')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case4_power_divider.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    # 创建GIF动画：Pareto前沿演化
    print("  正在生成Pareto演化动画...")
    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 12))
    
    n_frames = 50
    
    def init():
        for ax in [ax1, ax2, ax3, ax4]:
            ax.clear()
        return []
    
    def update(frame):
        for ax in [ax1, ax2, ax3, ax4]:
            ax.clear()
        
        progress = (frame + 1) / n_frames
        n_pareto_show = max(1, int(len(pareto_samples) * progress))
        
        # 子图1: 插入损耗 vs 隔离度
        ax1.scatter(f1, f2, c='lightgray', alpha=0.2, s=15)
        ax1.scatter(f1[pareto_mask][:n_pareto_show], f2[pareto_mask][:n_pareto_show], 
                   c='red', s=50, edgecolors='black')
        ax1.set_xlabel('Insertion Loss (dB)', fontsize=10)
        ax1.set_ylabel('Isolation (dB)', fontsize=10)
        ax1.set_title(f'Insertion Loss vs Isolation ({n_pareto_show}/{len(pareto_samples)})', fontsize=11, fontweight='bold')
        ax1.grid(True, alpha=0.3)
        
        # 子图2: 幅度平衡 vs 成本
        ax2.scatter(f3, f4, c='lightgray', alpha=0.2, s=15)
        ax2.scatter(f3[pareto_mask][:n_pareto_show], f4[pareto_mask][:n_pareto_show], 
                   c='blue', s=50, edgecolors='black')
        ax2.set_xlabel('Amplitude Balance (dB)', fontsize=10)
        ax2.set_ylabel('Cost (relative)', fontsize=10)
        ax2.set_title(f'Balance vs Cost ({n_pareto_show}/{len(pareto_samples)})', fontsize=11, fontweight='bold')
        ax2.grid(True, alpha=0.3)
        
        # 子图3: 设计变量空间
        ax3.scatter(samples[:, 0], samples[:, 2], c='lightgray', alpha=0.2, s=15)
        ax3.scatter(pareto_samples[:n_pareto_show, 0], pareto_samples[:n_pareto_show, 2], 
                   c='green', s=50, edgecolors='black')
        ax3.set_xlabel('Line Width (mm)', fontsize=10)
        ax3.set_ylabel('Substrate Thickness (mm)', fontsize=10)
        ax3.set_title(f'Design Space ({n_pareto_show}/{len(pareto_samples)})', fontsize=11, fontweight='bold')
        ax3.grid(True, alpha=0.3)
        
        # 子图4: 进度条
        ax4.barh(['Progress'], [progress*100], color='steelblue', edgecolor='black')
        ax4.set_xlim(0, 100)
        ax4.set_xlabel('Completion (%)', fontsize=10)
        ax4.set_title('Optimization Progress', fontsize=11, fontweight='bold')
        ax4.text(progress*100/2, 0, f'{progress*100:.0f}%', ha='center', va='center', fontsize=12, fontweight='bold')
        
        plt.tight_layout()
        return []
    
    anim = FuncAnimation(fig, update, init_func=init, frames=n_frames, 
                        interval=100, blit=False, repeat=True)
    anim.save(os.path.join(output_dir, 'case4_pareto_evolution.gif'), 
             writer=PillowWriter(fps=10), dpi=100)
    plt.close()
    
    print("  ✓ 案例4完成：功率分配器多目标优化")
    return True

# =============================================================================
# 案例5：多物理场灵敏度分析
# =============================================================================
def case5_sensitivity_analysis():
    """案例5：多物理场灵敏度分析"""
    print("\n案例5：多物理场灵敏度分析")
    print("-" * 60)
    
    # 定义多物理场响应函数
    def multiphysics_response(x):
        """多物理场响应：电磁、热、结构"""
        L, W, H, k = x  # 长度、宽度、高度、热导率
        
        # 电磁响应（谐振频率偏移）
        f_res = 2.4e9 * (1 - 0.01 * (L - 30) / 30 - 0.005 * (W - 20) / 20)
        
        # 热响应（最高温度）
        T_max = 80 - 5 * k + 2 * H
        
        # 结构响应（最大应力）
        sigma_max = 50 + 0.5 * L - 0.3 * W + 0.2 * H
        
        return np.array([f_res/1e9, T_max, sigma_max])
    
    # 基准设计点
    x0 = np.array([30.0, 20.0, 2.0, 200.0])  # L, W, H, k
    param_names = ['Length (mm)', 'Width (mm)', 'Height (mm)', 'Thermal Conductivity']
    response_names = ['Resonant Freq (GHz)', 'Max Temp (°C)', 'Max Stress (MPa)']
    
    # 局部灵敏度分析（有限差分法）
    print("  计算局部灵敏度...")
    dx = 0.01  # 1%扰动
    n_params = len(x0)
    n_responses = 3
    
    sensitivity_matrix = np.zeros((n_responses, n_params))
    
    for i in range(n_params):
        x_perturbed = x0.copy()
        x_perturbed[i] = x0[i] * (1 + dx)
        
        y_base = multiphysics_response(x0)
        y_perturbed = multiphysics_response(x_perturbed)
        
        sensitivity_matrix[:, i] = (y_perturbed - y_base) / (x0[i] * dx)
    
    # 全局灵敏度分析（Sobol方法简化版 - Morris筛选法）
    print("  计算全局灵敏度（Morris筛选法）...")
    
    np.random.seed(42)
    n_trajectories = 50
    delta = 0.1
    
    mu = np.zeros((n_responses, n_params))
    sigma = np.zeros((n_responses, n_params))
    
    for t in range(n_trajectories):
        # 随机起点
        x_base = np.array([np.random.uniform(25, 35), np.random.uniform(15, 25),
                          np.random.uniform(1.5, 2.5), np.random.uniform(150, 250)])
        
        for i in range(n_params):
            x_perturbed = x_base.copy()
            x_perturbed[i] += delta * x_base[i]
            
            y_base = multiphysics_response(x_base)
            y_perturbed = multiphysics_response(x_perturbed)
            
            ee = (y_perturbed - y_base) / delta  # 基本效应
            mu[:, i] += np.abs(ee)
            sigma[:, i] += ee**2
    
    mu /= n_trajectories
    sigma = np.sqrt(sigma / n_trajectories - mu**2)
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 灵敏度矩阵热力图
    ax1 = fig.add_subplot(2, 3, 1)
    im1 = ax1.imshow(sensitivity_matrix, cmap='RdBu_r', aspect='auto')
    ax1.set_xticks(range(n_params))
    ax1.set_yticks(range(n_responses))
    ax1.set_xticklabels(param_names, rotation=45, ha='right')
    ax1.set_yticklabels(response_names)
    ax1.set_title('Local Sensitivity Matrix', fontsize=12, fontweight='bold')
    plt.colorbar(im1, ax=ax1)
    
    # 添加数值标注
    for i in range(n_responses):
        for j in range(n_params):
            text = ax1.text(j, i, f'{sensitivity_matrix[i, j]:.2f}',
                           ha="center", va="center", color="black", fontsize=9)
    
    # Morris灵敏度 - 均值
    ax2 = fig.add_subplot(2, 3, 2)
    x_pos = np.arange(n_params)
    width = 0.25
    for i, resp_name in enumerate(response_names):
        ax2.bar(x_pos + i*width, mu[i], width, label=resp_name, alpha=0.8)
    ax2.set_xlabel('Design Parameter', fontsize=10)
    ax2.set_ylabel('Mean of Elementary Effects', fontsize=10)
    ax2.set_title('Morris Sensitivity (Mean)', fontsize=12, fontweight='bold')
    ax2.set_xticks(x_pos + width)
    ax2.set_xticklabels([p.split()[0] for p in param_names], rotation=45, ha='right')
    ax2.legend()
    ax2.grid(True, alpha=0.3, axis='y')
    
    # Morris灵敏度 - 标准差
    ax3 = fig.add_subplot(2, 3, 3)
    for i, resp_name in enumerate(response_names):
        ax3.bar(x_pos + i*width, sigma[i], width, label=resp_name, alpha=0.8)
    ax3.set_xlabel('Design Parameter', fontsize=10)
    ax3.set_ylabel('Std of Elementary Effects', fontsize=10)
    ax3.set_title('Morris Sensitivity (Std)', fontsize=12, fontweight='bold')
    ax3.set_xticks(x_pos + width)
    ax3.set_xticklabels([p.split()[0] for p in param_names], rotation=45, ha='right')
    ax3.legend()
    ax3.grid(True, alpha=0.3, axis='y')
    
    # Morris图（mu-sigma平面）
    ax4 = fig.add_subplot(2, 3, 4)
    colors = ['red', 'blue', 'green']
    markers = ['o', 's', '^']
    for i, (resp_name, color, marker) in enumerate(zip(response_names, colors, markers)):
        ax4.scatter(mu[i], sigma[i], c=color, marker=marker, s=100, label=resp_name, edgecolors='black')
        for j, param_name in enumerate(param_names):
            ax4.annotate(param_name.split()[0], (mu[i, j], sigma[i, j]), 
                        xytext=(5, 5), textcoords='offset points', fontsize=8)
    ax4.set_xlabel('Mean of Elementary Effects', fontsize=10)
    ax4.set_ylabel('Std of Elementary Effects', fontsize=10)
    ax4.set_title('Morris Plot (mu-sigma)', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 龙卷风图
    ax5 = fig.add_subplot(2, 3, 5)
    # 计算每个响应对各参数的灵敏度范围
    tornado_data = []
    for i in range(n_responses):
        low = sensitivity_matrix[i] - 0.5 * np.abs(sensitivity_matrix[i])
        high = sensitivity_matrix[i] + 0.5 * np.abs(sensitivity_matrix[i])
        for j in range(n_params):
            tornado_data.append({
                'Response': response_names[i],
                'Parameter': param_names[j],
                'Low': low[j],
                'High': high[j]
            })
    
    # 绘制龙卷风图（以第一个响应为例）
    y_pos = np.arange(n_params)
    low_vals = sensitivity_matrix[0] - 0.5 * np.abs(sensitivity_matrix[0])
    high_vals = sensitivity_matrix[0] + 0.5 * np.abs(sensitivity_matrix[0])
    ranges = high_vals - low_vals
    
    ax5.barh(y_pos, ranges, left=low_vals, color='steelblue', edgecolor='black', alpha=0.7)
    ax5.set_yticks(y_pos)
    ax5.set_yticklabels([p.split()[0] for p in param_names])
    ax5.set_xlabel('Sensitivity Range', fontsize=10)
    ax5.set_title(f'Tornado Chart: {response_names[0]}', fontsize=12, fontweight='bold')
    ax5.axvline(x=0, color='red', linestyle='--', linewidth=1)
    ax5.grid(True, alpha=0.3, axis='x')
    
    # 参数重要性排序
    ax6 = fig.add_subplot(2, 3, 6)
    importance = np.mean(np.abs(sensitivity_matrix), axis=0)
    sorted_idx = np.argsort(importance)[::-1]
    ax6.barh(range(n_params), importance[sorted_idx], color='coral', edgecolor='black')
    ax6.set_yticks(range(n_params))
    ax6.set_yticklabels([param_names[i].split()[0] for i in sorted_idx])
    ax6.set_xlabel('Average Absolute Sensitivity', fontsize=10)
    ax6.set_title('Parameter Importance Ranking', fontsize=12, fontweight='bold')
    ax6.grid(True, alpha=0.3, axis='x')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case5_sensitivity_analysis.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print("  灵敏度分析结果：")
    for i, resp in enumerate(response_names):
        print(f"    {resp}:")
        for j, param in enumerate(param_names):
            print(f"      {param}: {sensitivity_matrix[i, j]:.4f}")
    
    print("  ✓ 案例5完成：多物理场灵敏度分析")
    return True

# =============================================================================
# 案例6：代理模型辅助优化（含GIF动画）
# =============================================================================
def case6_surrogate_optimization():
    """案例6：代理模型辅助优化"""
    print("\n案例6：代理模型辅助优化")
    print("-" * 60)
    
    # 真实的高 fidelity 函数（模拟计算昂贵的多物理场仿真）
    def expensive_function(x):
        """昂贵的多物理场仿真"""
        x1, x2 = x
        # 复杂的非线性响应（模拟真实的多物理场行为）
        f = (x1 - 2)**2 + (x2 - 3)**2 + 0.5 * np.sin(3*x1) * np.cos(2*x2) + 0.1 * x1 * x2
        return f
    
    # Kriging代理模型（简化版）
    class KrigingSurrogate:
        def __init__(self):
            self.X = None
            self.y = None
            self.theta = 1.0
        
        def fit(self, X, y):
            self.X = np.array(X)
            self.y = np.array(y)
        
        def predict(self, X_new):
            X_new = np.atleast_2d(X_new)
            predictions = []
            for x in X_new:
                # 径向基函数插值
                distances = np.linalg.norm(self.X - x, axis=1)
                weights = np.exp(-self.theta * distances**2)
                if np.sum(weights) > 0:
                    pred = np.sum(weights * self.y) / np.sum(weights)
                else:
                    pred = np.mean(self.y)
                predictions.append(pred)
            return np.array(predictions)
    
    # 自适应采样优化
    print("  执行自适应代理模型优化...")
    
    np.random.seed(42)
    bounds = [(0, 5), (0, 5)]
    
    # 初始样本
    n_initial = 10
    X_train = np.random.rand(n_initial, 2) * 5
    y_train = np.array([expensive_function(x) for x in X_train])
    
    # 优化历史
    history = {
        'X': [X_train.copy()],
        'y': [y_train.copy()],
        'best': [y_train.min()],
        'n_evals': [n_initial]
    }
    
    n_iterations = 20
    
    for iteration in range(n_iterations):
        # 构建代理模型
        surrogate = KrigingSurrogate()
        surrogate.fit(X_train, y_train)
        
        # 在代理模型上优化（寻找EI最大的点）
        def expected_improvement(x):
            y_pred = surrogate.predict([x])[0]
            y_best = y_train.min()
            # 简化的EI计算
            if y_pred < y_best:
                return y_best - y_pred
            return 0
        
        # 网格搜索寻找下一个采样点
        grid_size = 50
        x1_grid = np.linspace(0, 5, grid_size)
        x2_grid = np.linspace(0, 5, grid_size)
        X1, X2 = np.meshgrid(x1_grid, x2_grid)
        
        ei_values = np.zeros_like(X1)
        for i in range(grid_size):
            for j in range(grid_size):
                ei_values[i, j] = expected_improvement([X1[i, j], X2[i, j]])
        
        # 选择EI最大的点（排除已有样本）
        max_ei = -1
        next_x = None
        for i in range(grid_size):
            for j in range(grid_size):
                x_candidate = np.array([X1[i, j], X2[i, j]])
                # 检查是否接近已有样本
                too_close = False
                for x_existing in X_train:
                    if np.linalg.norm(x_candidate - x_existing) < 0.1:
                        too_close = True
                        break
                if not too_close and ei_values[i, j] > max_ei:
                    max_ei = ei_values[i, j]
                    next_x = x_candidate
        
        if next_x is None:
            next_x = np.random.rand(2) * 5
        
        # 评估昂贵函数
        next_y = expensive_function(next_x)
        
        # 更新训练集
        X_train = np.vstack([X_train, next_x])
        y_train = np.append(y_train, next_y)
        
        # 记录历史
        history['X'].append(X_train.copy())
        history['y'].append(y_train.copy())
        history['best'].append(y_train.min())
        history['n_evals'].append(len(X_train))
    
    print(f"  总函数评估次数：{len(X_train)}")
    print(f"  最优值：{y_train.min():.4f}")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 10))
    
    # 真实函数曲面
    ax1 = fig.add_subplot(2, 3, 1, projection='3d')
    x1_fine = np.linspace(0, 5, 100)
    x2_fine = np.linspace(0, 5, 100)
    X1_fine, X2_fine = np.meshgrid(x1_fine, x2_fine)
    Z_fine = np.zeros_like(X1_fine)
    for i in range(100):
        for j in range(100):
            Z_fine[i, j] = expensive_function([X1_fine[i, j], X2_fine[i, j]])
    
    surf = ax1.plot_surface(X1_fine, X2_fine, Z_fine, cmap='viridis', alpha=0.7)
    ax1.scatter(X_train[:, 0], X_train[:, 1], y_train, c='red', s=50, edgecolors='black')
    ax1.set_xlabel('x1', fontsize=10)
    ax1.set_ylabel('x2', fontsize=10)
    ax1.set_zlabel('f(x)', fontsize=10)
    ax1.set_title('True Function with Samples', fontsize=12, fontweight='bold')
    
    # 代理模型预测曲面
    ax2 = fig.add_subplot(2, 3, 2, projection='3d')
    surrogate_final = KrigingSurrogate()
    surrogate_final.fit(X_train, y_train)
    Z_surrogate = np.zeros_like(X1_fine)
    for i in range(100):
        for j in range(100):
            Z_surrogate[i, j] = surrogate_final.predict([[X1_fine[i, j], X2_fine[i, j]]])[0]
    
    surf2 = ax2.plot_surface(X1_fine, X2_fine, Z_surrogate, cmap='plasma', alpha=0.7)
    ax2.scatter(X_train[:, 0], X_train[:, 1], y_train, c='red', s=50, edgecolors='black')
    ax2.set_xlabel('x1', fontsize=10)
    ax2.set_ylabel('x2', fontsize=10)
    ax2.set_zlabel('f(x)', fontsize=10)
    ax2.set_title('Surrogate Model Prediction', fontsize=12, fontweight='bold')
    
    # 采样点分布
    ax3 = fig.add_subplot(2, 3, 3)
    scatter = ax3.scatter(X_train[:, 0], X_train[:, 1], c=range(len(X_train)), 
                         cmap='coolwarm', s=100, edgecolors='black')
    ax3.plot(X_train[:, 0], X_train[:, 1], 'k--', alpha=0.3)
    ax3.set_xlabel('x1', fontsize=10)
    ax3.set_ylabel('x2', fontsize=10)
    ax3.set_title('Sampling Sequence', fontsize=12, fontweight='bold')
    plt.colorbar(scatter, ax=ax3, label='Sample Index')
    
    # 收敛曲线
    ax4 = fig.add_subplot(2, 3, 4)
    ax4.plot(history['n_evals'], history['best'], 'b-', linewidth=2, marker='o')
    ax4.set_xlabel('Function Evaluations', fontsize=10)
    ax4.set_ylabel('Best Objective Value', fontsize=10)
    ax4.set_title('Optimization Convergence', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    
    # 误差分析
    ax5 = fig.add_subplot(2, 3, 5)
    # 在测试点上评估误差
    n_test = 100
    X_test = np.random.rand(n_test, 2) * 5
    y_true = np.array([expensive_function(x) for x in X_test])
    y_pred = surrogate_final.predict(X_test)
    errors = np.abs(y_true - y_pred)
    
    ax5.scatter(y_true, y_pred, c='steelblue', alpha=0.6, edgecolors='black')
    ax5.plot([y_true.min(), y_true.max()], [y_true.min(), y_true.max()], 'r--', linewidth=2, label='Perfect')
    ax5.set_xlabel('True Value', fontsize=10)
    ax5.set_ylabel('Predicted Value', fontsize=10)
    ax5.set_title(f'Surrogate Accuracy (MAE={np.mean(errors):.3f})', fontsize=12, fontweight='bold')
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 误差直方图
    ax6 = fig.add_subplot(2, 3, 6)
    ax6.hist(errors, bins=20, color='coral', edgecolor='black', alpha=0.7)
    ax6.set_xlabel('Absolute Error', fontsize=10)
    ax6.set_ylabel('Frequency', fontsize=10)
    ax6.set_title('Prediction Error Distribution', fontsize=12, fontweight='bold')
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'case6_surrogate_optimization.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    # 创建GIF动画：代理模型优化过程
    print("  正在生成优化过程动画...")
    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 12))
    
    n_frames = min(30, len(history['X']))
    frame_indices = np.linspace(0, len(history['X'])-1, n_frames, dtype=int)
    
    def init():
        for ax in [ax1, ax2, ax3, ax4]:
            ax.clear()
        return []
    
    def update(frame_idx):
        for ax in [ax1, ax2, ax3, ax4]:
            ax.clear()
        
        idx = frame_indices[frame_idx]
        X_current = history['X'][idx]
        y_current = history['y'][idx]
        best_current = history['best'][idx]
        
        # 子图1: 当前样本分布
        ax1.scatter(X_current[:, 0], X_current[:, 1], c=y_current, cmap='viridis', 
                   s=100, edgecolors='black')
        ax1.set_xlabel('x1', fontsize=10)
        ax1.set_ylabel('x2', fontsize=10)
        ax1.set_title(f'Sample Distribution (n={len(X_current)})', fontsize=11, fontweight='bold')
        ax1.set_xlim(0, 5)
        ax1.set_ylim(0, 5)
        ax1.grid(True, alpha=0.3)
        
        # 子图2: 收敛曲线
        ax2.plot(history['n_evals'][:idx+1], history['best'][:idx+1], 'b-', linewidth=2, marker='o')
        ax2.set_xlabel('Function Evaluations', fontsize=10)
        ax2.set_ylabel('Best Objective', fontsize=10)
        ax2.set_title(f'Convergence (Best={best_current:.3f})', fontsize=11, fontweight='bold')
        ax2.grid(True, alpha=0.3)
        
        # 子图3: 代理模型预测
        if len(X_current) >= 5:
            surrogate_anim = KrigingSurrogate()
            surrogate_anim.fit(X_current, y_current)
            Z_anim = np.zeros_like(X1_fine)
            for i in range(50):
                for j in range(50):
                    Z_anim[i*2, j*2] = surrogate_anim.predict([[X1_fine[i*2, j*2], X2_fine[i*2, j*2]]])[0]
            
            im = ax3.contourf(X1_fine[::2, ::2], X2_fine[::2, ::2], Z_anim[::2, ::2], 
                             levels=15, cmap='plasma')
            ax3.scatter(X_current[:, 0], X_current[:, 1], c='red', s=50, edgecolors='black')
            ax3.set_xlabel('x1', fontsize=10)
            ax3.set_ylabel('x2', fontsize=10)
            ax3.set_title('Surrogate Model', fontsize=11, fontweight='bold')
        
        # 子图4: 进度
        progress = (frame_idx + 1) / n_frames
        ax4.barh(['Progress'], [progress*100], color='steelblue', edgecolor='black')
        ax4.set_xlim(0, 100)
        ax4.set_xlabel('Completion (%)', fontsize=10)
        ax4.set_title('Optimization Progress', fontsize=11, fontweight='bold')
        ax4.text(progress*100/2, 0, f'{progress*100:.0f}%', ha='center', va='center', 
                fontsize=12, fontweight='bold')
        
        plt.tight_layout()
        return []
    
    anim = FuncAnimation(fig, update, init_func=init, frames=n_frames, 
                        interval=150, blit=False, repeat=True)
    anim.save(os.path.join(output_dir, 'case6_surrogate_optimization.gif'), 
             writer=PillowWriter(fps=7), dpi=100)
    plt.close()
    
    print("  ✓ 案例6完成：代理模型辅助优化")
    return True


# =============================================================================
# 主程序
# =============================================================================
if __name__ == "__main__":
    results = {}
    
    # 运行所有案例
    results['case1'] = case1_power_amplifier_multiphysics()
    results['case2'] = case2_antenna_heatsink_design()
    results['case3'] = case3_connector_thermal_structural()
    results['case4'] = case4_power_divider_optimization()
    results['case5'] = case5_sensitivity_analysis()
    results['case6'] = case6_surrogate_optimization()
    
    # 汇总
    print("\n" + "=" * 70)
    print("仿真完成！")
    print("=" * 70)
    
    print("\n案例执行结果汇总：")
    for case_name, success in results.items():
        status = "成功" if success else "失败"
        print(f"  {case_name}: {status}")
    
    print(f"\n所有结果图片已保存到: {output_dir}")
    
    # 列出生成的文件
    print("\n生成的文件列表：")
    for file in sorted(os.listdir(output_dir)):
        if file.endswith(('.png', '.gif')):
            print(f"  - {file}")