主题099：多物理场优化设计

一、引言

1.1 多物理场优化设计的背景与意义

多物理场优化设计是现代工程仿真领域的前沿方向，它将多物理场耦合分析与优化算法相结合，实现复杂工程系统的自动设计。传统的设计方法往往采用"试错法"或基于经验的设计准则，难以在多个物理场耦合的复杂条件下找到最优设计方案。多物理场优化设计通过建立精确的数学模型，利用高效的优化算法，在满足各种约束条件的前提下，自动搜索最优设计参数，大大提高了设计效率和产品质量。

在航空航天领域，飞行器结构需要同时满足气动性能、热防护、结构强度等多重约束；在能源领域，核反应堆设计涉及中子物理、热工水力、结构力学等多个物理场的耦合；在电子领域，芯片封装需要考虑电磁兼容、热管理、机械可靠性等问题。这些复杂的工程问题都需要多物理场优化设计方法来解决。

1.2 多物理场优化设计的挑战

多物理场优化设计面临诸多挑战：

计算成本高：多物理场仿真本身计算量巨大，而优化过程需要反复调用仿真程序，导致总体计算成本极高。一个典型的多物理场优化问题可能需要数百甚至数千次完整的仿真计算。

强非线性：多物理场问题通常具有强非线性特征，包括材料非线性、几何非线性和边界非线性。这使得优化问题可能存在多个局部最优解，增加了全局寻优的难度。

多目标冲突：实际工程问题往往涉及多个相互冲突的目标，如重量与强度、效率与成本、性能与可靠性等。如何在Pareto前沿上找到满意的设计方案是一个复杂的多目标优化问题。

约束条件复杂：工程优化问题通常涉及大量约束条件，包括等式约束、不等式约束、显式约束和隐式约束（通过仿真计算获得）。处理这些约束需要专门的约束优化算法。

设计变量类型多样：设计变量可能是连续变量（如尺寸参数）、离散变量（如材料选择）或整数变量（如层数、单元数），这要求优化算法能够处理混合变量问题。

1.3 本主题内容概述

本主题将系统介绍多物理场优化设计的基本理论、数值方法和工程应用。主要内容包括：

• 优化理论基础：介绍优化问题的数学描述、分类和基本理论
• 单目标优化算法：详细讲解梯度法、遗传算法、粒子群优化等经典算法
• 多目标优化方法：介绍Pareto优化、NSGA-II等多目标优化技术
• 代理模型技术：讲解Kriging、神经网络等代理模型方法，用于降低计算成本
• 拓扑优化方法：介绍SIMP、水平集等结构拓扑优化方法
• 工程应用案例：通过6个典型案例展示多物理场优化设计的实际应用

二、优化理论基础

2.1 优化问题的数学描述

2.1.1 一般优化问题形式

多物理场优化问题可以统一表示为以下数学形式：

最小化: f(x)
约束条件:
    g_i(x) ≤ 0,  i = 1, 2, ..., m    (不等式约束)
    h_j(x) = 0,  j = 1, 2, ..., p    (等式约束)
    x_k^L ≤ x_k ≤ x_k^U, k = 1, 2, ..., n  (边界约束)

其中：

• x = [x₁, x₂, …, xₙ]ᵀ 是设计变量向量
• f(x) 是目标函数
• g_i(x) 是不等式约束函数
• h_j(x) 是等式约束函数
• x_k^L 和 x_k^U 是设计变量的上下界

2.1.2 多物理场优化问题的特点

在多物理场优化中，目标函数和约束函数通常需要通过数值仿真计算获得，无法表示为设计变量的显式函数。这带来了以下特点：

黑箱特性：仿真程序可以看作一个黑箱，输入设计变量，输出性能指标。优化算法只能通过输入输出关系来搜索最优解，无法直接获取梯度信息。

计算昂贵：每次仿真可能需要几分钟到几小时，甚至几天。因此，优化算法需要在有限的计算预算内找到尽可能好的解。

噪声影响：数值仿真结果可能存在数值误差，导致目标函数和约束函数存在噪声。优化算法需要具有一定的抗噪声能力。

隐式约束：某些约束条件（如收敛性约束）只能通过仿真过程隐式地体现，这增加了约束处理的难度。

2.2 优化问题分类

2.2.1 按目标函数数量分类

单目标优化：只有一个目标函数需要最小化或最大化。数学形式为：

min f(x)
s.t. x ∈ Ω

其中Ω是可行域。

多目标优化：需要同时优化多个相互冲突的目标。数学形式为：

min F(x) = [f₁(x), f₂(x), ..., fₖ(x)]
s.t. x ∈ Ω

多目标优化的结果是一组Pareto最优解，而非单一最优解。

2.2.2 按设计变量类型分类

连续优化：所有设计变量都是连续实数。这是最常见的优化问题类型，可以使用梯度-based方法和无梯度方法求解。

离散优化：设计变量只能取离散值。如材料选择问题，只能从有限的材料库中选择。

整数优化：设计变量必须是整数。如层合板的铺层数、加强筋的数量等。

混合整数优化：同时包含连续变量和整数/离散变量。这类问题最具挑战性，需要专门的混合整数优化算法。

2.2.3 按约束条件分类

无约束优化：没有任何约束条件，只需在定义域内寻找最优解。

等式约束优化：只包含等式约束。通常使用拉格朗日乘子法处理。

不等式约束优化：包含不等式约束。KKT条件是判断最优性的重要理论工具。

混合约束优化：同时包含等式和不等式约束。这是工程优化中最常见的情况。

2.3 最优性条件

2.3.1 无约束优化的最优性条件

对于无约束优化问题，局部最优解的必要条件是梯度为零：

∇f(x*) = 0

这称为驻点条件。但驻点可能是极小值点、极大值点或鞍点。

充分条件是Hessian矩阵正定：

∇²f(x*) ≻ 0

即Hessian矩阵的所有特征值都大于零。

2.3.2 约束优化的KKT条件

对于约束优化问题，Karush-Kuhn-Tucker (KKT) 条件是判断局部最优性的重要理论工具。

定义拉格朗日函数：

L(x, λ, μ) = f(x) + Σ λᵢ gᵢ(x) + Σ μⱼ hⱼ(x)

KKT条件包括：

1. 驻点条件：∇ₓL(x*, λ*, μ*) = 0
2. 原始可行性：gᵢ(x*) ≤ 0, hⱼ(x*) = 0
3. 对偶可行性：λᵢ ≥ 0
4. 互补松弛条件：λᵢ gᵢ(x*) = 0

对于凸优化问题，满足KKT条件的点一定是全局最优解。对于非凸问题，KKT条件只是局部最优的必要条件。

2.4 多物理场优化的特殊考虑

2.4.1 耦合约束处理

