电磁场仿真-主题077-拓扑优化
·
第077篇：拓扑优化
摘要
拓扑优化是结构优化设计领域最具革命性的技术之一，它能够在给定的设计空间内自动寻找最优的材料分布，无需预先设定结构形式。在电磁场仿真领域，拓扑优化技术已被广泛应用于天线设计、波导结构优化、超材料设计、微波器件设计等众多方向。本教程系统介绍电磁场拓扑优化的核心理论与方法，包括SIMP材料插值模型、密度过滤与投影技术、伴随法灵敏度分析、优化算法（OC方法与MMA方法），并通过波导优化、天线结构优化和多目标优化三个典型案例展示拓扑优化在电磁工程中的实际应用。通过本教程的学习，读者将掌握拓扑优化的基本原理和实现方法，能够独立开展电磁结构的拓扑优化设计工作。
关键词
拓扑优化，SIMP方法，伴随法灵敏度分析，密度过滤，波导优化，天线设计，多目标优化，电磁结构优化
1. 拓扑优化概述
1.1 什么是拓扑优化
拓扑优化（Topology Optimization）是一种数学方法，用于在给定的设计空间、载荷条件和约束条件下，寻找最优的材料分布方案。与传统的尺寸优化和形状优化不同，拓扑优化能够在设计过程中改变结构的拓扑形式——即材料的有无和连通方式，从而获得创新的结构设计。
在电磁工程领域，拓扑优化的目标是找到最优的介电常数或磁导率分布，使得电磁器件的性能达到最优。例如：
- 天线设计：优化天线辐射体的形状，实现指定的辐射方向图
- 波导设计：设计波导结构，实现高效的电磁波传输
- 超材料设计：构建具有特殊电磁响应的人工结构
- 滤波器设计：优化滤波器拓扑，获得理想的频率响应
1.2 拓扑优化的数学表述
拓扑优化问题可以表述为一个约束优化问题：
minρJ(ρ,u(ρ))\min_{\rho} J(\rho, u(\rho))ρminJ(ρ,u(ρ))
s.t.gj(ρ)≤0,j=1,2,...,m\text{s.t.} \quad g_j(\rho) \leq 0, \quad j = 1, 2, ..., ms.t.gj(ρ)≤0,j=1,2,...,m
ρmin≤ρ≤ρmax\rho_{\min} \leq \rho \leq \rho_{\max}ρmin≤ρ≤ρmax
其中：
- ρ\rhoρ 是设计变量（密度场）
- JJJ 是目标函数（如传输系数、反射系数等）
- u(ρ)u(\rho)u(ρ) 是状态变量（电磁场分布）
- gjg_jgj 是约束函数（如体积约束、带宽约束等）
1.3 拓扑优化的发展历程
拓扑优化的发展经历了以下重要阶段：
1988年：Bendsøe和Kikuchi提出均匀化方法（Homogenization Method），奠定了拓扑优化的理论基础。
1989年：Bendsøe提出SIMP方法（Solid Isotropic Material with Penalization），成为最广泛使用的拓扑优化方法。
1990年代：Sigmund等发展了密度过滤技术，解决了棋盘格现象和网格依赖性等问题。
2000年代：Guest等提出密度投影方法，实现了清晰的0/1解。
2010年代至今：拓扑优化在电磁领域的应用快速发展，出现了水平集方法、进化结构优化（ESO）、双向进化结构优化（BESO）等多种方法。
1.4 电磁拓扑优化的特点
电磁拓扑优化相比结构力学拓扑优化具有以下特点：
- 多物理场耦合：电磁问题通常涉及电场、磁场、热场等多物理场的耦合
- 频域特性：电磁响应具有强烈的频率依赖性，需要在宽频带内优化
- 共振现象：电磁结构常出现共振，导致目标函数的高度非线性
- 制造约束：电磁结构的制造需要考虑材料特性、加工工艺等约束
2. SIMP材料插值模型
2.1 SIMP方法的基本原理
SIMP（Solid Isotropic Material with Penalization）方法是拓扑优化中最常用的材料插值方法。其核心思想是通过引入惩罚因子，使中间密度（0 < ρ < 1）的材料属性变得不经济，从而推动优化结果向0/1的离散解收敛。
SIMP插值公式为：
ε(ρ)=εmin+(εmax−εmin)ρp\varepsilon(\rho) = \varepsilon_{\min} + (\varepsilon_{\max} - \varepsilon_{\min}) \rho^pε(ρ)=εmin+(εmax−εmin)ρp
其中：
- ε(ρ)\varepsilon(\rho)ε(ρ) 是等效介电常数
- ρ∈[0,1]\rho \in [0, 1]ρ∈[0,1] 是单元密度
- εmin\varepsilon_{\min}εmin 是空材料的介电常数（通常接近1）
- εmax\varepsilon_{\max}εmax 是实体材料的介电常数
- ppp 是惩罚因子（通常取3）
2.2 惩罚因子的作用
惩罚因子 ppp 的选择对优化结果有重要影响：
- p = 1：线性插值，优化结果会出现大量中间密度（"灰度"现象）
- p = 3：标准惩罚，能有效抑制中间密度
- p > 3：强惩罚，可能导致优化收敛困难
惩罚因子的物理意义是：当 p>1p > 1p>1 时，中间密度材料的"性价比"低于纯材料。例如，密度为0.5的材料，其有效介电常数为 0.53×εmax=0.125×εmax0.5^3 \times \varepsilon_{\max} = 0.125 \times \varepsilon_{\max}0.53×εmax=0.125×εmax，远低于线性插值的 0.5×εmax0.5 \times \varepsilon_{\max}0.5×εmax。
2.3 SIMP方法的灵敏度分析
SIMP插值对密度的导数为：
∂ε∂ρ=p(εmax−εmin)ρp−1\frac{\partial \varepsilon}{\partial \rho} = p (\varepsilon_{\max} - \varepsilon_{\min}) \rho^{p-1}∂ρ∂ε=p(εmax−εmin)ρp−1
这个导数在灵敏度分析中起着关键作用，它描述了目标函数对设计变量的敏感程度。
2.4 SIMP方法的优缺点
优点：
- 实现简单，计算效率高
- 与有限元方法天然兼容
- 经过大量工程验证，可靠性高
缺点：
- 存在网格依赖性（需要配合过滤技术）
- 可能产生棋盘格现象
- 对复杂的多材料问题处理能力有限
3. 密度过滤与投影技术
3.1 棋盘格现象与网格依赖性
在拓扑优化中，如果不采取特殊措施，优化结果会出现以下数值问题：
棋盘格现象（Checkerboard Pattern）：优化结果呈现棋盘状的黑白交替分布，这种结构在物理上是不合理的，也无法制造。
网格依赖性（Mesh Dependency）：优化结果严重依赖于网格划分，网格加密后可能得到完全不同的拓扑。
这些问题的根源在于：
- 数值计算中的伪模式
- 缺乏对最小特征尺寸的控制
- 优化问题的数学病态性
3.2 密度过滤方法
密度过滤是解决上述问题的有效方法。其基本思想是：每个单元的物理密度不是独立的设计变量，而是周围单元设计密度的加权平均。
过滤公式为：
ρ~e=∑i∈Neweiρi∑i∈Newei\tilde{\rho}_e = \frac{\sum_{i \in N_e} w_{ei} \rho_i}{\sum_{i \in N_e} w_{ei}}ρ~e=∑i∈Newei∑i∈Neweiρi
其中：
- ρ~e\tilde{\rho}_eρ~e 是单元e的过滤后密度
- NeN_eNe 是单元e的邻域
- weiw_{ei}wei 是权重函数，通常采用线性衰减：wei=max(0,rmin−dei)w_{ei} = \max(0, r_{\min} - d_{ei})wei=max(0,rmin−dei)
- rminr_{\min}rmin 是过滤半径
- deid_{ei}dei 是单元e和i之间的距离
过滤半径 rminr_{\min}rmin 的选择很重要：
- 太小：无法有效消除棋盘格
- 太大：过度平滑，细节丢失
- 经验值：3-5倍单元尺寸
3.3 密度投影方法
虽然密度过滤消除了棋盘格，但优化结果仍可能包含大量的中间密度（灰度）。密度投影技术通过非线性映射，将过滤后的密度投影到接近0或1的值。
常用的Heaviside投影函数为：
ρˉ=tanh(βη)+tanh(β(ρ~−η))tanh(βη)+tanh(β(1−η))\bar{\rho} = \frac{\tanh(\beta \eta) + \tanh(\beta (\tilde{\rho} - \eta))}{\tanh(\beta \eta) + \tanh(\beta (1 - \eta))}ρˉ=tanh(βη)+tanh(β(1−η))tanh(βη)+tanh(β(ρ~−η))
其中：
- ρˉ\bar{\rho}ρˉ 是投影后的密度
- ρ~\tilde{\rho}ρ~ 是过滤后的密度
- β\betaβ 是投影锐度参数（β越大，投影越锐利）
- η\etaη 是投影阈值（通常取0.5）
投影锐度参数 β\betaβ 通常采用Continuation策略：在优化初期使用较小的β（如1），随着优化进行逐渐增大（如到32），以获得清晰的0/1解。
3.4 灵敏度过滤
当使用密度过滤时，目标函数对设计变量的灵敏度也需要相应调整。根据链式法则：
∂J∂ρi=∑e∂J∂ρ~e∂ρ~e∂ρi\frac{\partial J}{\partial \rho_i} = \sum_{e} \frac{\partial J}{\partial \tilde{\rho}_e} \frac{\partial \tilde{\rho}_e}{\partial \rho_i}∂ρi∂J=e∑∂ρ~e∂J∂ρi∂ρ~e
其中：
∂ρ~e∂ρi=wei∑j∈Newej\frac{\partial \tilde{\rho}_e}{\partial \rho_i} = \frac{w_{ei}}{\sum_{j \in N_e} w_{ej}}∂ρi∂ρ~e=∑j∈Newejwei
这种灵敏度过滤确保了优化算法的正确收敛。
4. 伴随法灵敏度分析
4.1 灵敏度分析的重要性
在拓扑优化中，需要计算目标函数和约束函数对设计变量的导数（灵敏度）。对于大规模问题（设计变量数可达百万级），直接计算这些导数的计算成本极高。伴随法（Adjoint Method）提供了一种高效的灵敏度计算方法。
4.2 伴随法的基本原理
考虑一个一般形式的优化问题：
J(ρ)=J(ρ,u(ρ))J(\rho) = J(\rho, u(\rho))J(ρ)=J(ρ,u(ρ))
其中 uuu 是状态变量（电磁场），满足状态方程：
R(ρ,u)=0R(\rho, u) = 0R(ρ,u)=0
目标函数的全导数为：
dJdρ=∂J∂ρ+∂J∂u∂u∂ρ\frac{dJ}{d\rho} = \frac{\partial J}{\partial \rho} + \frac{\partial J}{\partial u} \frac{\partial u}{\partial \rho}dρdJ=∂ρ∂J+∂u∂J∂ρ∂u
直接计算 ∂u∂ρ\frac{\partial u}{\partial \rho}∂ρ∂u 需要求解多个状态方程（每个设计变量一个）。伴随法通过引入伴随变量 λ\lambdaλ，避免了这一昂贵的计算。
伴随方程为：
(∂R∂u)Tλ=−(∂J∂u)T\left(\frac{\partial R}{\partial u}\right)^T \lambda = -\left(\frac{\partial J}{\partial u}\right)^T(∂u∂R)Tλ=−(∂u∂J)T
灵敏度可以高效计算为：
dJdρ=∂J∂ρ+λT∂R∂ρ\frac{dJ}{d\rho} = \frac{\partial J}{\partial \rho} + \lambda^T \frac{\partial R}{\partial \rho}dρdJ=∂ρ∂J+λT∂ρ∂R
4.3 电磁场问题的伴随方程
对于电磁场问题，状态方程通常是Maxwell方程的离散形式。以频域电磁场为例：
∇×(1μ∇×E)−k02εE=−jk0Z0J\nabla \times \left(\frac{1}{\mu} \nabla \times \mathbf{E}\right) - k_0^2 \varepsilon \mathbf{E} = -j k_0 Z_0 \mathbf{J}∇×(μ1∇×E)−k02εE=−jk0Z0J
伴随方程为：
∇×(1μ∇×λ)−k02ελ=−∂J∂E\nabla \times \left(\frac{1}{\mu} \nabla \times \boldsymbol{\lambda}\right) - k_0^2 \varepsilon \boldsymbol{\lambda} = -\frac{\partial J}{\partial \mathbf{E}}∇×(μ1∇×λ)−k02ελ=−∂E∂J
注意到伴随方程与原始方程具有相同的形式，只是源项不同。这意味着可以使用相同的求解器来求解原始场和伴随场。
4.4 伴随法的计算效率
伴随法的计算效率优势在于：
- 直接法：计算n个设计变量的灵敏度需要求解n次状态方程
- 伴随法：无论设计变量数多少，只需求解1次原始方程 + 1次伴随方程
对于大规模拓扑优化问题（n ~ 10^6），伴随法可以将计算时间从数天缩短到数小时。
5. 优化算法
5.1 最优性准则方法（OC）
最优性准则（Optimality Criteria, OC）方法是拓扑优化中最常用的算法之一。它基于KKT最优性条件，通过启发式更新规则迭代改进设计。
OC更新公式为：
ρenew={max(ρmin,ρe−m)if ρeBeη≤max(ρmin,ρe−m)min(1,ρe+m)if ρeBeη≥min(1,ρe+m)ρeBeηotherwise\rho_e^{new} = \begin{cases} \max(\rho_{\min}, \rho_e - m) & \text{if } \rho_e B_e^\eta \leq \max(\rho_{\min}, \rho_e - m) \\ \min(1, \rho_e + m) & \text{if } \rho_e B_e^\eta \geq \min(1, \rho_e + m) \\ \rho_e B_e^\eta & \text{otherwise} \end{cases}ρenew=⎩ ⎨ ⎧max(ρmin,ρe−m)min(1,ρe+m)ρeBeηif ρeBeη≤max(ρmin,ρe−m)if ρeBeη≥min(1,ρe+m)otherwise
其中：
- Be=−∂J/∂ρeλ∂V/∂ρeB_e = -\frac{\partial J/\partial \rho_e}{\lambda \partial V/\partial \rho_e}Be=−λ∂V/∂ρe∂J/∂ρe
- η\etaη 是阻尼因子（通常取0.5）
- mmm 是移动限制（通常取0.2）
- λ\lambdaλ 是拉格朗日乘子，通过二分法确定以满足体积约束
OC方法的优点：
- 实现简单，计算效率高
- 收敛稳定，对初始设计不敏感
- 适用于单约束问题
5.2 移动渐近线方法（MMA）
移动渐近线方法（Method of Moving Asymptotes, MMA）是一种更通用的优化算法，适用于多约束问题。
MMA的核心思想是：在当前设计点处构建一个凸的近似子问题，通过移动渐近线来控制近似质量。
MMA子问题形式为：
minx∑i=1n(pi(k)Ui(k)−xi+qi(k)xi−Li(k))+r0(k)\min_{x} \sum_{i=1}^{n} \left(\frac{p_i^{(k)}}{U_i^{(k)} - x_i} + \frac{q_i^{(k)}}{x_i - L_i^{(k)}}\right) + r_0^{(k)}xmini=1∑n(Ui(k)−xipi(k)+xi−Li(k)qi(k))+r0(k)
s.t.∑i=1n(pij(k)Ui(k)−xi+qij(k)xi−Li(k))+rj(k)≤0,j=1,...,m\text{s.t.} \quad \sum_{i=1}^{n} \left(\frac{p_{ij}^{(k)}}{U_i^{(k)} - x_i} + \frac{q_{ij}^{(k)}}{x_i - L_i^{(k)}}\right) + r_j^{(k)} \leq 0, \quad j = 1, ..., ms.t.i=1∑n(Ui(k)−xipij(k)+xi−Li(k)qij(k))+rj(k)≤0,j=1,...,m
αi(k)≤xi≤βi(k),i=1,...,n\alpha_i^{(k)} \leq x_i \leq \beta_i^{(k)}, \quad i = 1, ..., nαi(k)≤xi≤βi(k),i=1,...,n
其中 Li(k)L_i^{(k)}Li(k) 和 Ui(k)U_i^{(k)}Ui(k) 是移动渐近线，根据优化历史动态调整。
MMA的优点：
- 适用于多约束问题
- 收敛速度快
- 数值稳定性好
5.3 其他优化算法
序列线性规划（SLP）：在每一步将问题线性化，求解线性规划子问题。
序列二次规划（SQP）：在每一步构建二次近似，求解二次规划子问题。
进化算法：如遗传算法、粒子群优化等，适用于高度非凸问题，但计算成本高。
6. 波导拓扑优化
6.1 波导优化的目标与约束
波导拓扑优化的典型目标是最大化特定频率下的功率传输系数，同时满足体积约束：
maxρT(ρ)=PoutPin\max_{\rho} T(\rho) = \frac{P_{out}}{P_{in}}ρmaxT(ρ)=PinPout
s.t.∫Ωρ dΩ∣Ω∣≤Vmax\text{s.t.} \quad \frac{\int_{\Omega} \rho \, d\Omega}{|\Omega|} \leq V_{\max}s.t.∣Ω∣∫ΩρdΩ≤Vmax
其中 TTT 是传输系数，VmaxV_{\max}Vmax 是最大允许体积分数。
6.2 波导优化的物理模型
波导中的电磁场传播可以用Helmholtz方程描述：
∇2E+k02ε(ρ)E=0\nabla^2 E + k_0^2 \varepsilon(\rho) E = 0∇2E+k02ε(ρ)E=0
其中 k0=ω/ck_0 = \omega/ck0=ω/c 是自由空间波数，ε(ρ)\varepsilon(\rho)ε(ρ) 是SIMP插值后的介电常数。
传输系数可以通过计算输出端口的功率流获得：
T=Re{∫SoutE×H∗⋅dS}Re{∫SinE×H∗⋅dS}T = \frac{\text{Re}\{\int_{S_{out}} \mathbf{E} \times \mathbf{H}^* \cdot d\mathbf{S}\}}{\text{Re}\{\int_{S_{in}} \mathbf{E} \times \mathbf{H}^* \cdot d\mathbf{S}\}}T=Re{∫SinE×H∗⋅dS}Re{∫SoutE×H∗⋅dS}
6.3 波导优化的灵敏度
使用伴随法，传输系数对设计变量的灵敏度为：
dTdρe=−k02∂ε∂ρe∫ΩeE⋅λ dΩ\frac{dT}{d\rho_e} = -k_0^2 \frac{\partial \varepsilon}{\partial \rho_e} \int_{\Omega_e} E \cdot \lambda \, d\OmegadρedT=−k02∂ρe∂ε∫ΩeE⋅λdΩ
其中 λ\lambdaλ 是伴随场，满足伴随方程：
∇2λ+k02ε(ρ)λ=−∂T∂E\nabla^2 \lambda + k_0^2 \varepsilon(\rho) \lambda = -\frac{\partial T}{\partial E}∇2λ+k02ε(ρ)λ=−∂E∂T
6.4 波导优化案例
考虑一个直波导的拓扑优化问题：
- 设计域：60×40网格
- 目标频率：3 GHz
- 目标体积分数：40%
- 惩罚因子：p = 3
- 过滤半径：3个单元
优化过程通常需要30-50次迭代收敛。优化结果会形成一个清晰的波导结构，中间为实体材料（ε = 12），周围为空气（ε = 1）。
7. 天线结构拓扑优化
7.1 天线优化的特殊考虑
天线拓扑优化相比波导优化有以下特殊考虑：
- 宽带特性：天线通常需要在宽频带内工作
- 辐射方向图：需要优化辐射方向图，而不仅仅是传输系数
- 阻抗匹配：需要同时优化阻抗匹配特性
- 制造约束：天线结构通常需要满足特定的制造约束
7.2 天线优化的目标函数
天线优化的目标函数可以包括：
带宽最大化：
J=−minf∈[fmin,fmax]∣S11(f)∣J = -\min_{f \in [f_{\min}, f_{\max}]} |S_{11}(f)|J=−f∈[fmin,fmax]min∣S11(f)∣
方向图优化：
J=∫θ∣F(θ)−Ftarget(θ)∣2dθJ = \int_{\theta} |F(\theta) - F_{target}(\theta)|^2 d\thetaJ=∫θ∣F(θ)−Ftarget(θ)∣2dθ
多目标优化：
J=w1Jbandwidth+w2Jpattern+w3JimpedanceJ = w_1 J_{bandwidth} + w_2 J_{pattern} + w_3 J_{impedance}J=w1Jbandwidth+w2Jpattern+w3Jimpedance
7.3 天线优化的设计域
天线拓扑优化通常采用以下设计策略：
- 辐射体优化：优化天线的辐射结构形状
- 地板优化：优化接地板的形状以改善辐射特性
- 馈电网络优化：优化馈电网络的结构
设计域通常从中心向外扩展，初始设计为简单的圆形或方形区域。
7.4 天线优化案例
考虑一个宽带天线的拓扑优化：
- 设计域：50×50网格
- 频率范围：2-6 GHz
- 目标体积分数：30%
- 优化目标：最大化带宽指标
优化结果会形成一个具有复杂形状的天线结构，这种结构通常难以通过传统设计方法获得。
8. 多目标拓扑优化
8.1 多目标优化的概念
实际工程问题通常涉及多个相互冲突的目标，例如：
- 最大化性能 vs 最小化材料用量
- 最大化带宽 vs 最小化尺寸
- 最大化增益 vs 简化制造
多目标优化旨在找到一组Pareto最优解，在这些解中，无法在不牺牲其他目标的情况下改进任何一个目标。
8.2 多目标优化方法
加权和方法：
J=∑i=1nwiJi,∑i=1nwi=1J = \sum_{i=1}^{n} w_i J_i, \quad \sum_{i=1}^{n} w_i = 1J=i=1∑nwiJi,i=1∑nwi=1
通过改变权重组合，可以获得Pareto前沿上的不同解。
ε-约束方法：将除一个目标外的所有目标转化为约束。
进化算法：如NSGA-II，可以直接搜索Pareto前沿。
8.3 Pareto前沿分析
Pareto前沿描述了不同目标之间的权衡关系。通过分析Pareto前沿，工程师可以：
- 理解目标之间的冲突程度
- 选择最符合工程需求的折衷方案
- 评估设计改进的潜力
8.4 多目标优化案例
考虑一个同时优化连通性和平滑度的双目标问题：
- 目标1：最大化结构连通性
- 目标2：最大化介电常数分布的平滑度
通过改变权重组合（如w1=0.9, w2=0.1到w1=0.1, w2=0.9），可以获得一系列不同的优化设计，形成Pareto前沿。
9. 拓扑优化的数值实现
9.1 实现流程
拓扑优化的数值实现通常遵循以下流程：
1. 初始化
- 设置网格和设计变量
- 初始化密度场
- 设置优化参数
2. 迭代优化
for iter = 1 to max_iter:
a. 过滤和投影密度场
b. 计算材料属性
c. 求解状态方程（电磁场）
d. 计算目标函数和约束
e. 求解伴随方程（如需要）
f. 计算灵敏度
g. 过滤灵敏度
h. 更新设计变量
i. 检查收敛性
3. 输出结果
- 保存最终设计
- 绘制收敛历史
- 后处理和可视化
9.2 收敛准则
常用的收敛准则包括：
设计变量变化：
∥ρ(k+1)−ρ(k)∥∥ρ(k)∥<ε1\frac{\|\rho^{(k+1)} - \rho^{(k)}\|}{\|\rho^{(k)}\|} < \varepsilon_1∥ρ(k)∥∥ρ(k+1)−ρ(k)∥<ε1
目标函数变化：
∣J(k+1)−J(k)∣∣J(k)∣<ε2\frac{|J^{(k+1)} - J^{(k)}|}{|J^{(k)}|} < \varepsilon_2∣J(k)∣∣J(k+1)−J(k)∣<ε2
最大迭代次数：达到预设的最大迭代次数
9.3 数值稳定性
拓扑优化中常见的数值问题及解决方法：
NaN值：
- 原因：数值不稳定、除零、不合理的材料参数
- 解决：添加正则化项、检查边界条件、限制设计变量范围
振荡：
- 原因：移动限制过大、目标函数高度非线性
- 解决：减小移动限制、增加阻尼、使用过滤技术
不收敛：
- 原因：问题病态、约束冲突
- 解决：检查问题设置、放松约束、使用鲁棒优化算法
10. 拓扑优化的工程应用
10.1 微波器件设计
拓扑优化在微波器件设计中的应用包括：
- 滤波器设计：优化滤波器拓扑，实现陡峭的截止特性
- 耦合器设计：设计定向耦合器，实现精确的功率分配
- 功分器设计：优化功分器结构，实现宽频带工作
10.2 天线阵列设计
拓扑优化在天线阵列设计中的应用：
- 阵元布局优化：优化阵元位置，降低副瓣电平
- 馈电网络优化：设计紧凑的馈电网络
- 超表面设计：设计具有特殊功能的超表面结构
10.3 超材料设计
拓扑优化在超材料设计中的应用：
- 负折射率材料：设计具有负折射率的结构
- 隐身斗篷：优化隐身结构，实现电磁波绕射
- 完美吸收体：设计宽带完美吸收结构
10.4 制造约束的处理
实际制造约束需要在拓扑优化中考虑：
最小特征尺寸：通过过滤半径控制
连通性约束：确保结构连通，避免孤岛
对称性约束：施加对称性以减少制造复杂度
各向异性制造：考虑制造工艺的各向异性特性
11. Python仿真代码实现
以下是完整的拓扑优化仿真代码，实现了本章介绍的所有核心算法：
# -*- coding: utf-8 -*-
"""
主题077：拓扑优化仿真
Topology Optimization for Electromagnetic Applications
本程序实现电磁场仿真中的拓扑优化算法，包括：
1. SIMP材料插值模型
2. 密度过滤与投影
3. 伴随法灵敏度分析
4. 优化算法（OC、MMA）
5. 电磁波导优化
6. 天线结构优化
7. 多目标拓扑优化
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.sparse import diags, csr_matrix
from scipy.sparse.linalg import spsolve
from scipy.ndimage import gaussian_filter
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("=" * 60)
print("拓扑优化仿真程序")
print("=" * 60)
# =============================================================================
# 第一部分：SIMP材料插值模型
# =============================================================================
print("\n[1] SIMP材料插值模型...")
class SIMPMaterial:
"""SIMP (Solid Isotropic Material with Penalization) 材料模型"""
def __init__(self, epsilon_min=1.0, epsilon_max=10.0, p=3.0):
"""
参数:
epsilon_min: 空材料介电常数 (通常接近1)
epsilon_max: 实体材料介电常数
p: 惩罚因子
"""
self.epsilon_min = epsilon_min
self.epsilon_max = epsilon_max
self.p = p
def interpolate(self, rho):
"""
SIMP插值
参数:
rho: 密度场 (0到1之间)
返回:
epsilon: 等效介电常数
"""
return self.epsilon_min + (self.epsilon_max - self.epsilon_min) * rho**self.p
def derivative(self, rho):
"""
SIMP插值对密度的导数
参数:
rho: 密度场
返回:
depsilon_drho: 导数
"""
return self.p * (self.epsilon_max - self.epsilon_min) * rho**(self.p - 1)
# 演示SIMP插值
simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=10.0, p=3.0)
rho_demo = np.linspace(0, 1, 100)
epsilon_demo = simp.interpolate(rho_demo)
def plot_simp_model():
"""绘制SIMP模型"""
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# SIMP插值曲线
for p in [1, 2, 3, 4, 5]:
simp_temp = SIMPMaterial(epsilon_min=1.0, epsilon_max=10.0, p=p)
epsilon_temp = simp_temp.interpolate(rho_demo)
axes[0].plot(rho_demo, epsilon_temp, linewidth=2, label=f'p={p}')
axes[0].axhline(y=1.0, color='k', linestyle='--', alpha=0.3)
axes[0].axhline(y=10.0, color='k', linestyle='--', alpha=0.3)
axes[0].set_xlabel('Density (rho)', fontsize=11)
axes[0].set_ylabel('Relative Permittivity', fontsize=11)
axes[0].set_title('SIMP Interpolation Model', fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 惩罚效果示意
rho_binary = np.where(rho_demo > 0.5, 1.0, 0.0)
axes[1].plot(rho_demo, rho_demo, 'b--', linewidth=2, label='Linear (p=1)')
axes[1].plot(rho_demo, rho_demo**3, 'r-', linewidth=2, label='SIMP (p=3)')
axes[1].plot(rho_demo, rho_binary, 'g:', linewidth=2, label='Binary')
axes[1].set_xlabel('Physical Density', fontsize=11)
axes[1].set_ylabel('Effective Density', fontsize=11)
axes[1].set_title('Penalization Effect', fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/simp_model.png', dpi=150, bbox_inches='tight')
plt.close()
plot_simp_model()
print("✓ SIMP模型图已保存")
# =============================================================================
# 第二部分：密度过滤与投影
# =============================================================================
print("\n[2] 密度过滤与投影...")
class DensityFilter:
"""密度过滤器 - 消除棋盘格现象"""
def __init__(self, nx, ny, rmin):
"""
参数:
nx, ny: 网格数
rmin: 过滤半径（网格单元数）
"""
self.nx = nx
self.ny = ny
self.rmin = rmin
# 构建过滤权重矩阵
self.H, self.Hs = self._build_filter_matrix()
def _build_filter_matrix(self):
"""构建过滤矩阵"""
n = self.nx * self.ny
H = sparse.lil_matrix((n, n))
for i in range(self.nx):
for j in range(self.ny):
e1 = i * self.ny + j
for k in range(max(0, i - int(self.rmin)), min(self.nx, i + int(self.rmin) + 1)):
for l in range(max(0, j - int(self.rmin)), min(self.ny, j + int(self.rmin) + 1)):
e2 = k * self.ny + l
dist = np.sqrt((i - k)**2 + (j - l)**2)
H[e1, e2] = max(0, self.rmin - dist)
H = H.tocsr()
Hs = np.array(H.sum(axis=1)).flatten()
return H, Hs
def filter_density(self, rho):
"""过滤密度场"""
rho_flat = rho.flatten()
rho_filtered = (self.H @ rho_flat) / self.Hs
return np.array(rho_filtered).flatten().reshape(self.nx, self.ny)
def filter_sensitivity(self, d_obj_drho_filtered, rho):
"""过滤灵敏度"""
n = self.nx * self.ny
d_obj_drho = np.zeros(n)
d_obj_drho_filtered_flat = d_obj_drho_filtered.flatten()
for i in range(n):
for j in range(n):
if self.H[i, j] > 0:
d_obj_drho[j] += float(self.H[i, j]) * d_obj_drho_filtered_flat[i] / self.Hs[i]
return d_obj_drho.reshape(self.nx, self.ny)
class DensityProjection:
"""密度投影 - 获得清晰的0/1解"""
def __init__(self, beta=1.0, eta=0.5):
"""
参数:
beta: 投影锐度参数
eta: 投影阈值
"""
self.beta = beta
self.eta = eta
def project(self, rho_tilde):
"""
Heaviside投影
参数:
rho_tilde: 过滤后的密度
返回:
rho: 投影后的密度
"""
return (np.tanh(self.beta * self.eta) + np.tanh(self.beta * (rho_tilde - self.eta))) / \
(np.tanh(self.beta * self.eta) + np.tanh(self.beta * (1 - self.eta)))
def derivative(self, rho_tilde):
"""投影函数的导数"""
numerator = self.beta * (1 - np.tanh(self.beta * (rho_tilde - self.eta))**2)
denominator = np.tanh(self.beta * self.eta) + np.tanh(self.beta * (1 - self.eta))
return numerator / denominator
# 演示过滤与投影
def demo_filter_projection():
"""演示密度过滤和投影效果"""
nx, ny = 60, 40
# 创建带有棋盘格的初始密度场
np.random.seed(42)
rho_raw = np.random.rand(nx, ny)
# 添加棋盘格模式
for i in range(nx):
for j in range(ny):
if (i + j) % 2 == 0:
rho_raw[i, j] = 0.9
else:
rho_raw[i, j] = 0.1
# 应用过滤
df = DensityFilter(nx, ny, rmin=3.0)
rho_filtered = df.filter_density(rho_raw)
# 应用投影
proj = DensityProjection(beta=8.0, eta=0.5)
rho_projected = proj.project(rho_filtered)
# 绘制
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
im0 = axes[0, 0].imshow(rho_raw.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 0].set_title('Raw Density (with checkerboard)', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
plt.colorbar(im0, ax=axes[0, 0])
im1 = axes[0, 1].imshow(rho_filtered.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 1].set_title('Filtered Density', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('y')
plt.colorbar(im1, ax=axes[0, 1])
im2 = axes[1, 0].imshow(rho_projected.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[1, 0].set_title('Projected Density (beta=8)', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('x')
axes[1, 0].set_ylabel('y')
plt.colorbar(im2, ax=axes[1, 0])
# 投影函数曲线
rho_test = np.linspace(0, 1, 100)
for beta in [1, 4, 8, 16]:
proj_temp = DensityProjection(beta=beta, eta=0.5)
rho_proj_test = proj_temp.project(rho_test)
axes[1, 1].plot(rho_test, rho_proj_test, linewidth=2, label=f'beta={beta}')
axes[1, 1].plot([0, 1], [0, 1], 'k--', alpha=0.3)
axes[1, 1].set_xlabel('Filtered Density', fontsize=11)
axes[1, 1].set_ylabel('Projected Density', fontsize=11)
axes[1, 1].set_title('Projection Function', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/density_filter_projection.png', dpi=150, bbox_inches='tight')
plt.close()
demo_filter_projection()
print("✓ 密度过滤与投影图已保存")
# =============================================================================
# 第三部分：伴随法灵敏度分析
# =============================================================================
print("\n[3] 伴随法灵敏度分析...")
class AdjointSensitivity:
"""伴随法灵敏度分析"""
def __init__(self, nx, ny, dx, dy, omega, epsilon):
"""
参数:
nx, ny: 网格数
dx, dy: 网格间距
omega: 角频率
epsilon: 介电常数分布
"""
self.nx = nx
self.ny = ny
self.dx = dx
self.dy = dy
self.omega = omega
self.epsilon = epsilon
self.k0 = omega / 3e8 # 自由空间波数
def build_fem_matrix(self, epsilon_field):
"""构建FEM系统矩阵"""
n = self.nx * self.ny
# 构建离散的Helmholtz算子
# (nabla^2 + k^2) E = 0
row = []
col = []
data = []
for i in range(self.nx):
for j in range(self.ny):
idx = i * self.ny + j
# 对角元素
k2 = (self.k0**2) * epsilon_field[i, j]
coeff = 2/self.dx**2 + 2/self.dy**2 - k2
row.append(idx)
col.append(idx)
data.append(coeff)
# x方向邻居
if i > 0:
row.append(idx)
col.append((i-1) * self.ny + j)
data.append(-1/self.dx**2)
if i < self.nx - 1:
row.append(idx)
col.append((i+1) * self.ny + j)
data.append(-1/self.dx**2)
# y方向邻居
if j > 0:
row.append(idx)
col.append(i * self.ny + (j-1))
data.append(-1/self.dy**2)
if j < self.ny - 1:
row.append(idx)
col.append(i * self.ny + (j+1))
data.append(-1/self.dy**2)
A = csr_matrix((data, (row, col)), shape=(n, n))
return A
def solve_field(self, epsilon_field, source):
"""求解电磁场"""
A = self.build_fem_matrix(epsilon_field)
n = self.nx * self.ny
# 添加简单的边界条件（Dirichlet）
# 简化处理：直接求解
E = spsolve(A, source.flatten())
return E.reshape(self.nx, self.ny)
def compute_sensitivity(self, epsilon_field, E, objective_type='energy'):
"""
计算灵敏度
参数:
epsilon_field: 介电常数场
E: 电场分布
objective_type: 目标函数类型
返回:
sensitivity: 灵敏度场
"""
if objective_type == 'energy':
# 目标函数：电场能量
# J = integral(|E|^2)
# dJ/depsilon = -omega^2 * |E|^2
sensitivity = -self.omega**2 * np.abs(E)**2
elif objective_type == 'transmission':
# 目标函数：传输系数
# 需要伴随场
n = self.nx * self.ny
# 构建伴随源（在输出端）
adj_source = np.zeros((self.nx, self.ny))
adj_source[-5:, self.ny//2-2:self.ny//2+2] = 1.0
# 求解伴随场
A = self.build_fem_matrix(epsilon_field)
lambda_field = spsolve(A.T, adj_source.flatten())
lambda_field = lambda_field.reshape(self.nx, self.ny)
# 灵敏度：dJ/depsilon = -omega^2 * E * lambda
sensitivity = -self.omega**2 * np.real(E * lambda_field)
else:
sensitivity = np.zeros_like(epsilon_field)
return sensitivity
# 演示伴随法灵敏度
def demo_adjoint_sensitivity():
"""演示伴随法灵敏度计算"""
nx, ny = 40, 30
dx = dy = 0.01
omega = 2 * np.pi * 3e9 # 3 GHz
# 初始密度场
rho = np.ones((nx, ny)) * 0.5
simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=10.0, p=3.0)
epsilon_field = simp.interpolate(rho)
# 创建伴随灵敏度分析器
adj = AdjointSensitivity(nx, ny, dx, dy, omega, epsilon_field)
# 创建源
source = np.zeros((nx, ny))
source[5, ny//2-2:ny//2+2] = 1.0 # 输入端口
# 求解场
E = adj.solve_field(epsilon_field, source)
# 计算灵敏度
sens_energy = adj.compute_sensitivity(epsilon_field, E, 'energy')
sens_trans = adj.compute_sensitivity(epsilon_field, E, 'transmission')
# 绘制
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
im0 = axes[0, 0].imshow(np.abs(E).T, origin='lower', cmap='hot')
axes[0, 0].set_title('Electric Field |E|', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
plt.colorbar(im0, ax=axes[0, 0])
im1 = axes[0, 1].imshow(np.angle(E).T, origin='lower', cmap='hsv')
axes[0, 1].set_title('Phase of E', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('y')
plt.colorbar(im1, ax=axes[0, 1])
im2 = axes[1, 0].imshow(sens_energy.T, origin='lower', cmap='RdBu_r')
axes[1, 0].set_title('Sensitivity (Energy)', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('x')
axes[1, 0].set_ylabel('y')
plt.colorbar(im2, ax=axes[1, 0])
im3 = axes[1, 1].imshow(sens_trans.T, origin='lower', cmap='RdBu_r')
axes[1, 1].set_title('Sensitivity (Transmission)', fontsize=12, fontweight='bold')
axes[1, 1].set_xlabel('x')
axes[1, 1].set_ylabel('y')
plt.colorbar(im3, ax=axes[1, 1])
plt.tight_layout()
plt.savefig(f'{output_dir}/adjoint_sensitivity.png', dpi=150, bbox_inches='tight')
plt.close()
demo_adjoint_sensitivity()
print("✓ 伴随法灵敏度分析图已保存")
# =============================================================================
# 第四部分：优化算法
# =============================================================================
print("\n[4] 优化算法...")
class OptimalityCriteria:
"""最优性准则（OC）方法"""
def __init__(self, volfrac, move=0.2):
"""
参数:
volfrac: 目标体积分数
move: 密度变化限制
"""
self.volfrac = volfrac
self.move = move
def update(self, rho, dc, dv):
"""
OC更新
参数:
rho: 当前密度
dc: 目标函数灵敏度
dv: 体积约束灵敏度
返回:
rho_new: 更新后的密度
"""
n = rho.size
rho_flat = rho.flatten()
dc_flat = dc.flatten()
dv_flat = dv.flatten()
# 二分法寻找拉格朗日乘子
l1, l2 = 1e-9, 1e9
for _ in range(100): # 最大迭代次数
lmid = 0.5 * (l1 + l2)
# OC更新公式
rho_new = np.maximum(0.0, np.maximum(rho_flat - self.move,
np.minimum(1.0, np.minimum(rho_flat + self.move,
rho_flat * np.sqrt(-dc_flat / (dv_flat * lmid + 1e-10))))))
# 检查体积约束
if np.sum(rho_new) > self.volfrac * n:
l1 = lmid
else:
l2 = lmid
if (l2 - l1) / (l1 + l2 + 1e-10) < 1e-4:
break
return rho_new.reshape(rho.shape)
class MethodOfMovingAsymptotes:
"""移动渐近线方法（MMA）"""
def __init__(self, n, m, move=0.2):
"""
参数:
n: 设计变量数
m: 约束数
move: 移动限制
"""
self.n = n
self.m = m
self.move = move
# MMA参数
self.asy_init = 0.5
self.asy_incr = 1.2
self.asy_decr = 0.7
# 渐近线
self.low = None
self.upp = None
# 历史值
self.xold1 = None
self.xold2 = None
def update_asymptotes(self, x, xold1, xold2):
"""更新渐近线"""
if self.low is None:
self.low = x - self.asy_init
self.upp = x + self.asy_init
else:
# 根据变化趋势调整
for i in range(self.n):
if x[i] > xold1[i] and xold1[i] > xold2[i]:
# 同向变化，扩大渐近线
self.low[i] = x[i] - self.asy_incr * (xold1[i] - self.low[i])
self.upp[i] = x[i] + self.asy_incr * (self.upp[i] - xold1[i])
elif x[i] < xold1[i] and xold1[i] < xold2[i]:
# 反向变化，缩小渐近线
self.low[i] = x[i] - self.asy_decr * (xold1[i] - self.low[i])
self.upp[i] = x[i] + self.asy_decr * (self.upp[i] - xold1[i])
else:
# 混合变化
self.low[i] = x[i] - (xold1[i] - self.low[i])
self.upp[i] = x[i] + (self.upp[i] - xold1[i])
# 限制渐近线范围
self.low = np.maximum(self.low, x - 10)
self.low = np.minimum(self.low, x - 0.01)
self.upp = np.minimum(self.upp, x + 10)
self.upp = np.maximum(self.upp, x + 0.01)
def optimize(self, x, f, df, g, dg, xmin, xmax):
"""
MMA优化步骤
参数:
x: 当前设计变量
f: 目标函数值
df: 目标函数梯度
g: 约束函数值
dg: 约束函数梯度
xmin, xmax: 变量边界
返回:
x_new: 新设计变量
"""
# 更新历史
if self.xold1 is not None:
self.xold2 = self.xold1.copy()
self.xold1 = x.copy()
# 更新渐近线
if self.xold1 is not None and self.xold2 is not None:
self.update_asymptotes(x, self.xold1, self.xold2)
# 简化：使用OC方法作为替代
# 实际MMA实现较复杂，这里使用简化版本
oc = OptimalityCriteria(volfrac=np.mean(x), move=self.move)
dv = np.ones_like(x)
x_new = oc.update(x.reshape(1, -1), df, dv.reshape(1, -1))
return x_new.flatten()
# =============================================================================
# 第五部分：波导拓扑优化
# =============================================================================
print("\n[5] 波导拓扑优化...")
class WaveguideOptimizer:
"""波导拓扑优化器"""
def __init__(self, nx, ny, target_frequency, volfrac=0.4):
"""
参数:
nx, ny: 设计域网格数
target_frequency: 目标频率 (Hz)
volfrac: 目标体积分数
"""
self.nx = nx
self.ny = ny
self.target_freq = target_frequency
self.volfrac = volfrac
# 物理参数
self.c = 3e8
self.omega = 2 * np.pi * target_frequency
self.wavelength = self.c / target_frequency
# 网格间距
self.dx = self.wavelength / 20
self.dy = self.wavelength / 20
# SIMP模型
self.simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=12.0, p=3.0)
# 过滤
rmin = 3.0
self.filter = DensityFilter(nx, ny, rmin)
self.projection = DensityProjection(beta=1.0, eta=0.5)
# 优化器
self.oc = OptimalityCriteria(volfrac, move=0.2)
def initialize(self):
"""初始化密度场"""
return np.ones((self.nx, self.ny)) * self.volfrac
def analyze(self, rho):
"""
分析当前设计
返回:
objective: 目标函数值（传输系数）
constraint: 约束函数值（体积）
"""
# 过滤和投影
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
# 计算介电常数
epsilon = self.simp.interpolate(rho_phys)
# 简化：计算传输系数（基于波导模式匹配）
# 实际应求解完整的电磁场
# 模拟传输系数计算
# 传输与密度分布的平滑度相关
transmission = self._compute_transmission(epsilon)
# 目标函数：最大化传输（最小化负传输）
objective = -transmission
# 体积约束
constraint = np.mean(rho_phys) - self.volfrac
return objective, constraint, rho_phys
def _compute_transmission(self, epsilon):
"""计算传输系数（简化模型）"""
# 简化的传输模型
# 传输与介电常数分布的均匀性和连通性相关
# 计算介电常数的空间变化
grad_x = np.gradient(epsilon, axis=0)
grad_y = np.gradient(epsilon, axis=1)
grad_mag = np.sqrt(grad_x**2 + grad_y**2)
# 平滑度越高，传输越好
smoothness = 1.0 / (1.0 + np.mean(grad_mag))
# 检查连通性（简化）
center_line = epsilon[self.nx//2, :]
connectivity = np.mean(center_line > np.mean(epsilon))
transmission = 0.5 * smoothness + 0.5 * connectivity
return transmission
def compute_sensitivities(self, rho):
"""计算灵敏度"""
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
epsilon = self.simp.interpolate(rho_phys)
# 数值微分计算灵敏度（使用更稳定的采样方法）
dc = np.zeros_like(rho)
eps = 1e-4
obj0, _, _ = self.analyze(rho)
# 采样计算灵敏度以加速
sample_step = max(1, self.nx // 10)
for i in range(0, self.nx, sample_step):
for j in range(0, self.ny, sample_step):
rho_perturbed = rho.copy()
rho_perturbed[i, j] = min(rho_perturbed[i, j] + eps, 1.0)
obj_perturbed, _, _ = self.analyze(rho_perturbed)
if not np.isnan(obj_perturbed) and not np.isnan(obj0):
dc[i, j] = (obj_perturbed - obj0) / eps
# 插值填充
from scipy.ndimage import zoom
if sample_step > 1:
dc = zoom(dc, (self.nx/dc.shape[0], self.ny/dc.shape[1]), order=1)
dc = dc[:self.nx, :self.ny] # 确保尺寸匹配
# 过滤灵敏度
dc_filtered = self.filter.filter_sensitivity(dc, rho)
# 投影灵敏度链式法则
dproj_drho = self.projection.derivative(rho_filtered)
dc_phys = dc_filtered * dproj_drho
# 处理NaN值
dc_phys = np.nan_to_num(dc_phys, nan=0.0, posinf=0.0, neginf=0.0)
# 体积约束灵敏度
dv = np.ones_like(rho) / (self.nx * self.ny)
return dc_phys, dv
def optimize(self, n_iterations=50):
"""运行优化"""
rho = self.initialize()
history = {
'objective': [],
'volume': [],
'transmission': []
}
for iter in range(n_iterations):
# 分析
obj, vol, rho_phys = self.analyze(rho)
# 计算灵敏度
dc, dv = self.compute_sensitivities(rho)
# OC更新
rho_new = self.oc.update(rho, dc, dv)
# 投影锐度更新
if iter > 0 and iter % 10 == 0:
self.projection.beta = min(self.projection.beta * 1.2, 32.0)
# 记录历史
history['objective'].append(obj)
history['volume'].append(np.mean(rho_phys))
history['transmission'].append(-obj)
# 打印进度
if iter % 10 == 0:
print(f" Iter {iter}: Obj={obj:.4f}, Vol={np.mean(rho_phys):.4f}, Trans={-obj:.4f}")
rho = rho_new
return rho, history
# 运行波导优化
print(" 运行波导拓扑优化...")
wave_opt = WaveguideOptimizer(nx=60, ny=40, target_frequency=3e9, volfrac=0.4)
rho_final, history = wave_opt.optimize(n_iterations=30)
# 绘制波导优化结果
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 最终设计
im0 = axes[0, 0].imshow(rho_final.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 0].set_title('Optimized Waveguide Design', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
plt.colorbar(im0, ax=axes[0, 0])
# 目标函数历史
axes[0, 1].plot(history['objective'], 'b-', linewidth=2)
axes[0, 1].set_xlabel('Iteration', fontsize=11)
axes[0, 1].set_ylabel('Objective Function', fontsize=11)
axes[0, 1].set_title('Optimization Convergence', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3)
# 体积历史
axes[1, 0].plot(history['volume'], 'r-', linewidth=2)
axes[1, 0].axhline(y=wave_opt.volfrac, color='k', linestyle='--', label='Target')
axes[1, 0].set_xlabel('Iteration', fontsize=11)
axes[1, 0].set_ylabel('Volume Fraction', fontsize=11)
axes[1, 0].set_title('Volume Constraint', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# 传输系数历史
axes[1, 1].plot(history['transmission'], 'g-', linewidth=2)
axes[1, 1].set_xlabel('Iteration', fontsize=11)
axes[1, 1].set_ylabel('Transmission Coefficient', fontsize=11)
axes[1, 1].set_title('Transmission Performance', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/waveguide_optimization.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 波导拓扑优化图已保存")
# =============================================================================
# 第六部分：天线结构优化
# =============================================================================
print("\n[6] 天线结构拓扑优化...")
class AntennaOptimizer:
"""天线结构拓扑优化"""
def __init__(self, nx, ny, freq_min, freq_max, volfrac=0.3):
"""
参数:
nx, ny: 设计域网格数
freq_min, freq_max: 频率范围
volfrac: 目标体积分数
"""
self.nx = nx
self.ny = ny
self.freq_min = freq_min
self.freq_max = freq_max
self.volfrac = volfrac
# 频率采样点
self.freqs = np.linspace(freq_min, freq_max, 5)
# SIMP模型
self.simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=4.0, p=3.0)
# 过滤
self.filter = DensityFilter(nx, ny, rmin=2.5)
self.projection = DensityProjection(beta=1.0, eta=0.5)
# 优化器
self.oc = OptimalityCriteria(volfrac, move=0.2)
def initialize(self):
"""初始化 - 从中心开始"""
rho = np.zeros((self.nx, self.ny))
cx, cy = self.nx // 2, self.ny // 2
# 创建初始种子
for i in range(self.nx):
for j in range(self.ny):
dist = np.sqrt((i - cx)**2 + (j - cy)**2)
if dist < min(self.nx, self.ny) // 4:
rho[i, j] = self.volfrac
return rho
def compute_bandwidth(self, rho):
"""计算带宽性能（简化模型）"""
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
epsilon = self.simp.interpolate(rho_phys)
# 简化的带宽评估
# 带宽与结构的对称性和连通性相关
# 对称性
left = epsilon[:, :self.ny//2]
right = epsilon[:, self.ny//2:]
if right.shape[1] > left.shape[1]:
right = right[:, :-1]
symmetry = 1.0 - np.mean(np.abs(left - right[:, ::-1])) / np.mean(epsilon)
# 连通性
center_mass = np.mean(np.where(epsilon > np.mean(epsilon)))
connectivity = 1.0 - abs(center_mass - self.nx//2) / (self.nx//2)
bandwidth = 0.6 * symmetry + 0.4 * connectivity
return bandwidth
def analyze(self, rho):
"""分析天线性能"""
bandwidth = self.compute_bandwidth(rho)
# 目标：最大化带宽
objective = -bandwidth
# 体积约束
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
constraint = np.mean(rho_phys) - self.volfrac
return objective, constraint, rho_phys
def optimize(self, n_iterations=40):
"""运行优化"""
rho = self.initialize()
history = {
'bandwidth': [],
'volume': []
}
for iter in range(n_iterations):
obj, vol, rho_phys = self.analyze(rho)
# 数值灵敏度
dc = np.zeros_like(rho)
eps = 1e-6
obj0 = obj
for i in range(0, self.nx, 3): # 采样计算以加速
for j in range(0, self.ny, 3):
rho_perturbed = rho.copy()
rho_perturbed[i, j] = min(rho_perturbed[i, j] + eps, 1.0)
obj_perturbed, _, _ = self.analyze(rho_perturbed)
dc[i, j] = (obj_perturbed - obj0) / eps
# 插值填充
from scipy.ndimage import zoom
if self.nx > 10 and self.ny > 10:
dc = zoom(dc, (self.nx/dc.shape[0], self.ny/dc.shape[1]), order=1)
dv = np.ones_like(rho) / (self.nx * self.ny)
# OC更新
rho_new = self.oc.update(rho, dc, dv)
# 更新投影锐度
if iter > 0 and iter % 10 == 0:
self.projection.beta = min(self.projection.beta * 1.3, 32.0)
history['bandwidth'].append(-obj)
history['volume'].append(np.mean(rho_phys))
if iter % 10 == 0:
print(f" Iter {iter}: BW={-obj:.4f}, Vol={np.mean(rho_phys):.4f}")
rho = rho_new
return rho, history
# 运行天线优化
print(" 运行天线拓扑优化...")
ant_opt = AntennaOptimizer(nx=50, ny=50, freq_min=2e9, freq_max=6e9, volfrac=0.3)
rho_antenna, history_ant = ant_opt.optimize(n_iterations=30)
# 绘制天线优化结果
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 最终设计
im0 = axes[0, 0].imshow(rho_antenna.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 0].set_title('Optimized Antenna Structure', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
plt.colorbar(im0, ax=axes[0, 0])
# 初始设计对比
rho_init = ant_opt.initialize()
im1 = axes[0, 1].imshow(rho_init.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 1].set_title('Initial Design', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('y')
plt.colorbar(im1, ax=axes[0, 1])
# 带宽收敛
axes[1, 0].plot(history_ant['bandwidth'], 'b-', linewidth=2)
axes[1, 0].set_xlabel('Iteration', fontsize=11)
axes[1, 0].set_ylabel('Bandwidth Metric', fontsize=11)
axes[1, 0].set_title('Bandwidth Optimization', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3)
# 体积收敛
axes[1, 1].plot(history_ant['volume'], 'r-', linewidth=2)
axes[1, 1].axhline(y=ant_opt.volfrac, color='k', linestyle='--', label='Target')
axes[1, 1].set_xlabel('Iteration', fontsize=11)
axes[1, 1].set_ylabel('Volume Fraction', fontsize=11)
axes[1, 1].set_title('Volume Constraint', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/antenna_optimization.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 天线拓扑优化图已保存")
12. 结果分析与讨论
12.1 SIMP插值模型分析
从SIMP插值曲线可以看出：
- 当惩罚因子p=1时，插值呈线性关系，无法有效抑制中间密度
- 当p=3时，中间密度的材料属性被显著惩罚，有利于获得0/1解
- 当p>3时，惩罚过强可能导致优化收敛困难
12.2 密度过滤效果
密度过滤有效消除了棋盘格现象：
- 原始密度场呈现明显的棋盘格模式
- 经过过滤后，密度分布变得平滑
- 投影后获得接近0/1的清晰结构
12.3 优化收敛性
波导和天线优化案例显示：
- 目标函数在迭代初期快速下降
- 体积约束得到有效满足
- 优化通常在30-50次迭代后收敛
12.4 多目标权衡
Pareto前沿分析表明：
- 不同权重组合产生不同的优化设计
- 连通性和平滑度之间存在权衡关系
- 工程师可以根据实际需求选择合适的设计方案
13. 拓扑优化的前沿发展
13.1 多尺度拓扑优化
多尺度拓扑优化同时优化宏观结构和微观结构，实现材料-结构一体化设计。这种方法可以：
- 获得更优异的性能
- 实现超材料设计
- 满足多物理场需求
13.2 鲁棒拓扑优化
考虑不确定性的鲁棒拓扑优化：
- 制造不确定性
- 材料参数不确定性
- 载荷不确定性
13.3 机器学习辅助优化
机器学习技术在拓扑优化中的应用：
- 代理模型加速
- 生成式设计
- 实时优化
14. 总结
本教程系统介绍了电磁场仿真中的拓扑优化技术，包括：
- 理论基础：SIMP材料插值模型、密度过滤与投影、伴随法灵敏度分析
- 优化算法：最优性准则方法（OC）、移动渐近线方法（MMA）
- 应用案例：波导优化、天线结构优化、多目标优化
- 工程实践：数值实现流程、收敛准则、制造约束处理
通过本教程的学习，读者应能够：
- 理解拓扑优化的基本原理
- 掌握SIMP方法和密度过滤技术
- 实现伴随法灵敏度分析
- 开展电磁结构的拓扑优化设计
拓扑优化作为结构优化设计的前沿技术，在电磁工程领域具有广阔的应用前景。随着计算能力的提升和算法的改进，拓扑优化将在天线设计、微波器件、超材料等领域发挥越来越重要的作用。
参考文献
-
Bendsøe, M. P., & Sigmund, O. (2003). Topology Optimization: Theory, Methods, and Applications. Springer.
-
Sigmund, O. (2007). Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization, 33(4-5), 401-424.
-
Guest, J. K., Prévost, J. H., & Belytschko, T. (2004). Achieving minimum length scale in topology optimization using nodal design variables and projection functions. International Journal for Numerical Methods in Engineering, 61(2), 238-254.
-
Jensen, J. S., & Sigmund, O. (2011). Topology optimization for nano-photonics. Laser & Photonics Reviews, 5(2), 308-321.
-
Aage, N., Andreassen, E., & Lazarov, B. S. (2015). Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework. Structural and Multidisciplinary Optimization, 51(3), 565-572.
附录：完整代码清单
# -*- coding: utf-8 -*-
"""
主题077：拓扑优化仿真
Topology Optimization for Electromagnetic Applications
本程序实现电磁场仿真中的拓扑优化算法，包括：
1. SIMP材料插值模型
2. 密度过滤与投影
3. 伴随法灵敏度分析
4. 优化算法（OC、MMA）
5. 电磁波导优化
6. 天线结构优化
7. 多目标拓扑优化
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.sparse import diags, csr_matrix
from scipy.sparse.linalg import spsolve
from scipy.ndimage import gaussian_filter
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("=" * 60)
print("拓扑优化仿真程序")
print("=" * 60)
# =============================================================================
# 第一部分：SIMP材料插值模型
# =============================================================================
print("\n[1] SIMP材料插值模型...")
class SIMPMaterial:
"""SIMP (Solid Isotropic Material with Penalization) 材料模型"""
def __init__(self, epsilon_min=1.0, epsilon_max=10.0, p=3.0):
"""
参数:
epsilon_min: 空材料介电常数 (通常接近1)
epsilon_max: 实体材料介电常数
p: 惩罚因子
"""
self.epsilon_min = epsilon_min
self.epsilon_max = epsilon_max
self.p = p
def interpolate(self, rho):
"""
SIMP插值
参数:
rho: 密度场 (0到1之间)
返回:
epsilon: 等效介电常数
"""
return self.epsilon_min + (self.epsilon_max - self.epsilon_min) * rho**self.p
def derivative(self, rho):
"""
SIMP插值对密度的导数
参数:
rho: 密度场
返回:
depsilon_drho: 导数
"""
return self.p * (self.epsilon_max - self.epsilon_min) * rho**(self.p - 1)
# 演示SIMP插值
simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=10.0, p=3.0)
rho_demo = np.linspace(0, 1, 100)
epsilon_demo = simp.interpolate(rho_demo)
def plot_simp_model():
"""绘制SIMP模型"""
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# SIMP插值曲线
for p in [1, 2, 3, 4, 5]:
simp_temp = SIMPMaterial(epsilon_min=1.0, epsilon_max=10.0, p=p)
epsilon_temp = simp_temp.interpolate(rho_demo)
axes[0].plot(rho_demo, epsilon_temp, linewidth=2, label=f'p={p}')
axes[0].axhline(y=1.0, color='k', linestyle='--', alpha=0.3)
axes[0].axhline(y=10.0, color='k', linestyle='--', alpha=0.3)
axes[0].set_xlabel('Density (rho)', fontsize=11)
axes[0].set_ylabel('Relative Permittivity', fontsize=11)
axes[0].set_title('SIMP Interpolation Model', fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 惩罚效果示意
rho_binary = np.where(rho_demo > 0.5, 1.0, 0.0)
axes[1].plot(rho_demo, rho_demo, 'b--', linewidth=2, label='Linear (p=1)')
axes[1].plot(rho_demo, rho_demo**3, 'r-', linewidth=2, label='SIMP (p=3)')
axes[1].plot(rho_demo, rho_binary, 'g:', linewidth=2, label='Binary')
axes[1].set_xlabel('Physical Density', fontsize=11)
axes[1].set_ylabel('Effective Density', fontsize=11)
axes[1].set_title('Penalization Effect', fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/simp_model.png', dpi=150, bbox_inches='tight')
plt.close()
plot_simp_model()
print("✓ SIMP模型图已保存")
# =============================================================================
# 第二部分：密度过滤与投影
# =============================================================================
print("\n[2] 密度过滤与投影...")
class DensityFilter:
"""密度过滤器 - 消除棋盘格现象"""
def __init__(self, nx, ny, rmin):
"""
参数:
nx, ny: 网格数
rmin: 过滤半径（网格单元数）
"""
self.nx = nx
self.ny = ny
self.rmin = rmin
# 构建过滤权重矩阵
self.H, self.Hs = self._build_filter_matrix()
def _build_filter_matrix(self):
"""构建过滤矩阵"""
n = self.nx * self.ny
H = sparse.lil_matrix((n, n))
for i in range(self.nx):
for j in range(self.ny):
e1 = i * self.ny + j
for k in range(max(0, i - int(self.rmin)), min(self.nx, i + int(self.rmin) + 1)):
for l in range(max(0, j - int(self.rmin)), min(self.ny, j + int(self.rmin) + 1)):
e2 = k * self.ny + l
dist = np.sqrt((i - k)**2 + (j - l)**2)
H[e1, e2] = max(0, self.rmin - dist)
H = H.tocsr()
Hs = np.array(H.sum(axis=1)).flatten()
return H, Hs
def filter_density(self, rho):
"""过滤密度场"""
rho_flat = rho.flatten()
rho_filtered = (self.H @ rho_flat) / self.Hs
return np.array(rho_filtered).flatten().reshape(self.nx, self.ny)
def filter_sensitivity(self, d_obj_drho_filtered, rho):
"""过滤灵敏度"""
n = self.nx * self.ny
d_obj_drho = np.zeros(n)
d_obj_drho_filtered_flat = d_obj_drho_filtered.flatten()
for i in range(n):
for j in range(n):
if self.H[i, j] > 0:
d_obj_drho[j] += float(self.H[i, j]) * d_obj_drho_filtered_flat[i] / self.Hs[i]
return d_obj_drho.reshape(self.nx, self.ny)
class DensityProjection:
"""密度投影 - 获得清晰的0/1解"""
def __init__(self, beta=1.0, eta=0.5):
"""
参数:
beta: 投影锐度参数
eta: 投影阈值
"""
self.beta = beta
self.eta = eta
def project(self, rho_tilde):
"""
Heaviside投影
参数:
rho_tilde: 过滤后的密度
返回:
rho: 投影后的密度
"""
return (np.tanh(self.beta * self.eta) + np.tanh(self.beta * (rho_tilde - self.eta))) / \
(np.tanh(self.beta * self.eta) + np.tanh(self.beta * (1 - self.eta)))
def derivative(self, rho_tilde):
"""投影函数的导数"""
numerator = self.beta * (1 - np.tanh(self.beta * (rho_tilde - self.eta))**2)
denominator = np.tanh(self.beta * self.eta) + np.tanh(self.beta * (1 - self.eta))
return numerator / denominator
# 演示过滤与投影
def demo_filter_projection():
"""演示密度过滤和投影效果"""
nx, ny = 60, 40
# 创建带有棋盘格的初始密度场
np.random.seed(42)
rho_raw = np.random.rand(nx, ny)
# 添加棋盘格模式
for i in range(nx):
for j in range(ny):
if (i + j) % 2 == 0:
rho_raw[i, j] = 0.9
else:
rho_raw[i, j] = 0.1
# 应用过滤
df = DensityFilter(nx, ny, rmin=3.0)
rho_filtered = df.filter_density(rho_raw)
# 应用投影
proj = DensityProjection(beta=8.0, eta=0.5)
rho_projected = proj.project(rho_filtered)
# 绘制
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
im0 = axes[0, 0].imshow(rho_raw.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 0].set_title('Raw Density (with checkerboard)', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
plt.colorbar(im0, ax=axes[0, 0])
im1 = axes[0, 1].imshow(rho_filtered.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 1].set_title('Filtered Density', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('y')
plt.colorbar(im1, ax=axes[0, 1])
im2 = axes[1, 0].imshow(rho_projected.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[1, 0].set_title('Projected Density (beta=8)', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('x')
axes[1, 0].set_ylabel('y')
plt.colorbar(im2, ax=axes[1, 0])
# 投影函数曲线
rho_test = np.linspace(0, 1, 100)
for beta in [1, 4, 8, 16]:
proj_temp = DensityProjection(beta=beta, eta=0.5)
rho_proj_test = proj_temp.project(rho_test)
axes[1, 1].plot(rho_test, rho_proj_test, linewidth=2, label=f'beta={beta}')
axes[1, 1].plot([0, 1], [0, 1], 'k--', alpha=0.3)
axes[1, 1].set_xlabel('Filtered Density', fontsize=11)
axes[1, 1].set_ylabel('Projected Density', fontsize=11)
axes[1, 1].set_title('Projection Function', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/density_filter_projection.png', dpi=150, bbox_inches='tight')
plt.close()
demo_filter_projection()
print("✓ 密度过滤与投影图已保存")
# =============================================================================
# 第三部分：伴随法灵敏度分析
# =============================================================================
print("\n[3] 伴随法灵敏度分析...")
class AdjointSensitivity:
"""伴随法灵敏度分析"""
def __init__(self, nx, ny, dx, dy, omega, epsilon):
"""
参数:
nx, ny: 网格数
dx, dy: 网格间距
omega: 角频率
epsilon: 介电常数分布
"""
self.nx = nx
self.ny = ny
self.dx = dx
self.dy = dy
self.omega = omega
self.epsilon = epsilon
self.k0 = omega / 3e8 # 自由空间波数
def build_fem_matrix(self, epsilon_field):
"""构建FEM系统矩阵"""
n = self.nx * self.ny
# 构建离散的Helmholtz算子
# (nabla^2 + k^2) E = 0
row = []
col = []
data = []
for i in range(self.nx):
for j in range(self.ny):
idx = i * self.ny + j
# 对角元素
k2 = (self.k0**2) * epsilon_field[i, j]
coeff = 2/self.dx**2 + 2/self.dy**2 - k2
row.append(idx)
col.append(idx)
data.append(coeff)
# x方向邻居
if i > 0:
row.append(idx)
col.append((i-1) * self.ny + j)
data.append(-1/self.dx**2)
if i < self.nx - 1:
row.append(idx)
col.append((i+1) * self.ny + j)
data.append(-1/self.dx**2)
# y方向邻居
if j > 0:
row.append(idx)
col.append(i * self.ny + (j-1))
data.append(-1/self.dy**2)
if j < self.ny - 1:
row.append(idx)
col.append(i * self.ny + (j+1))
data.append(-1/self.dy**2)
A = csr_matrix((data, (row, col)), shape=(n, n))
return A
def solve_field(self, epsilon_field, source):
"""求解电磁场"""
A = self.build_fem_matrix(epsilon_field)
n = self.nx * self.ny
# 添加简单的边界条件（Dirichlet）
# 简化处理：直接求解
E = spsolve(A, source.flatten())
return E.reshape(self.nx, self.ny)
def compute_sensitivity(self, epsilon_field, E, objective_type='energy'):
"""
计算灵敏度
参数:
epsilon_field: 介电常数场
E: 电场分布
objective_type: 目标函数类型
返回:
sensitivity: 灵敏度场
"""
if objective_type == 'energy':
# 目标函数：电场能量
# J = integral(|E|^2)
# dJ/depsilon = -omega^2 * |E|^2
sensitivity = -self.omega**2 * np.abs(E)**2
elif objective_type == 'transmission':
# 目标函数：传输系数
# 需要伴随场
n = self.nx * self.ny
# 构建伴随源（在输出端）
adj_source = np.zeros((self.nx, self.ny))
adj_source[-5:, self.ny//2-2:self.ny//2+2] = 1.0
# 求解伴随场
A = self.build_fem_matrix(epsilon_field)
lambda_field = spsolve(A.T, adj_source.flatten())
lambda_field = lambda_field.reshape(self.nx, self.ny)
# 灵敏度：dJ/depsilon = -omega^2 * E * lambda
sensitivity = -self.omega**2 * np.real(E * lambda_field)
else:
sensitivity = np.zeros_like(epsilon_field)
return sensitivity
# 演示伴随法灵敏度
def demo_adjoint_sensitivity():
"""演示伴随法灵敏度计算"""
nx, ny = 40, 30
dx = dy = 0.01
omega = 2 * np.pi * 3e9 # 3 GHz
# 初始密度场
rho = np.ones((nx, ny)) * 0.5
simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=10.0, p=3.0)
epsilon_field = simp.interpolate(rho)
# 创建伴随灵敏度分析器
adj = AdjointSensitivity(nx, ny, dx, dy, omega, epsilon_field)
# 创建源
source = np.zeros((nx, ny))
source[5, ny//2-2:ny//2+2] = 1.0 # 输入端口
# 求解场
E = adj.solve_field(epsilon_field, source)
# 计算灵敏度
sens_energy = adj.compute_sensitivity(epsilon_field, E, 'energy')
sens_trans = adj.compute_sensitivity(epsilon_field, E, 'transmission')
# 绘制
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
im0 = axes[0, 0].imshow(np.abs(E).T, origin='lower', cmap='hot')
axes[0, 0].set_title('Electric Field |E|', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
plt.colorbar(im0, ax=axes[0, 0])
im1 = axes[0, 1].imshow(np.angle(E).T, origin='lower', cmap='hsv')
axes[0, 1].set_title('Phase of E', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('y')
plt.colorbar(im1, ax=axes[1, 1])
im2 = axes[1, 0].imshow(sens_energy.T, origin='lower', cmap='RdBu_r')
axes[1, 0].set_title('Sensitivity (Energy)', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('x')
axes[1, 0].set_ylabel('y')
plt.colorbar(im2, ax=axes[1, 0])
im3 = axes[1, 1].imshow(sens_trans.T, origin='lower', cmap='RdBu_r')
axes[1, 1].set_title('Sensitivity (Transmission)', fontsize=12, fontweight='bold')
axes[1, 1].set_xlabel('x')
axes[1, 1].set_ylabel('y')
plt.colorbar(im3, ax=axes[1, 1])
plt.tight_layout()
plt.savefig(f'{output_dir}/adjoint_sensitivity.png', dpi=150, bbox_inches='tight')
plt.close()
demo_adjoint_sensitivity()
print("✓ 伴随法灵敏度分析图已保存")
# =============================================================================
# 第四部分：优化算法
# =============================================================================
print("\n[4] 优化算法...")
class OptimalityCriteria:
"""最优性准则（OC）方法"""
def __init__(self, volfrac, move=0.2):
"""
参数:
volfrac: 目标体积分数
move: 密度变化限制
"""
self.volfrac = volfrac
self.move = move
def update(self, rho, dc, dv):
"""
OC更新
参数:
rho: 当前密度
dc: 目标函数灵敏度
dv: 体积约束灵敏度
返回:
rho_new: 更新后的密度
"""
n = rho.size
rho_flat = rho.flatten()
dc_flat = dc.flatten()
dv_flat = dv.flatten()
# 二分法寻找拉格朗日乘子
l1, l2 = 1e-9, 1e9
for _ in range(100): # 最大迭代次数
lmid = 0.5 * (l1 + l2)
# OC更新公式
rho_new = np.maximum(0.0, np.maximum(rho_flat - self.move,
np.minimum(1.0, np.minimum(rho_flat + self.move,
rho_flat * np.sqrt(-dc_flat / (dv_flat * lmid + 1e-10))))))
# 检查体积约束
if np.sum(rho_new) > self.volfrac * n:
l1 = lmid
else:
l2 = lmid
if (l2 - l1) / (l1 + l2 + 1e-10) < 1e-4:
break
return rho_new.reshape(rho.shape)
class MethodOfMovingAsymptotes:
"""移动渐近线方法（MMA）"""
def __init__(self, n, m, move=0.2):
"""
参数:
n: 设计变量数
m: 约束数
move: 移动限制
"""
self.n = n
self.m = m
self.move = move
# MMA参数
self.asy_init = 0.5
self.asy_incr = 1.2
self.asy_decr = 0.7
# 渐近线
self.low = None
self.upp = None
# 历史值
self.xold1 = None
self.xold2 = None
def update_asymptotes(self, x, xold1, xold2):
"""更新渐近线"""
if self.low is None:
self.low = x - self.asy_init
self.upp = x + self.asy_init
else:
# 根据变化趋势调整
for i in range(self.n):
if x[i] > xold1[i] and xold1[i] > xold2[i]:
# 同向变化，扩大渐近线
self.low[i] = x[i] - self.asy_incr * (xold1[i] - self.low[i])
self.upp[i] = x[i] + self.asy_incr * (self.upp[i] - xold1[i])
elif x[i] < xold1[i] and xold1[i] < xold2[i]:
# 反向变化，缩小渐近线
self.low[i] = x[i] - self.asy_decr * (xold1[i] - self.low[i])
self.upp[i] = x[i] + self.asy_decr * (self.upp[i] - xold1[i])
else:
# 混合变化
self.low[i] = x[i] - (xold1[i] - self.low[i])
self.upp[i] = x[i] + (self.upp[i] - xold1[i])
# 限制渐近线范围
self.low = np.maximum(self.low, x - 10)
self.low = np.minimum(self.low, x - 0.01)
self.upp = np.minimum(self.upp, x + 10)
self.upp = np.maximum(self.upp, x + 0.01)
def optimize(self, x, f, df, g, dg, xmin, xmax):
"""
MMA优化步骤
参数:
x: 当前设计变量
f: 目标函数值
df: 目标函数梯度
g: 约束函数值
dg: 约束函数梯度
xmin, xmax: 变量边界
返回:
x_new: 新设计变量
"""
# 更新历史
if self.xold1 is not None:
self.xold2 = self.xold1.copy()
self.xold1 = x.copy()
# 更新渐近线
if self.xold1 is not None and self.xold2 is not None:
self.update_asymptotes(x, self.xold1, self.xold2)
# 简化：使用OC方法作为替代
# 实际MMA实现较复杂，这里使用简化版本
oc = OptimalityCriteria(volfrac=np.mean(x), move=self.move)
dv = np.ones_like(x)
x_new = oc.update(x.reshape(1, -1), df, dv.reshape(1, -1))
return x_new.flatten()
# =============================================================================
# 第五部分：波导拓扑优化
# =============================================================================
print("\n[5] 波导拓扑优化...")
class WaveguideOptimizer:
"""波导拓扑优化器"""
def __init__(self, nx, ny, target_frequency, volfrac=0.4):
"""
参数:
nx, ny: 设计域网格数
target_frequency: 目标频率 (Hz)
volfrac: 目标体积分数
"""
self.nx = nx
self.ny = ny
self.target_freq = target_frequency
self.volfrac = volfrac
# 物理参数
self.c = 3e8
self.omega = 2 * np.pi * target_frequency
self.wavelength = self.c / target_frequency
# 网格间距
self.dx = self.wavelength / 20
self.dy = self.wavelength / 20
# SIMP模型
self.simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=12.0, p=3.0)
# 过滤
rmin = 3.0
self.filter = DensityFilter(nx, ny, rmin)
self.projection = DensityProjection(beta=1.0, eta=0.5)
# 优化器
self.oc = OptimalityCriteria(volfrac, move=0.2)
def initialize(self):
"""初始化密度场"""
return np.ones((self.nx, self.ny)) * self.volfrac
def analyze(self, rho):
"""
分析当前设计
返回:
objective: 目标函数值（传输系数）
constraint: 约束函数值（体积）
"""
# 过滤和投影
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
# 计算介电常数
epsilon = self.simp.interpolate(rho_phys)
# 简化：计算传输系数（基于波导模式匹配）
# 实际应求解完整的电磁场
# 模拟传输系数计算
# 传输与密度分布的平滑度相关
transmission = self._compute_transmission(epsilon)
# 目标函数：最大化传输（最小化负传输）
objective = -transmission
# 体积约束
constraint = np.mean(rho_phys) - self.volfrac
return objective, constraint, rho_phys
def _compute_transmission(self, epsilon):
"""计算传输系数（简化模型）"""
# 简化的传输模型
# 传输与介电常数分布的均匀性和连通性相关
# 计算介电常数的空间变化
grad_x = np.gradient(epsilon, axis=0)
grad_y = np.gradient(epsilon, axis=1)
grad_mag = np.sqrt(grad_x**2 + grad_y**2)
# 平滑度越高，传输越好
smoothness = 1.0 / (1.0 + np.mean(grad_mag))
# 检查连通性（简化）
center_line = epsilon[self.nx//2, :]
connectivity = np.mean(center_line > np.mean(epsilon))
transmission = 0.5 * smoothness + 0.5 * connectivity
return transmission
def compute_sensitivities(self, rho):
"""计算灵敏度"""
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
epsilon = self.simp.interpolate(rho_phys)
# 数值微分计算灵敏度（使用更稳定的采样方法）
dc = np.zeros_like(rho)
eps = 1e-4
obj0, _, _ = self.analyze(rho)
# 采样计算灵敏度以加速
sample_step = max(1, self.nx // 10)
for i in range(0, self.nx, sample_step):
for j in range(0, self.ny, sample_step):
rho_perturbed = rho.copy()
rho_perturbed[i, j] = min(rho_perturbed[i, j] + eps, 1.0)
obj_perturbed, _, _ = self.analyze(rho_perturbed)
if not np.isnan(obj_perturbed) and not np.isnan(obj0):
dc[i, j] = (obj_perturbed - obj0) / eps
# 插值填充
from scipy.ndimage import zoom
if sample_step > 1:
dc = zoom(dc, (self.nx/dc.shape[0], self.ny/dc.shape[1]), order=1)
dc = dc[:self.nx, :self.ny] # 确保尺寸匹配
# 过滤灵敏度
dc_filtered = self.filter.filter_sensitivity(dc, rho)
# 投影灵敏度链式法则
dproj_drho = self.projection.derivative(rho_filtered)
dc_phys = dc_filtered * dproj_drho
# 处理NaN值
dc_phys = np.nan_to_num(dc_phys, nan=0.0, posinf=0.0, neginf=0.0)
# 体积约束灵敏度
dv = np.ones_like(rho) / (self.nx * self.ny)
return dc_phys, dv
def optimize(self, n_iterations=50):
"""运行优化"""
rho = self.initialize()
history = {
'objective': [],
'volume': [],
'transmission': []
}
for iter in range(n_iterations):
# 分析
obj, vol, rho_phys = self.analyze(rho)
# 计算灵敏度
dc, dv = self.compute_sensitivities(rho)
# OC更新
rho_new = self.oc.update(rho, dc, dv)
# 投影锐度更新
if iter > 0 and iter % 10 == 0:
self.projection.beta = min(self.projection.beta * 1.2, 32.0)
# 记录历史
history['objective'].append(obj)
history['volume'].append(np.mean(rho_phys))
history['transmission'].append(-obj)
# 打印进度
if iter % 10 == 0:
print(f" Iter {iter}: Obj={obj:.4f}, Vol={np.mean(rho_phys):.4f}, Trans={-obj:.4f}")
rho = rho_new
return rho, history
# 运行波导优化
print(" 运行波导拓扑优化...")
wave_opt = WaveguideOptimizer(nx=60, ny=40, target_frequency=3e9, volfrac=0.4)
rho_final, history = wave_opt.optimize(n_iterations=30)
# 绘制波导优化结果
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 最终设计
im0 = axes[0, 0].imshow(rho_final.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 0].set_title('Optimized Waveguide Design', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
plt.colorbar(im0, ax=axes[0, 0])
# 目标函数历史
axes[0, 1].plot(history['objective'], 'b-', linewidth=2)
axes[0, 1].set_xlabel('Iteration', fontsize=11)
axes[0, 1].set_ylabel('Objective Function', fontsize=11)
axes[0, 1].set_title('Optimization Convergence', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3)
# 体积历史
axes[1, 0].plot(history['volume'], 'r-', linewidth=2)
axes[1, 0].axhline(y=wave_opt.volfrac, color='k', linestyle='--', label='Target')
axes[1, 0].set_xlabel('Iteration', fontsize=11)
axes[1, 0].set_ylabel('Volume Fraction', fontsize=11)
axes[1, 0].set_title('Volume Constraint', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# 传输系数历史
axes[1, 1].plot(history['transmission'], 'g-', linewidth=2)
axes[1, 1].set_xlabel('Iteration', fontsize=11)
axes[1, 1].set_ylabel('Transmission Coefficient', fontsize=11)
axes[1, 1].set_title('Transmission Performance', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/waveguide_optimization.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 波导拓扑优化图已保存")
# =============================================================================
# 第六部分：天线结构优化
# =============================================================================
print("\n[6] 天线结构拓扑优化...")
class AntennaOptimizer:
"""天线结构拓扑优化"""
def __init__(self, nx, ny, freq_min, freq_max, volfrac=0.3):
"""
参数:
nx, ny: 设计域网格数
freq_min, freq_max: 频率范围
volfrac: 目标体积分数
"""
self.nx = nx
self.ny = ny
self.freq_min = freq_min
self.freq_max = freq_max
self.volfrac = volfrac
# 频率采样点
self.freqs = np.linspace(freq_min, freq_max, 5)
# SIMP模型
self.simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=4.0, p=3.0)
# 过滤
self.filter = DensityFilter(nx, ny, rmin=2.5)
self.projection = DensityProjection(beta=1.0, eta=0.5)
# 优化器
self.oc = OptimalityCriteria(volfrac, move=0.2)
def initialize(self):
"""初始化 - 从中心开始"""
rho = np.zeros((self.nx, self.ny))
cx, cy = self.nx // 2, self.ny // 2
# 创建初始种子
for i in range(self.nx):
for j in range(self.ny):
dist = np.sqrt((i - cx)**2 + (j - cy)**2)
if dist < min(self.nx, self.ny) // 4:
rho[i, j] = self.volfrac
return rho
def compute_bandwidth(self, rho):
"""计算带宽性能（简化模型）"""
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
epsilon = self.simp.interpolate(rho_phys)
# 简化的带宽评估
# 带宽与结构的对称性和连通性相关
# 对称性
left = epsilon[:, :self.ny//2]
right = epsilon[:, self.ny//2:]
if right.shape[1] > left.shape[1]:
right = right[:, :-1]
symmetry = 1.0 - np.mean(np.abs(left - right[:, ::-1])) / np.mean(epsilon)
# 连通性
center_mass = np.mean(np.where(epsilon > np.mean(epsilon)))
connectivity = 1.0 - abs(center_mass - self.nx//2) / (self.nx//2)
bandwidth = 0.6 * symmetry + 0.4 * connectivity
return bandwidth
def analyze(self, rho):
"""分析天线性能"""
bandwidth = self.compute_bandwidth(rho)
# 目标：最大化带宽
objective = -bandwidth
# 体积约束
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
constraint = np.mean(rho_phys) - self.volfrac
return objective, constraint, rho_phys
def optimize(self, n_iterations=40):
"""运行优化"""
rho = self.initialize()
history = {
'bandwidth': [],
'volume': []
}
for iter in range(n_iterations):
obj, vol, rho_phys = self.analyze(rho)
# 数值灵敏度
dc = np.zeros_like(rho)
eps = 1e-6
obj0 = obj
for i in range(0, self.nx, 3): # 采样计算以加速
for j in range(0, self.ny, 3):
rho_perturbed = rho.copy()
rho_perturbed[i, j] = min(rho_perturbed[i, j] + eps, 1.0)
obj_perturbed, _, _ = self.analyze(rho_perturbed)
dc[i, j] = (obj_perturbed - obj0) / eps
# 插值填充
from scipy.ndimage import zoom
if self.nx > 10 and self.ny > 10:
dc = zoom(dc, (self.nx/dc.shape[0], self.ny/dc.shape[1]), order=1)
dv = np.ones_like(rho) / (self.nx * self.ny)
# OC更新
rho_new = self.oc.update(rho, dc, dv)
# 更新投影锐度
if iter > 0 and iter % 10 == 0:
self.projection.beta = min(self.projection.beta * 1.3, 32.0)
history['bandwidth'].append(-obj)
history['volume'].append(np.mean(rho_phys))
if iter % 10 == 0:
print(f" Iter {iter}: BW={-obj:.4f}, Vol={np.mean(rho_phys):.4f}")
rho = rho_new
return rho, history
# 运行天线优化
print(" 运行天线拓扑优化...")
ant_opt = AntennaOptimizer(nx=50, ny=50, freq_min=2e9, freq_max=6e9, volfrac=0.3)
rho_antenna, history_ant = ant_opt.optimize(n_iterations=30)
# 绘制天线优化结果
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 最终设计
im0 = axes[0, 0].imshow(rho_antenna.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 0].set_title('Optimized Antenna Structure', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
plt.colorbar(im0, ax=axes[0, 0])
# 初始设计对比
rho_init = ant_opt.initialize()
im1 = axes[0, 1].imshow(rho_init.T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[0, 1].set_title('Initial Design', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('y')
plt.colorbar(im1, ax=axes[0, 1])
# 带宽收敛
axes[1, 0].plot(history_ant['bandwidth'], 'b-', linewidth=2)
axes[1, 0].set_xlabel('Iteration', fontsize=11)
axes[1, 0].set_ylabel('Bandwidth Metric', fontsize=11)
axes[1, 0].set_title('Bandwidth Optimization', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3)
# 体积收敛
axes[1, 1].plot(history_ant['volume'], 'r-', linewidth=2)
axes[1, 1].axhline(y=ant_opt.volfrac, color='k', linestyle='--', label='Target')
axes[1, 1].set_xlabel('Iteration', fontsize=11)
axes[1, 1].set_ylabel('Volume Fraction', fontsize=11)
axes[1, 1].set_title('Volume Constraint', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/antenna_optimization.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 天线拓扑优化图已保存")
# =============================================================================
# 第七部分：多目标拓扑优化
# =============================================================================
print("\n[7] 多目标拓扑优化...")
class MultiObjectiveOptimizer:
"""多目标拓扑优化（加权和方法）"""
def __init__(self, nx, ny, volfrac=0.35):
self.nx = nx
self.ny = ny
self.volfrac = volfrac
self.simp = SIMPMaterial(epsilon_min=1.0, epsilon_max=8.0, p=3.0)
self.filter = DensityFilter(nx, ny, rmin=2.5)
self.projection = DensityProjection(beta=1.0, eta=0.5)
self.oc = OptimalityCriteria(volfrac, move=0.15)
def initialize(self):
return np.ones((self.nx, self.ny)) * self.volfrac
def objective_1(self, rho_phys):
"""目标1：最大化刚度/连通性"""
# 基于结构连通性的度量
from scipy.ndimage import label
labeled, num_features = label(rho_phys > 0.5)
if num_features == 0:
return 0.0
# 最大连通区域占比
sizes = [np.sum(labeled == i) for i in range(1, num_features + 1)]
max_size = max(sizes) if sizes else 0
connectivity = max_size / (self.nx * self.ny)
return connectivity
def objective_2(self, rho_phys):
"""目标2：最小化介电常数变化梯度（平滑度）"""
epsilon = self.simp.interpolate(rho_phys)
grad_x = np.gradient(epsilon, axis=0)
grad_y = np.gradient(epsilon, axis=1)
smoothness = 1.0 / (1.0 + np.mean(grad_x**2 + grad_y**2))
return smoothness
def combined_objective(self, rho, w1=0.5, w2=0.5):
"""组合目标函数"""
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
obj1 = self.objective_1(rho_phys)
obj2 = self.objective_2(rho_phys)
# 加权组合（最小化问题）
combined = -(w1 * obj1 + w2 * obj2)
return combined, np.mean(rho_phys)
def optimize(self, weights_list, n_iterations=25):
"""对不同权重组合进行优化"""
results = []
for w1, w2 in weights_list:
print(f" Optimizing with w1={w1:.2f}, w2={w2:.2f}")
rho = self.initialize()
for iter in range(n_iterations):
# 数值计算灵敏度
dc = np.zeros_like(rho)
eps = 1e-6
obj0, _ = self.combined_objective(rho, w1, w2)
for i in range(0, self.nx, 4):
for j in range(0, self.ny, 4):
rho_perturbed = rho.copy()
rho_perturbed[i, j] = min(rho_perturbed[i, j] + eps, 1.0)
obj_perturbed, _ = self.combined_objective(rho_perturbed, w1, w2)
dc[i, j] = (obj_perturbed - obj0) / eps
from scipy.ndimage import zoom
if dc.shape[0] < self.nx:
dc = zoom(dc, (self.nx/dc.shape[0], self.ny/dc.shape[1]), order=1)
dv = np.ones_like(rho) / (self.nx * self.ny)
rho_new = self.oc.update(rho, dc, dv)
if iter > 0 and iter % 8 == 0:
self.projection.beta = min(self.projection.beta * 1.2, 24.0)
rho = rho_new
# 最终评估
rho_filtered = self.filter.filter_density(rho)
rho_phys = self.projection.project(rho_filtered)
obj1 = self.objective_1(rho_phys)
obj2 = self.objective_2(rho_phys)
results.append({
'rho': rho,
'rho_phys': rho_phys,
'w1': w1,
'w2': w2,
'obj1': obj1,
'obj2': obj2
})
return results
# 运行多目标优化
print(" 运行多目标拓扑优化...")
multi_opt = MultiObjectiveOptimizer(nx=40, ny=40, volfrac=0.35)
weights_list = [(0.9, 0.1), (0.7, 0.3), (0.5, 0.5), (0.3, 0.7), (0.1, 0.9)]
results = multi_opt.optimize(weights_list, n_iterations=20)
# 绘制多目标优化结果
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
# Pareto前沿
obj1_values = [r['obj1'] for r in results]
obj2_values = [r['obj2'] for r in results]
axes[0, 0].scatter(obj1_values, obj2_values, s=200, c='red', zorder=5)
for i, (w1, w2) in enumerate(weights_list):
axes[0, 0].annotate(f'({w1:.1f},{w2:.1f})', (obj1_values[i], obj2_values[i]),
xytext=(5, 5), textcoords='offset points', fontsize=9)
axes[0, 0].plot(obj1_values, obj2_values, 'b--', alpha=0.5)
axes[0, 0].set_xlabel('Connectivity', fontsize=11)
axes[0, 0].set_ylabel('Smoothness', fontsize=11)
axes[0, 0].set_title('Pareto Front', fontsize=12, fontweight='bold')
axes[0, 0].grid(True, alpha=0.3)
# 不同权重下的设计
for idx, (i, j) in enumerate([(0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]):
if idx < len(results):
r = results[idx]
im = axes[i, j].imshow(r['rho_phys'].T, origin='lower', cmap='gray_r', vmin=0, vmax=1)
axes[i, j].set_title(f'w1={r["w1"]:.1f}, w2={r["w2"]:.1f}\nC={r["obj1"]:.3f}, S={r["obj2"]:.3f}',
fontsize=10, fontweight='bold')
axes[i, j].set_xlabel('x')
axes[i, j].set_ylabel('y')
plt.colorbar(im, ax=axes[i, j])
plt.tight_layout()
plt.savefig(f'{output_dir}/multiobjective_optimization.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 多目标拓扑优化图已保存")
# =============================================================================
# 第八部分：优化结果对比与总结
# =============================================================================
print("\n[8] 优化结果对比与总结...")
def create_summary_plot():
"""创建优化方法对比总结图"""
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
# 1. SIMP惩罚因子影响
p_values = [1, 2, 3, 4, 5]
grayness = []
for p in p_values:
# 模拟中间密度比例
rho_test = np.random.rand(1000)
epsilon = 1 + 9 * rho_test**p
# 计算"灰度"程度
gray = np.sum((rho_test > 0.1) & (rho_test < 0.9)) / len(rho_test)
grayness.append(gray)
axes[0, 0].plot(p_values, grayness, 'bo-', linewidth=2, markersize=8)
axes[0, 0].set_xlabel('Penalization Factor p', fontsize=11)
axes[0, 0].set_ylabel('Intermediate Density Ratio', fontsize=11)
axes[0, 0].set_title('Effect of Penalization Factor', fontsize=12, fontweight='bold')
axes[0, 0].grid(True, alpha=0.3)
# 2. 过滤半径影响
rmin_values = [1, 2, 3, 4, 5]
min_feature_size = [r * 2 for r in rmin_values] # 最小特征尺寸
axes[0, 1].plot(rmin_values, min_feature_size, 'rs-', linewidth=2, markersize=8)
axes[0, 1].set_xlabel('Filter Radius rmin', fontsize=11)
axes[0, 1].set_ylabel('Minimum Feature Size (cells)', fontsize=11)
axes[0, 1].set_title('Filter Radius vs Feature Size', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3)
# 3. 体积分数vs性能
volfracs = np.linspace(0.1, 0.6, 10)
# 模拟性能曲线
performance = 1 - np.exp(-5 * volfracs) # 饱和曲线
axes[0, 2].plot(volfracs, performance, 'g^-', linewidth=2, markersize=8)
axes[0, 2].set_xlabel('Volume Fraction', fontsize=11)
axes[0, 2].set_ylabel('Normalized Performance', fontsize=11)
axes[0, 2].set_title('Performance vs Material Usage', fontsize=12, fontweight='bold')
axes[0, 2].grid(True, alpha=0.3)
# 4. 收敛速度对比
iterations = np.arange(0, 50)
# OC方法
oc_conv = 1.0 * np.exp(-iterations / 10) + 0.05
# MMA方法（模拟）
mma_conv = 1.0 * np.exp(-iterations / 8) + 0.03
axes[1, 0].semilogy(iterations, oc_conv, 'b-', linewidth=2, label='OC Method')
axes[1, 0].semilogy(iterations, mma_conv, 'r--', linewidth=2, label='MMA Method')
axes[1, 0].set_xlabel('Iteration', fontsize=11)
axes[1, 0].set_ylabel('Objective Value (log)', fontsize=11)
axes[1, 0].set_title('Convergence Comparison', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
# 5. 计算成本分析
grid_sizes = [20, 40, 60, 80, 100]
# 计算成本大致与网格数成正比
compute_cost = [n**2 / 1000 for n in grid_sizes]
axes[1, 1].plot(grid_sizes, compute_cost, 'mv-', linewidth=2, markersize=8)
axes[1, 1].set_xlabel('Grid Size (NxN)', fontsize=11)
axes[1, 1].set_ylabel('Relative Compute Cost', fontsize=11)
axes[1, 1].set_title('Computational Cost Scaling', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)
# 6. 优化算法特性对比
methods = ['OC', 'MMA', 'SLP', 'SQP']
convergence = [8, 9, 7, 9] # 收敛速度评分
robustness = [9, 8, 7, 8] # 鲁棒性评分
complexity = [3, 7, 5, 8] # 实现复杂度
x = np.arange(len(methods))
width = 0.25
axes[1, 2].bar(x - width, convergence, width, label='Convergence', color='skyblue')
axes[1, 2].bar(x, robustness, width, label='Robustness', color='lightgreen')
axes[1, 2].bar(x + width, complexity, width, label='Complexity', color='lightcoral')
axes[1, 2].set_xlabel('Method', fontsize=11)
axes[1, 2].set_ylabel('Score (1-10)', fontsize=11)
axes[1, 2].set_title('Algorithm Comparison', fontsize=12, fontweight='bold')
axes[1, 2].set_xticks(x)
axes[1, 2].set_xticklabels(methods)
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig(f'{output_dir}/optimization_summary.png', dpi=150, bbox_inches='tight')
plt.close()
create_summary_plot()
print("✓ 优化结果对比总结图已保存")
# =============================================================================
# 第九部分：总结与输出
# =============================================================================
print("\n" + "=" * 60)
print("拓扑优化仿真完成！")
print("=" * 60)
print("\n生成的图表：")
print(" 1. simp_model.png - SIMP材料插值模型")
print(" 2. density_filter_projection.png - 密度过滤与投影")
print(" 3. adjoint_sensitivity.png - 伴随法灵敏度分析")
print(" 4. waveguide_optimization.png - 波导拓扑优化")
print(" 5. antenna_optimization.png - 天线结构优化")
print(" 6. multiobjective_optimization.png - 多目标拓扑优化")
print(" 7. optimization_summary.png - 优化结果对比总结")
print("\n所有图表已保存到:", output_dir)
print("=" * 60)
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
4
0- 0
已为社区贡献99条内容
所有评论(0)