在多物理场优化中，不同物理场之间可能存在耦合约束。例如，在热-结构耦合问题中，温度场和位移场相互影响，需要同时满足热平衡方程和力平衡方程。

处理耦合约束的方法包括：

松耦合方法：将多物理场问题分解为多个单物理场子问题，通过迭代求解。优化时分别处理各物理场的约束。

紧耦合方法：将所有物理场的控制方程联立求解，形成统一的约束系统。这种方法精度高但计算成本大。

降阶模型方法：通过模型降阶技术减少耦合系统的自由度，在保证精度的前提下提高计算效率。

2.4.2 多尺度优化

多物理场问题往往涉及多个时间和空间尺度。例如，在复合材料优化中，微观纤维排列影响宏观力学性能。

多尺度优化策略：

分层优化：在不同尺度上分别建立优化模型，通过尺度间的信息传递实现协同优化。

并发优化：同时优化所有尺度的设计变量，建立多尺度耦合的优化模型。

代理模型方法：在细观尺度建立代理模型，将微观结构-性能关系快速映射到宏观优化中。

三、单目标优化算法

3.1 基于梯度的优化算法

3.1.1 梯度下降法

梯度下降法是最基本的优化算法，其迭代公式为：

x^{k+1} = x^k - αₖ ∇f(x^k)

其中αₖ是步长，可以通过线搜索确定。

算法特点：

• 实现简单，计算效率高
• 对于强凸问题收敛速度快
• 容易陷入局部最优
• 对于病态问题收敛慢

步长选择策略：

• 固定步长：简单但可能不稳定
• 精确线搜索：计算成本高但收敛性好
• 回溯线搜索：平衡计算成本和收敛性

3.1.2 牛顿法

牛顿法利用二阶信息（Hessian矩阵）来加速收敛：

x^{k+1} = x^k - [∇²f(x^k)]⁻¹ ∇f(x^k)

算法特点：

• 在极值点附近具有二次收敛速度
• 需要计算和存储Hessian矩阵
• 对于大规模问题计算成本高
• 可能收敛到鞍点

拟牛顿法：通过近似Hessian矩阵或其逆矩阵来降低计算成本，如BFGS算法、DFP算法等。

3.1.3 共轭梯度法

共轭梯度法是介于梯度下降法和牛顿法之间的算法，不需要存储Hessian矩阵，但收敛速度优于梯度下降法。

算法流程：

1. 初始化：x⁰，计算g⁰ = ∇f(x⁰)，设d⁰ = -g⁰
2. 线搜索：确定步长αₖ
3. 更新：x^{k+1} = x^k + αₖ d^k
4. 计算新梯度：g^{k+1} = ∇f(x^{k+1})
5. 计算共轭方向：d^{k+1} = -g^{k+1} + βₖ d^k
6. 重复步骤2-5直到收敛

应用范围：特别适合大规模稀疏优化问题，如有限元优化。

3.2 无梯度优化算法

3.2.1 遗传算法

遗传算法(Genetic Algorithm, GA)是一种模拟生物进化过程的启发式优化算法。

基本流程：

1. 初始化：随机生成初始种群
2. 适应度评估：计算每个个体的适应度（目标函数值）
3. 选择：根据适应度选择优秀个体进入下一代
4. 交叉：通过交叉操作产生新个体
5. 变异：对个体进行随机变异，增加种群多样性
6. 终止判断：如果满足终止条件则停止，否则返回步骤2

关键操作：

选择算子：

• 轮盘赌选择：按适应度比例选择
• 锦标赛选择：随机选取若干个体，选择最优者
• 排序选择：按适应度排序后选择

交叉算子：

• 单点交叉：随机选择交叉点，交换后半部分
• 多点交叉：多个交叉点进行基因交换
• 均匀交叉：每个基因独立决定是否交换

变异算子：

• 位变异：对二进制编码的个体翻转某些位
• 高斯变异：对实数编码的个体添加高斯噪声

算法参数：

• 种群规模：通常50-200
• 交叉概率：0.6-0.9
• 变异概率：0.01-0.1
• 终止代数：100-1000

优点：

• 不依赖梯度信息，适用于黑箱优化
• 具有全局搜索能力，不易陷入局部最优
• 可以处理离散变量和混合变量

缺点：

• 计算成本高，需要大量函数评估
• 参数设置敏感，需要经验调整
• 收敛速度慢，特别是后期收敛

3.2.2 粒子群优化

粒子群优化(Particle Swarm Optimization, PSO)模拟鸟群觅食行为，通过粒子间的协作寻找最优解。

算法原理：

每个粒子具有位置和速度两个属性，在搜索空间中飞行。粒子根据以下信息更新自己的速度和位置：

• 自身历史最优位置（个体经验）
• 种群历史最优位置（群体经验）

速度更新公式：

vᵢ^{k+1} = w·vᵢ^k + c₁·r₁·(pbestᵢ - xᵢ^k) + c₂·r₂·(gbest - xᵢ^k)

位置更新公式：

xᵢ^{k+1} = xᵢ^k + vᵢ^{k+1}

其中：

• w：惯性权重，控制历史速度的影响
• c₁, c₂：学习因子，控制个体和群体经验的影响
• r₁, r₂：[0,1]范围内的随机数

算法流程：

1. 初始化粒子群的位置和速度
2. 评估每个粒子的适应度
3. 更新个体最优和全局最优
4. 按公式更新粒子速度和位置
5. 检查终止条件，未满足则返回步骤2

参数设置：

• 种群规模：20-50
• 惯性权重：0.4-0.9（可线性递减）
• 学习因子：c₁ = c₂ = 2.0（常用设置）
• 最大速度：通常设为搜索范围的10%-20%

改进策略：

惯性权重调整：

• 线性递减：w从0.9递减到0.4
• 自适应调整：根据收敛情况动态调整

拓扑结构：

• 全局PSO：所有粒子共享全局最优
• 局部PSO：粒子只与邻域内粒子交流

混合策略：

• 与模拟退火结合，增强跳出局部最优的能力
• 与差分进化结合，提高搜索效率

3.2.3 差分进化算法

差分进化(Differential Evolution, DE)是一种基于向量差分的进化算法，特别适合连续优化问题。

变异操作：

vᵢ = x_{r1} + F·(x_{r2} - x_{r3})

其中r1, r2, r3是随机选择的互不相同的个体，F是缩放因子。

交叉操作：

uᵢⱼ = { vᵢⱼ  if rand(0,1) ≤ CR or j = j_rand
       { xᵢⱼ  otherwise

CR是交叉概率，j_rand确保至少有一个维度来自变异向量。

选择操作：

xᵢ^{new} = { uᵢ  if f(uᵢ) ≤ f(xᵢ)
           { xᵢ  otherwise

采用贪婪选择，只有更好的个体才能进入下一代。

算法特点：

• 结构简单，易于实现
• 收敛速度快，特别是前期
• 控制参数少，易于调整
• 在连续优化问题上表现优异

参数设置：

• 种群规模：5D-10D（D为维度）
• 缩放因子F：0.5-1.0
• 交叉概率CR：0.1-1.0

变异策略：

• DE/rand/1：随机选择基向量
• DE/best/1：选择当前最优作为基向量
• DE/current-to-best/1：结合当前个体和最优个体
• DE/best/2：使用两个差分向量

3.3 约束处理技术

3.3.1 罚函数法

罚函数法通过将约束 violation 加入目标函数，将有约束问题转化为无约束问题。

外点罚函数：

P(x, r) = f(x) + r·Σ[max(0, gᵢ(x))]² + r·Σ[hⱼ(x)]²

当r→∞时，罚函数的最优解收敛于原问题的最优解。

内点罚函数：

P(x, r) = f(x) - r·Σ 1/gᵢ(x)

要求初始点在可行域内，随着r→0收敛于最优解。

增广拉格朗日法：

L(x, λ, r) = f(x) + Σ λᵢ·gᵢ(x) + r·Σ[max(0, gᵢ(x))]²

结合拉格朗日乘子和罚函数，收敛性更好。

3.3.2 可行性规则

在进化算法中，常用的约束处理方法是可行性规则：

1. 可行解总是优于不可行解
2. 两个可行解比较，目标函数值更优者胜出
3. 两个不可行解比较，约束 violation 更小者胜出

优点：

• 简单直观，易于实现
• 不需要设置罚因子
• 能较好地平衡可行性和最优性

缺点：

• 可能导致种群多样性丧失
• 对于高度约束问题效果不佳

3.3.3 多目标化方法

将约束 violation 作为额外的目标函数，将约束优化转化为多目标优化：

min [f(x), Σ max(0, gᵢ(x)), Σ |hⱼ(x)|]

优点：

• 避免罚因子的设置问题
• 可以获得一系列不同权衡的解
• 有助于保持种群多样性

缺点：

• 增加了问题的复杂度
• 需要多目标优化算法

四、多目标优化方法

4.1 Pareto最优性理论

4.1.1 Pareto支配关系

在多目标优化中，由于目标之间可能存在冲突，不存在一个解能在所有目标上都最优。因此引入Pareto支配的概念：

定义：解x₁ Pareto支配解x₂（记作x₁ ≺ x₂），当且仅当：

• 对所有目标，x₁不劣于x₂：fᵢ(x₁) ≤ fᵢ(x₂), ∀i
• 至少对一个目标，x₁严格优于x₂：∃j, fⱼ(x₁) < fⱼ(x₂)

Pareto最优解：如果不存在其他解能够支配它，则该解称为Pareto最优解。

Pareto前沿：所有Pareto最优解在目标空间中的映射形成的曲面（或曲线）。

4.1.2 Pareto前沿特性

理想点：各目标单独优化时的最优值组成的点，通常不可行。

z* = [min f₁(x), min f₂(x), ..., min fₖ(x)]

最低点：各目标单独优化时的最差值组成的点。

z^{nad} = [max f₁(x), max f₂(x), ..., max fₖ(x)]

（在Pareto前沿上取最大值）

乌托邦点：理想点的一个近似，用于归一化目标函数。

Pareto前沿形状：

• 凸前沿：加权和法可以找到所有Pareto最优解
• 非凸前沿：需要ε-约束法或其他特殊方法
• 不连续前沿：某些区域没有Pareto最优解

4.2 经典多目标优化方法

4.2.1 加权和法

将多目标问题转化为单目标问题：

min Σ wᵢ·fᵢ(x)
s.t. Σ wᵢ = 1, wᵢ ≥ 0

优点：

• 简单直观，可以使用任何单目标优化算法
• 对于凸Pareto前沿，通过调整权重可以获得所有Pareto最优解

缺点：

• 对于非凸前沿，无法获得某些区域的解
• 权重设置困难，缺乏物理意义
• 不同目标量级差异大时需要归一化

归一化处理：

fᵢ^{norm}(x) = [fᵢ(x) - fᵢ^{min}] / [fᵢ^{max} - fᵢ^{min}]

4.2.2 ε-约束法

保留一个主要目标，将其他目标转化为约束：

min fₖ(x)
s.t. fᵢ(x) ≤ εᵢ, i ≠ k

优点：

• 可以处理非凸Pareto前沿
• 决策者可以控制各目标的满意水平
• 获得的解都是Pareto最优的

缺点：

• 需要合理设置ε值
• 约束问题可能无解
• 需要求解多个单目标问题

ε值选择：

• 通过网格搜索系统性地改变ε值
• 根据Pareto前沿的密度自适应调整

4.2.3 目标规划法

为各目标设定目标值，最小化与目标的偏差：

min Σ |fᵢ(x) - tᵢ|

或

min Σ wᵢ·|fᵢ(x) - tᵢ|

优点：

• 符合工程决策习惯
• 可以处理目标优先级
• 容易理解，便于与决策者交互

缺点：

• 目标值设置困难
• 可能得到非Pareto最优解
• 对于复杂问题可能无解

4.3 进化多目标优化算法

4.3.1 NSGA-II算法

非支配排序遗传算法II (Non-dominated Sorting Genetic Algorithm II) 是最流行的多目标进化算法之一。

关键机制：

快速非支配排序：

• 将种群按支配关系分层
• 第一层：非支配解（Pareto前沿）
• 第二层：只被第一层支配的解
• 依此类推…

拥挤度计算：

• 对同一层的解计算拥挤度
• 拥挤度反映解周围个体的密度
• 优先选择拥挤度大的解，保持多样性

精英保留策略：

• 将父代和子代合并
• 按非支配排序和拥挤度选择下一代
• 保证优秀解不被丢失

算法流程：

1. 初始化种群P₀，计算目标函数值
2. 对P₀进行非支配排序
3. 使用遗传算子（选择、交叉、变异）产生子代Q₀
4. 合并父代和子代：Rₜ = Pₜ ∪ Qₜ
5. 对Rₜ进行非支配排序
6. 按排序层级和拥挤度选择N个个体组成P_{t+1}
7. 产生新的子代Q_{t+1}
8. 重复步骤4-7直到满足终止条件

优点：

• 收敛性好，能快速逼近真实Pareto前沿
• 多样性保持机制有效
• 计算复杂度低：O(MN²)，M为目标数，N为种群规模

缺点：

• 对于高维目标问题（M>3）效果下降
• 需要较大的种群规模

4.3.2 MOEA/D算法

基于分解的多目标进化算法 (Multi-Objective Evolutionary Algorithm based on Decomposition) 将多目标问题分解为多个单目标子问题，同时优化。

分解方法：

加权和法：

g^{ws}(x|λ) = Σ λᵢ·fᵢ(x)

Tchebycheff法：

g^{te}(x|λ, z*) = max λᵢ·|fᵢ(x) - z*ᵢ|

基于惩罚的边界交叉法(PBI)：

g^{pbi}(x|λ, z*) = d₁ + θ·d₂

其中d₁是到权重向量的距离，d₂是垂直距离。

邻域关系：

• 每个子问题对应一个权重向量
• 距离近的权重向量定义的子问题为邻域
• 子问题只在邻域内协作，降低计算成本

算法流程：

1. 初始化权重向量和邻域结构
2. 初始化种群，计算目标函数值
3. 对每个子问题i：
   • 从邻域内随机选择两个父代
   • 使用遗传算子产生新解y
   • 更新邻域内的解：如果g(y|λʲ) < g(xʲ|λʲ)，则用y替换xʲ
4. 重复步骤3直到满足终止条件

优点：

• 计算效率高，适合高维目标问题
• 可以控制Pareto前沿上解的分布
• 子问题间信息共享，加速收敛

缺点：

• 权重向量设计影响性能
• 对于非凸前沿效果不佳

4.3.3 性能评价指标

收敛性指标：

世代距离(GD)：

GD = √[Σ dᵢ²] / n

其中dᵢ是第i个解到真实Pareto前沿的最小距离。

反转世代距离(IGD)：

IGD = √[Σ d'ᵢ²] / n'

其中d’ᵢ是真实Pareto前沿上的点到获得的前沿的最小距离。IGD同时评价收敛性和多样性。

多样性指标：

Spacing(SP)：

SP = √[Σ (d̄ - dᵢ)² / n]

其中dᵢ是第i个解到其最近邻的距离，d̄是平均距离。

超体积(HV)：

HV = volume(∪ Sᵢ)

其中Sᵢ是以参考点和第i个解为对角顶点的超立方体。HV同时评价收敛性和多样性，值越大越好。

五、代理模型技术

5.1 代理模型概述

5.1.1 为什么需要代理模型

多物理场仿真计算成本高昂，直接使用仿真程序进行优化需要大量计算资源。代理模型（Surrogate Model）或元模型（Metamodel）通过近似原始仿真模型，在保持一定精度的前提下大幅降低计算成本。

应用场景：

• 优化设计：用代理模型替代昂贵的仿真
• 不确定性量化：快速进行蒙特卡洛模拟
• 可视化：生成响应面、等高线图
• 敏感性分析：快速计算灵敏度

5.1.2 代理模型构建流程

步骤1：试验设计(Design of Experiments, DOE)

• 确定采样点和样本数量
• 常用方法：拉丁超立方采样(LHS)、正交设计、均匀设计

步骤2：样本点仿真

• 在采样点处运行高精度仿真
• 获得输入-输出数据对

步骤3：模型拟合

• 选择代理模型类型
• 用样本数据训练模型参数

步骤4：模型验证

• 用额外的验证点评估模型精度
• 常用指标：R²、RMSE、MAE

步骤5：模型应用

• 在优化、UQ等任务中使用代理模型
• 必要时进行模型更新（自适应代理模型）

5.2 常用代理模型

5.2.1 多项式响应面法(RSM)

用多项式函数近似响应：

y = β₀ + Σ βᵢxᵢ + Σ βᵢᵢxᵢ² + Σ Σ βᵢⱼxᵢxⱼ + ε

线性模型：

y = β₀ + Σ βᵢxᵢ

二次模型：

y = β₀ + Σ βᵢxᵢ + Σ βᵢᵢxᵢ² + Σ Σ βᵢⱼxᵢxⱼ

参数估计：
使用最小二乘法估计系数：

β = (XᵀX)⁻¹ Xᵀy

优点：

• 形式简单，易于理解和解释
• 计算成本低，适合快速近似
• 可以分析各因素的影响

缺点：

• 对于高度非线性问题精度低
• 高维问题需要大量样本
• 外推能力差

5.2.2 Kriging模型

Kriging（高斯过程回归）是一种基于统计理论的插值方法，不仅给出预测值，还给出预测不确定性。

模型形式：

y(x) = f(x)ᵀβ + z(x)

其中f(x)ᵀβ是趋势函数（通常为常数或低阶多项式），z(x)是零均值高斯过程。

相关函数：

Corr[z(xⁱ), z(xʲ)] = R(θ, xⁱ, xʲ)

常用相关函数：

• 高斯相关：R = exp(-θ·d²)
• 指数相关：R = exp(-θ·|d|)
• Matérn相关：更灵活的形式

预测公式：

ŷ(x) = f(x)ᵀβ̂ + r(x)ᵀ R⁻¹(y - Fβ̂)

预测方差：

s²(x) = σ²[1 - r(x)ᵀ R⁻¹ r(x) + u(x)ᵀ(Fᵀ R⁻¹ F)⁻¹ u(x)]

优点：

• 插值精度高，通过所有样本点
• 提供预测不确定性估计
• 适合高度非线性问题
• 可以指导自适应采样

缺点：

• 计算成本高，需要矩阵求逆
• 超参数估计复杂
• 高维问题效果下降

5.2.3 径向基函数(RBF)

RBF用基函数的线性组合近似响应：

y(x) = Σ wᵢ·φ(||x - xᵢ||) + Σ βⱼ·pⱼ(x)

常用基函数：

• 高斯函数：φ® = exp(-εr²)
• 多二次函数：φ® = √(1 + (εr)²)
• 逆多二次函数：φ® = 1/√(1 + (εr)²)
• 薄板样条：φ® = r² log®

参数确定：
通过求解线性方程组确定权重wᵢ：

[Φ  P][w]   [y]
[Pᵀ 0][β] = [0]

优点：

• 形式灵活，可以逼近任意连续函数
• 插值精度高
• 计算效率较高

缺点：

• 形状参数选择困难
• 对于大规模问题矩阵求逆困难
• 外推能力有限

5.2.4 神经网络

人工神经网络，特别是深度学习模型，可以学习复杂的非线性映射关系。

多层感知机(MLP)：

y = W₂·σ(W₁·x + b₁) + b₂

激活函数：

• ReLU：σ(x) = max(0, x)
• Sigmoid：σ(x) = 1/(1 + exp(-x))
• Tanh：σ(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x))

训练过程：

1. 前向传播：计算网络输出
2. 计算损失：如MSE = Σ(y - ŷ)²/n
3. 反向传播：计算梯度
4. 参数更新：使用梯度下降或其变种

优点：

• 表达能力极强，可以逼近任意函数
• 适合高维问题
• 可以处理大规模数据

缺点：

• 需要大量训练数据
• 训练时间长，需要调参
• 黑箱特性，可解释性差
• 可能过拟合

5.3 自适应代理模型

5.3.1 基于期望改进的加点准则

在优化过程中，自适应地选择新的采样点，平衡探索和利用。

期望改进(EI)：

EI(x) = E[max(y_{min} - Y(x), 0)]

对于Kriging模型：

EI(x) = (y_{min} - ŷ(x))·Φ(Z) + s(x)·φ(Z)

其中Z = (y_{min} - ŷ(x))/s(x)，Φ和φ是标准正态分布的CDF和PDF。

加点策略：

1. 在当前代理模型上找到EI最大的点
2. 在该点运行高精度仿真
3. 更新代理模型
4. 重复直到收敛或预算耗尽

5.3.2 多精度代理模型

结合不同精度的仿真模型，在保证精度的前提下降低成本。

协同Kriging(Co-Kriging)：

• 同时建模高精度(HF)和低精度(LF)数据
• LF数据提供趋势信息
• HF数据修正局部误差

分层Kriging：

• 先建立LF模型
• 将LF模型作为趋势函数
• 用HF数据拟合残差

空间映射：

• 建立LF和HF之间的映射关系
• 通过映射校正LF预测

六、拓扑优化方法

6.1 拓扑优化概述

6.1.1 结构优化的三个层次

尺寸优化：在给定的结构形状和拓扑下，优化截面尺寸参数。如杆件的截面积、板的厚度等。

形状优化：在给定的拓扑结构下，优化结构的边界形状。如孔洞的形状、圆角的半径等。

拓扑优化：在给定的设计域内，确定材料的最优分布。可以创造新的孔洞、改变连接关系，设计空间最大。

6.1.2 拓扑优化的数学描述

一般形式：

min C(ρ) = FᵀU
s.t. K(ρ)U = F
     V(ρ) = Σ ρₑvₑ ≤ V̄
     0 < ρ_min ≤ ρₑ ≤ 1

其中：

• ρₑ是单元e的密度（设计变量）
• C是柔度（目标函数，最小化相当于最大化刚度）
• K是总体刚度矩阵
• F是载荷向量
• U是位移向量
• V是材料体积
• V̄是允许的最大体积

6.2 SIMP方法

6.2.1 密度-刚度插值

固体各向同性材料惩罚法(Solid Isotropic Material with Penalization, SIMP)是最常用的拓扑优化方法。

材料插值模型：

Eₑ = E_min + ρₑ^p (E₀ - E_min)

其中：

• Eₑ是单元e的弹性模量
• E₀是实体材料的弹性模量
• E_min是很小的数，避免奇异
• ρₑ ∈ [0,1]是单元密度
• p是惩罚因子（通常p=3）

惩罚机制：

• 当p>1时，中间密度(0<ρ<1)的材料效率低于线性插值
• 这驱使优化结果趋向0-1分布（实体或空洞）
• p越大，惩罚越强，但可能导致收敛困难

6.2.2 灵敏度分析

目标函数灵敏度：

∂C/∂ρₑ = -Uₑᵀ (∂Kₑ/∂ρₑ) Uₑ = -p·ρₑ^{p-1} (E₀ - E_min) Uₑᵀ K₀ₑ Uₑ

体积约束灵敏度：

∂V/∂ρₑ = vₑ

伴随法：
通过求解伴随方程，可以高效计算梯度：

Kλ = -∂C/∂U
∂C/∂ρₑ = ∂C/∂ρₑ + λᵀ (∂K/∂ρₑ U - ∂F/∂ρₑ)

6.2.3 数值稳定性技术

棋盘格现象：
优化结果出现材料与空洞交替的棋盘格模式，这是数值不稳定性导致的。

解决方法：

滤波技术：

ρ̃ₑ = Σ w(rₑᵢ)ρᵢ / Σ w(rₑᵢ)

其中w是权重函数，rₑᵢ是单元e和i之间的距离。

常用滤波函数：

• 线性滤波：w® = max(0, r_min - r)
• 高斯滤波：w® = exp(-r²/(2σ²))

投影方法：

ρ̄ₑ = tanh(βη) + tanh(β(ρ̃ₑ - η)) / (tanh(βη) + tanh(β(1-η)))

通过Heaviside投影获得清晰的0-1分布。

周长约束：
限制结构边界的总长度，避免过于复杂的边界。

6.3 水平集方法

6.3.1 水平集函数

水平集方法用隐式函数描述结构边界：

Ω = {x | φ(x) < 0}      (材料区域)
∂Ω = {x | φ(x) = 0}     (边界)
D\Ω = {x | φ(x) > 0}    (空洞区域)

符号距离函数：

|∇φ| = 1

保持水平集函数为符号距离函数有利于数值稳定性。

6.3.2 水平集演化

通过求解Hamilton-Jacobi方程演化水平集函数：

∂φ/∂t + Vₙ|∇φ| = 0

其中Vₙ是法向速度，通常取为目标函数的负灵敏度：

Vₙ = -G(x) = -(-∂C/∂Ω + λ)

速度扩展：
只在边界上定义速度，需要扩展到整个设计域。

6.3.3 SIMP vs 水平集

特性	SIMP	水平集
设计变量	单元密度	水平集函数
边界描述	模糊的	清晰的
拓扑变化	容易（孔洞可生灭）	困难（需重新初始化）
几何特征	可能出现孤岛	光滑连续
多物理场	容易扩展	较复杂

6.4 多物理场拓扑优化

6.4.1 热-结构耦合拓扑优化

同时考虑热传导和结构力学性能：

目标函数：

min C = w₁·C_thermal + w₂·C_mechanical

或

min C_mechanical
s.t. T_max ≤ T̄

耦合关系：

• 温度场影响材料属性（如弹性模量随温度变化）
• 热应力是结构载荷的一部分
• 材料分布同时影响热传导和结构刚度

6.4.2 流-固耦合拓扑优化

优化流体通道和固体结构的布局：

Navier-Stokes方程：

ρ(u·∇)u = -∇p + μ∇²u + f
∇·u = 0

Darcy定律近似：
在多孔介质区域，用Brinkman模型近似：

-μ∇²u + α(ρ)u + ∇p = f

其中α是渗透率，与密度相关。

拓扑优化模型：

min Δp 或 max Q
s.t. 流动方程
     体积约束

6.4.3 电磁-热耦合拓扑优化

优化电磁器件的散热结构：

Maxwell方程：

∇×E = -∂B/∂t
∇×H = J + ∂D/∂t

热生成：

q = J·E = σ|E|²  (焦耳热)

优化目标：

• 最小化热点温度
• 最大化散热效率
• 最小化电磁损耗

七、工程应用案例

7.1 案例1：热沉结构拓扑优化

7.1.1 问题描述

设计一个电子设备的热沉，在给定体积约束下最大化散热效率。

物理模型：

• 稳态热传导方程：∇·(k∇T) + q = 0
• 对流边界：-k∇T·n = h(T - T∞)

设计域：100mm × 60mm × 5mm的矩形区域

边界条件：

• 底部：固定热流输入 q = 10⁵ W/m²
• 其他表面：对流换热 h = 50 W/(m²·K), T∞ = 25°C

优化目标：最小化平均温度
约束条件：材料体积分数 ≤ 30%

7.1.2 Python实现

"""
案例1：热沉结构拓扑优化
使用SIMP方法进行热传导拓扑优化
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import LinearSegmentedColormap

print("="*60)
print("案例1：热沉结构拓扑优化")
print("="*60)

# 设计域参数
nelx, nely = 100, 60  # 网格数
Lx, Ly = 100e-3, 60e-3  # 尺寸 (m)
dx, dy = Lx/nelx, Ly/nely

# 材料参数
k0 = 200  # 固体导热系数 (W/(m·K))
k_min = 1e-3  # 空洞等效导热系数
h_conv = 50  # 对流换热系数 (W/(m²·K))
T_inf = 25  # 环境温度 (°C)
q_in = 1e5  # 热流输入 (W/m²)

# 优化参数
volfrac = 0.3  # 体积分数约束
p = 3  # 惩罚因子
rmin = 3  # 滤波半径
max_iter = 100  # 最大迭代数

# 初始化设计变量
rho = np.ones((nely, nelx)) * volfrac

# 热传导矩阵组装函数
def assemble_K(rho):
    """组装全局热传导矩阵"""
    n_dof = (nelx+1) * (nely+1)
    K = np.zeros((n_dof, n_dof))
    F = np.zeros(n_dof)
    
    # 单元热传导矩阵（简化的一阶近似）
    ke = k0 * np.array([[ 2/3, -1/6, -1/3, -1/6],
                        [-1/6,  2/3, -1/6, -1/3],
                        [-1/3, -1/6,  2/3, -1/6],
                        [-1/6, -1/3, -1/6,  2/3]])
    
    for ely in range(nely):
        for elx in range(nelx):
            # 单元节点编号
            n1 = ely * (nelx+1) + elx
            n2 = n1 + 1
            n3 = n1 + nelx + 1
            n4 = n3 + 1
            nodes = [n1, n2, n4, n3]
            
            # 插值导热系数
            k_e = k_min + rho[ely, elx]**p * (k0 - k_min)
            
            # 组装
            for i in range(4):
                for j in range(4):
                    K[nodes[i], nodes[j]] += k_e * ke[i, j] / k0
    
    # 对流边界条件（上表面）
    for elx in range(nelx):
        node1 = nely * (nelx+1) + elx
        node2 = node1 + 1
        K[node1, node1] += h_conv * dx / 2
        K[node2, node2] += h_conv * dx / 2
        F[node1] += h_conv * dx / 2 * T_inf
        F[node2] += h_conv * dx / 2 * T_inf
    
    # 热流输入（底面中心区域）
    center_start = nelx // 3
    center_end = 2 * nelx // 3
    for elx in range(center_start, center_end):
        node1 = elx
        node2 = elx + 1
        F[node1] += q_in * dx / 2
        F[node2] += q_in * dx / 2
    
    return K, F

# 求解温度场
def solve_temperature(K, F):
    """求解温度场"""
    n_dof = len(F)
    
    # 处理边界条件（固定边界温度）
    fixed_dofs = []
    # 四个角点固定温度
    fixed_dofs.extend([0, nelx, nely*(nelx+1), (nely+1)*(nelx+1)-1])
    
    free_dofs = [i for i in range(n_dof) if i not in fixed_dofs]
    
    # 求解
    T = np.ones(n_dof) * T_inf
    K_reduced = K[np.ix_(free_dofs, free_dofs)]
    F_reduced = F[free_dofs]
    
    try:
        T_reduced = np.linalg.solve(K_reduced, F_reduced)
        T[free_dofs] = T_reduced
    except:
        print("求解失败，使用迭代法")
        T_reduced = np.linalg.lstsq(K_reduced, F_reduced, rcond=None)[0]
        T[free_dofs] = T_reduced
    
    return T

# 灵敏度分析
def compute_sensitivity(T, rho):
    """计算目标函数对设计变量的灵敏度"""
    sensitivity = np.zeros((nely, nelx))
    
    for ely in range(nely):
        for elx in range(nelx):
            # 单元节点
            n1 = ely * (nelx+1) + elx
            n2 = n1 + 1
            n3 = n1 + nelx + 1
            n4 = n3 + 1
            
            T_e = np.array([T[n1], T[n2], T[n4], T[n3]])
            
            # 简化灵敏度计算
            ke = np.array([[ 2/3, -1/6, -1/3, -1/6],
                          [-1/6,  2/3, -1/6, -1/3],
                          [-1/3, -1/6,  2/3, -1/6],
                          [-1/6, -1/3, -1/6,  2/3]])
            
            dk_drho = p * rho[ely, elx]**(p-1) * (k0 - k_min)
            sensitivity[ely, elx] = -dk_drho * np.dot(T_e, np.dot(ke, T_e)) / k0
    
    return sensitivity

# 密度滤波
def density_filter(rho, sensitivity):
    """应用密度滤波"""
    rho_filtered = np.zeros_like(rho)
    sensitivity_filtered = np.zeros_like(sensitivity)
    
    for ely in range(nely):
        for elx in range(nelx):
            sum_weight = 0
            weighted_rho = 0
            weighted_sens = 0
            
            for dy in range(-int(rmin), int(rmin)+1):
                for dx in range(-int(rmin), int(rmin)+1):
                    ny, nx = ely + dy, elx + dx
                    if 0 <= ny < nely and 0 <= nx < nelx:
                        dist = np.sqrt(dx**2 + dy**2)
                        if dist < rmin:
                            weight = rmin - dist
                            sum_weight += weight
                            weighted_rho += weight * rho[ny, nx]
                            weighted_sens += weight * sensitivity[ny, nx]
            
            rho_filtered[ely, elx] = weighted_rho / sum_weight
            sensitivity_filtered[ely, elx] = weighted_sens / sum_weight
    
    return rho_filtered, sensitivity_filtered

# 优化循环
print("\n开始优化...")
history = {'obj': [], 'vol': []}

for iter in range(max_iter):
    # 有限元分析
    K, F = assemble_K(rho)
    T = solve_temperature(K, F)
    
    # 计算目标函数（平均温度）
    obj = np.mean(T)
    vol = np.mean(rho)
    
    history['obj'].append(obj)
    history['vol'].append(vol)
    
    if iter % 10 == 0:
        print(f"迭代 {iter}: 目标 = {obj:.2f}°C, 体积分数 = {vol:.3f}")
    
    # 灵敏度分析
    sensitivity = compute_sensitivity(T, rho)
    
    # 密度滤波
    rho_filtered, sensitivity_filtered = density_filter(rho, sensitivity)
    
    # 更新设计变量（OC方法）
    move = 0.2
    rho_new = np.zeros_like(rho)
    
    # 二分法寻找拉格朗日乘子
    l1, l2 = 0, 1e6
    while (l2 - l1) / (l1 + l2) > 1e-4:
        lmid = 0.5 * (l1 + l2)
        rho_new = np.maximum(0.001, np.maximum(rho - move, 
                            np.minimum(1.0, np.minimum(rho + move, 
                            rho * np.sqrt(-sensitivity_filtered / lmid)))))
        
        if np.mean(rho_new) > volfrac:
            l1 = lmid
        else:
            l2 = lmid
    
    rho = rho_new

print("\n优化完成!")
print(f"最终目标函数值: {history['obj'][-1]:.2f}°C")
print(f"最终体积分数: {history['vol'][-1]:.3f}")

# 可视化结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 优化后的结构
ax1 = axes[0, 0]
im1 = ax1.imshow(rho, cmap='gray_r', origin='lower', extent=[0, Lx*1000, 0, Ly*1000])
ax1.set_xlabel('X (mm)')
ax1.set_ylabel('Y (mm)')
ax1.set_title('Optimized Structure')
plt.colorbar(im1, ax=ax1, label='Density')

# 温度场
T_grid = T.reshape((nely+1, nelx+1))
ax2 = axes[0, 1]
im2 = ax2.imshow(T_grid, cmap='hot', origin='lower', extent=[0, Lx*1000, 0, Ly*1000])
ax2.set_xlabel('X (mm)')
ax2.set_ylabel('Y (mm)')
ax2.set_title('Temperature Field')
plt.colorbar(im2, ax=ax2, label='Temperature (°C)')

# 收敛历史
ax3 = axes[1, 0]
ax3.plot(history['obj'], 'b-', linewidth=2)
ax3.set_xlabel('Iteration')
ax3.set_ylabel('Average Temperature (°C)')
ax3.set_title('Objective Convergence')
ax3.grid(True, alpha=0.3)

# 体积分数历史
ax4 = axes[1, 1]
ax4.plot(history['vol'], 'r-', linewidth=2)
ax4.axhline(y=volfrac, color='k', linestyle='--', label='Target')
ax4.set_xlabel('Iteration')
ax4.set_ylabel('Volume Fraction')
ax4.set_title('Volume Constraint')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case1_heatsink_topology.png', dpi=150)
print("\n结果已保存到 case1_heatsink_topology.png")

7.2 案例2：电磁作动器多目标优化

7.2.1 问题描述

设计一个电磁作动器，在最大化推力的同时最小化能耗。

设计变量：

• 线圈匝数 N ∈ [100, 1000]
• 线圈电流 I ∈ [0.1, 5.0] A
• 气隙长度 g ∈ [0.1, 2.0] mm
• 铁芯半径 r ∈ [5, 50] mm

目标函数：

• f₁：最大化推力 F
• f₂：最小化功耗 P = I²R

约束条件：

• 温升 ΔT ≤ 50°C
• 线圈外径 ≤ 100 mm

7.2.2 Python实现

"""
案例2：电磁作动器多目标优化
使用NSGA-II算法进行多目标优化
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle, FancyBboxPatch

print("="*60)
print("案例2：电磁作动器多目标优化")
print("="*60)

# 物理常数
mu_0 = 4 * np.pi * 1e-7  # 真空磁导率 (H/m)
rho_cu = 1.68e-8  # 铜电阻率 (Ω·m)
k_cu = 400  # 铜导热系数 (W/(m·K))
h_conv = 25  # 对流换热系数 (W/(m²·K))
T_ambient = 25  # 环境温度 (°C)

# 设计变量边界
bounds = {
    'N': (100, 1000),      # 匝数
    'I': (0.1, 5.0),       # 电流 (A)
    'g': (0.1e-3, 2.0e-3), # 气隙 (m)
    'r': (5e-3, 50e-3)     # 铁芯半径 (m)
}

def evaluate_actuator(x):
    """
    评估电磁作动器性能
    x = [N, I, g, r]
    返回: [推力, 功耗, 温升, 外径]
    """
    N, I, g, r = x
    
    # 几何参数
    A = np.pi * r**2  # 铁芯截面积
    coil_thickness = 5e-3  # 线圈厚度 (m)
    coil_length = 30e-3    # 线圈长度 (m)
    wire_diameter = 0.5e-3  # 线径 (m)
    
    # 磁路计算（简化）
    # 磁阻
    R_gap = g / (mu_0 * A)
    R_iron = 0.1 * R_gap  # 铁芯磁阻（假设为气隙的10%）
    R_total = R_gap + R_iron
    
    # 磁动势
    MMF = N * I
    
    # 磁通
    Phi = MMF / R_total
    B = Phi / A  # 磁密
    
    # 电磁推力
    F = B**2 * A / (2 * mu_0)
    
    # 线圈电阻
    wire_length = N * 2 * np.pi * (r + coil_thickness/2)
    wire_area = np.pi * (wire_diameter/2)**2
    R_coil = rho_cu * wire_length / wire_area
    
    # 功耗
    P = I**2 * R_coil
    
    # 温升（简化计算）
    coil_surface = 2 * np.pi * (r + coil_thickness) * coil_length
    delta_T = P / (h_conv * coil_surface)
    
    # 外径
    outer_diameter = 2 * (r + coil_thickness)
    
    return F, P, delta_T, outer_diameter

def check_constraints(x):
    """检查约束条件"""
    _, _, delta_T, outer_dia = evaluate_actuator(x)
    
    g1 = delta_T - 50  # 温升约束
    g2 = outer_dia - 100e-3  # 外径约束
    
    return g1 <= 0 and g2 <= 0

# NSGA-II算法实现
class NSGA2:
    def __init__(self, pop_size=100, n_gen=200, n_var=4, n_obj=2):
        self.pop_size = pop_size
        self.n_gen = n_gen
        self.n_var = n_var
        self.n_obj = n_obj
        
    def initialize_population(self):
        """初始化种群"""
        pop = np.zeros((self.pop_size, self.n_var))
        pop[:, 0] = np.random.uniform(bounds['N'][0], bounds['N'][1], self.pop_size)
        pop[:, 1] = np.random.uniform(bounds['I'][0], bounds['I'][1], self.pop_size)
        pop[:, 2] = np.random.uniform(bounds['g'][0], bounds['g'][1], self.pop_size)
        pop[:, 3] = np.random.uniform(bounds['r'][0], bounds['r'][1], self.pop_size)
        return pop
    
    def evaluate_population(self, pop):
        """评估种群"""
        objectives = np.zeros((len(pop), self.n_obj))
        constraints = np.zeros(len(pop))
        
        for i, x in enumerate(pop):
            F, P, delta_T, outer_dia = evaluate_actuator(x)
            objectives[i, 0] = -F  # 最大化推力
            objectives[i, 1] = P   # 最小化功耗
            
            # 约束 violation
            g1 = max(0, delta_T - 50)
            g2 = max(0, outer_dia - 100e-3)
            constraints[i] = g1 + g2
        
        return objectives, constraints
    
    def non_dominated_sort(self, objectives, constraints):
        """非支配排序"""
        n = len(objectives)
        domination_count = np.zeros(n)
        dominated_solutions = [[] for _ in range(n)]
        fronts = [[]]
        
        for i in range(n):
            for j in range(i+1, n):
                # 考虑约束的支配关系
                if constraints[i] < constraints[j]:
                    dominated_solutions[i].append(j)
                    domination_count[j] += 1
                elif constraints[i] > constraints[j]:
                    dominated_solutions[j].append(i)
                    domination_count[i] += 1
                else:
                    # 约束相同，比较目标
                    obj_i = objectives[i]
                    obj_j = objectives[j]
                    
                    dominates_i = all(obj_i <= obj_j) and any(obj_i < obj_j)
                    dominates_j = all(obj_j <= obj_i) and any(obj_j < obj_i)
                    
                    if dominates_i:
                        dominated_solutions[i].append(j)
                        domination_count[j] += 1
                    elif dominates_j:
                        dominated_solutions[j].append(i)
                        domination_count[i] += 1
            
            if domination_count[i] == 0:
                fronts[0].append(i)
        
        i = 0
        while len(fronts[i]) > 0:
            next_front = []
            for p in fronts[i]:
                for q in dominated_solutions[p]:
                    domination_count[q] -= 1
                    if domination_count[q] == 0:
                        next_front.append(q)
            i += 1
            fronts.append(next_front)
        
        fronts.pop()
        return fronts
    
    def crowding_distance(self, objectives, front):
        """计算拥挤距离"""
        if len(front) <= 2:
            return np.inf * np.ones(len(front))
        
        distances = np.zeros(len(front))
        
        for m in range(self.n_obj):
            sorted_indices = np.argsort(objectives[front, m])
            f_max = objectives[front[sorted_indices[-1]], m]
            f_min = objectives[front[sorted_indices[0]], m]
            
            distances[sorted_indices[0]] = distances[sorted_indices[-1]] = np.inf
            
            for i in range(1, len(front)-1):
                distances[sorted_indices[i]] += (
                    objectives[front[sorted_indices[i+1]], m] - 
                    objectives[front[sorted_indices[i-1]], m]
                ) / (f_max - f_min + 1e-10)
        
        return distances
    
    def tournament_selection(self, pop, objectives, constraints):
        """锦标赛选择"""
        selected = []
        for _ in range(self.pop_size):
            i, j = np.random.randint(0, len(pop), 2)
            
            # 首先比较约束 violation
            if constraints[i] < constraints[j]:
                selected.append(i)
            elif constraints[j] < constraints[i]:
                selected.append(j)
            else:
                # 约束相同，比较支配关系
                obj_i = objectives[i]
                obj_j = objectives[j]
                
                dominates_i = all(obj_i <= obj_j) and any(obj_i < obj_j)
                dominates_j = all(obj_j <= obj_i) and any(obj_j < obj_i)
                
                if dominates_i:
                    selected.append(i)
                elif dominates_j:
                    selected.append(j)
                else:
                    selected.append(i if np.random.random() < 0.5 else j)
        
        return pop[selected]
    
    def crossover(self, parent1, parent2):
        """模拟二进制交叉(SBX)"""
        eta = 20
        child1, child2 = np.copy(parent1), np.copy(parent2)
        
        for i in range(self.n_var):
            if np.random.random() < 0.9:
                if abs(parent1[i] - parent2[i]) > 1e-14:
                    if parent1[i] < parent2[i]:
                        y1, y2 = parent1[i], parent2[i]
                    else:
                        y1, y2 = parent2[i], parent1[i]
                    
                    beta = 1.0 + (2.0 * (y1 - [bounds['N'][0], bounds['I'][0], 
                                           bounds['g'][0], bounds['r'][0]][i]) / (y2 - y1))
                    alpha = 2.0 - beta**(-(eta + 1))
                    
                    rand = np.random.random()
                    if rand <= 1.0 / alpha:
                        beta_q = (rand * alpha)**(1.0 / (eta + 1))
                    else:
                        beta_q = (1.0 / (2.0 - rand * alpha))**(1.0 / (eta + 1))
                    
                    c1 = 0.5 * ((y1 + y2) - beta_q * (y2 - y1))
                    c2 = 0.5 * ((y1 + y2) + beta_q * (y2 - y1))
                    
                    child1[i] = c1
                    child2[i] = c2
        
        return child1, child2
    
    def mutate(self, individual):
        """多项式变异"""
        eta_m = 20
        
        for i in range(self.n_var):
            if np.random.random() < 1.0 / self.n_var:
                y = individual[i]
                yl, yu = [bounds['N'], bounds['I'], bounds['g'], bounds['r']][i]
                
                delta1 = (y - yl) / (yu - yl)
                delta2 = (yu - y) / (yu - yl)
                
                rand = np.random.random()
                mut_pow = 1.0 / (eta_m + 1.0)
                
                if rand <= 0.5:
                    xy = 1.0 - delta1
                    val = 2.0 * rand + (1.0 - 2.0 * rand) * (xy**(eta_m + 1))
                    delta_q = val**mut_pow - 1.0
                else:
                    xy = 1.0 - delta2
                    val = 2.0 * (1.0 - rand) + 2.0 * (rand - 0.5) * (xy**(eta_m + 1))
                    delta_q = 1.0 - val**mut_pow
                
                y = y + delta_q * (yu - yl)
                individual[i] = np.clip(y, yl, yu)
        
        return individual
    
    def optimize(self):
        """主优化循环"""
        print("\n开始NSGA-II优化...")
        
        # 初始化
        pop = self.initialize_population()
        objectives, constraints = self.evaluate_population(pop)
        
        for gen in range(self.n_gen):
            # 选择
            mating_pool = self.tournament_selection(pop, objectives, constraints)
            
            # 交叉和变异
            offspring = []
            for i in range(0, self.pop_size, 2):
                parent1 = mating_pool[i]
                parent2 = mating_pool[(i+1) % self.pop_size]
                
                child1, child2 = self.crossover(parent1, parent2)
                child1 = self.mutate(child1)
                child2 = self.mutate(child2)
                
                offspring.extend([child1, child2])