はじめに
本記事では、Ansys Lumerical の varFDTDと逆設計LumOptを用いて、
クロス光導波路（waveguide crossing）を最適化する手法をご紹介します。
詳細は以下
https://optics.ansys.com/hc/en-us/articles/360042305314-Inverse-design-of-waveguide-crossing
従来のパラメータ掃引とは異なり、勾配ベース最適化を用いることで、
透過損失やクロストークを同時に低減する形状を効率的に探索できます。
まずは 2.5D（varFDTD） を用いて計算コストを抑えつつ設計空間を探索し、
得られた形状を必要に応じて 3D FDTD に展開する流れになります。
この検討では下の図に示すように上下左右にベースとなる導波路を設置して、クロス導波路を組んだ際にoutputが最も大きくなる様にクロス導波路の形状最適化を実施しています。
最適化された形状がこちら
今回の検討内容
上の図を見るとクロス導波路の根本部分に切れ目のようなものができていることや、構造自体が揺らいだ形状になっていることが確認できます。
こういった形状を見ると生成された形状は素子の作製プロセスに発生する形状のずれに対してロバスト性を持つのかが気になります。
そこで今回は以下のようなFOMを定義して、最適化ルーチンを回してみました。
図に示すように上下左右の光導波路の初期幅に対して0 nm(変動なし), ± 40 nm(変動アリ)とした場合について、それぞれのFOMを求めてその平均値が最大になるようにルーチンを組みなおしてみます。
実際のsample作成で40 nmもズレることはあり得ませんが、今回はテストということで。
スクリプト全体
lsfファイルのスクリプト(varFDTD_crossing_r1.lsf)
#############################################################################
# Scriptfile: varFDTD_crossing.lsf
#
# Description:
# Base simulation for adjoint shape-based optimization of SOI waveguide crossing in 2D
#
# Added:
# - Fabrication bias parameter delta_w (half-width bias)
# Access waveguide width is set by: wg_w = wg_w0 + 2*delta_w
##############################################################################
switchtolayout;
selectall;
delete;
## SIM PARAMS
size_x=5e-6;
size_y=5e-6;
mesh_x=20e-9;
mesh_y=20e-9;
finer_mesh_size=4.24e-6;
mesh_accuracy=4;
lam_c = 1.550e-6;
## MATERIAL
opt_material=addmaterial('Dielectric');
setmaterial(opt_material,'name','Si: non-dispersive');
n_opt = getindex('Si (Silicon) - Palik',c/lam_c);
setmaterial('Si: non-dispersive','Refractive Index',n_opt);
sub_material=addmaterial('Dielectric');
setmaterial(sub_material,'name','SiO2: non-dispersive');
n_sub = getindex('SiO2 (Glass) - Palik',c/lam_c);
setmaterial('SiO2: non-dispersive','Refractive Index',n_sub);
setmaterial('SiO2: non-dispersive',"color",[0, 0, 0, 0]);
## FAB VARIATION (added)
# delta_w is HALF-WIDTH bias:
# delta_w > 0 => wider
# delta_w < 0 => narrower
wg_w0 = 0.5e-6; # nominal full width
delta_w = 0; # <--
wg_w = wg_w0 + 2*delta_w;
## GEOMETRY
# INPUT WAVEGUIDE WEST (horizontal)
addrect;
set('name','input wg west');
set('x span',3.0e-6);
set('y span',wg_w);
set('z span',220e-9);
set('x',-3.5e-6);
set('y',0);
set('z',0);
set('material','Si: non-dispersive');
# OUTPUT WAVEGUIDE EAST (horizontal)
addrect;
set('name','output wg east');
set('x span',3.0e-6);
set('y span',wg_w);
set('z span',220e-9);
set('x',3.5e-6);
set('y',0);
set('z',0);
set('material','Si: non-dispersive');
# INPUT WAVEGUIDE SOUTH (vertical)
addrect;
set('name','input wg south');
set('x span',wg_w);
set('y span',3.0e-6);
set('z span',220e-9);
set('x',0);
set('y',3.5e-6);
set('z',0);
set('material','Si: non-dispersive');
# OUTPUT WAVEGUIDE NORTH (vertical)
addrect;
set('name','input wg north');
set('x span',wg_w);
set('y span',3.0e-6);
set('z span',220e-9);
set('x',0);
set('y',-3.5e-6);
set('z',0);
set('material','Si: non-dispersive');
# SUBSTRATE
addrect;
set('name','sub');
set('x span',10e-6);
set('y span',10e-6);
set('z span',4e-6);
set('y',0);
set('x',0);
set('z',0);
set('material','SiO2: non-dispersive');
set('override mesh order from material database',1);
set('mesh order',3);
set('alpha',0.8);
## varFDTD
addvarfdtd;
set('mesh accuracy',mesh_accuracy);
set('x min',-size_x/2);
set('x max',size_x/2);
set('y min',-size_y/2);
set('y max',size_y/2);
set('force symmetric y mesh',1);
set('y min bc','Anti-Symmetric');
set('z',0);
set('effective index method','variational');
set('can optimize mesh algorithm for extruded structures',1);
set('clamp values to physical material properties',1);
set('x0',-2.4e-6);
set('number of test points',3);
set('test points',[0,0; 0.0, 2.4; 2.4, 0]*1e-6);
## SOURCE
addmodesource;
set('direction','Forward');
set('injection axis','x-axis');
set('y',0);
set("y span",size_y);
set('x',-2.1e-6);
set('center wavelength',lam_c);
set('wavelength span',0);
set('mode selection','fundamental mode');
## MESH IN OPTIMIZABLE REGION
addmesh;
set('x',0);
set('x span',finer_mesh_size);
set('y',0);
set('y span',finer_mesh_size);
set('dx',mesh_x);
set('dy',mesh_y);
## OPTIMIZATION FIELDS MONITOR IN OPTIMIZABLE REGION
adddftmonitor;
set('name','opt_fields');
set('monitor type','2D Z-normal');
set('x',0);
set('x span',finer_mesh_size);
set('y',0);
set('y span',finer_mesh_size);
set('z',0);
## FOM FIELDS
adddftmonitor;
set('name','fom');
set('monitor type','Linear Y');
set('x',2.1e-6);
set('y',0);
set('y span',size_y);
set('z',0);
#save;
Pythonスクリプト(crossing_opt_varFDTD_mean_robust.py)
# -*- coding: utf-8 -*-
"""
Mean-robust optimization for SOI waveguide crossing (2D varFDTD + LumOpt)
- Fabrication variation modeled as uniform width bias: delta_w = {-20nm, 0, +20nm}
- Objective: maximize expected FOM (average over cases)
Output:
- Save nominal (delta_w=0) MODE project file (.lms) with the optimized polygon placed.
"""
import os, sys, re, time
import numpy as np
import scipy as sp
from scipy.optimize import minimize
# ====== 0) Project path ======
PROJECT_DIR = os.path.dirname(os.path.abspath(__file__))
os.chdir(PROJECT_DIR)
# ====== 1) Lumerical API path ======
LUM_ROOT = r"C:\Program Files\Lumerical\v252"
LUM_API_PY = os.path.join(LUM_ROOT, "api", "python")
try:
os.add_dll_directory(LUM_ROOT)
except Exception:
pass
if LUM_API_PY not in sys.path:
sys.path.insert(0, LUM_API_PY)
import lumapi
# ====== 2) LumOpt import ======
from lumopt.utilities.wavelengths import Wavelengths
from lumopt.geometries.polygon import FunctionDefinedPolygon
from lumopt.utilities.materials import Material
from lumopt.figures_of_merit.modematch import ModeMatch
from lumopt.optimizers.generic_optimizers import ScipyOptimizers
from lumopt.optimization import Optimization
# ====== 3) Robust LSF loader (avoid cp932 decode issues) ======
def load_lsf_text(path: str) -> str:
"""
Read LSF robustly on Windows:
try utf-8 then cp932, fallback to utf-8 with replacement.
Remove '#' comments similar to load_from_lsf.
"""
lines = None
for enc in ("utf-8", "cp932"):
try:
with open(path, "r", encoding=enc) as f:
lines = f.readlines()
break
except UnicodeDecodeError:
pass
if lines is None:
with open(path, "r", encoding="utf-8", errors="replace") as f:
lines = f.readlines()
cleaned = [line.strip().split("#", 1)[0] for line in lines]
return "\n".join([l for l in cleaned if l != ""])
# ====== 4) Helpers: inject delta_w, and tag all LSF names (optional) ======
def base_script_with_delta(base_script: str, delta_w: float) -> str:
"""
Replace the exact token 'delta_w = 0;' in the LSF with 'delta_w = <value>;'
(only the first occurrence is replaced)
"""
key = "delta_w = 0;"
rep = f"delta_w = {delta_w};"
if key not in base_script:
raise RuntimeError("LSF does not contain the exact token 'delta_w = 0;'. Please verify varFDTD_crossing_r1.lsf.")
return base_script.replace(key, rep, 1)
def add_cell_tag_to_names(lsf_text: str, tag: str) -> str:
"""
Append '_<tag>' to all object names in LSF:
set('name','input wg west') -> set('name','input wg west_<tag>')
This is mainly for GDS cell-name uniqueness. (Not required for saving .lms)
"""
pattern = r"set\('name','([^']+)'\)"
def repl(m: re.Match) -> str:
name = m.group(1)
return f"set('name','{name}_{tag}')"
return re.sub(pattern, repl, lsf_text)
# ====== 5) Read LSF ======
LSF_PATH = os.path.join(PROJECT_DIR, "varFDTD_crossing_r1.lsf")
crossing_base = load_lsf_text(LSF_PATH)
# ====== 6) Wavelengths ======
wavelengths = Wavelengths(start=1300e-9, stop=1800e-9, points=25)
# ====== 7) Fabrication variation cases (half-width bias) ======
DELTA_WS = np.array([-20e-9, 0.0, +20e-9], dtype=float)
WEIGHTS = np.array([1/3, 1/3, 1/3], dtype=float)
# ====== 8) Geometry function with delta_w ======
def cross(params, delta_w=0.0):
y_end = params[-1]
x_end = 0 - y_end
points_x = np.concatenate(([-2.01e-6], np.linspace(-2e-6, x_end, 10)))
points_y = np.concatenate(([0.25e-6], params))
# fabrication bias on half-width
points_y = points_y + delta_w
# keep feasible
points_y = np.clip(points_y, 0.05e-6, 1.0e-6)
n_interpolation_points = 50
polygon_points_x = np.linspace(min(points_x), max(points_x), n_interpolation_points)
interpolator = sp.interpolate.interp1d(points_x, points_y, kind='cubic')
polygon_points_y = interpolator(polygon_points_x)
polygon_points_y = np.clip(polygon_points_y, -1e-6, 1e-6)
pplu = [(x, y) for x, y in zip(polygon_points_x, polygon_points_y)]
ppld = [(x, -y) for x, y in zip(polygon_points_x, polygon_points_y)]
ppdl = [(-y, x) for x, y in zip(polygon_points_x, polygon_points_y)]
ppdr = [(y, x) for x, y in zip(polygon_points_x, polygon_points_y)]
pprd = [(-x, -y) for x, y in zip(polygon_points_x, polygon_points_y)]
ppru = [(-x, y) for x, y in zip(polygon_points_x, polygon_points_y)]
ppur = [(y, -x) for x, y in zip(polygon_points_x, polygon_points_y)]
ppul = [(-y, -x) for x, y in zip(polygon_points_x, polygon_points_y)]
polygon_points = np.array(
pplu[::-1] + ppld[:-1] + ppdl[::-1] + ppdr[:-1] +
pprd[::-1] + ppru[:-1] + ppur[::-1] + ppul[:-1]
)
return polygon_points
# ====== 9) Design variables / bounds ======
initial_params = np.linspace(0.25e-6, 0.6e-6, 10)
bounds = [(0.2e-6, 1e-6)] * initial_params.size
# ====== 10) Materials / thickness ======
eps_in = Material(name='Si: non-dispersive', mesh_order=2)
eps_out = Material(name='SiO2: non-dispersive', mesh_order=3)
depth = 220e-9
# ====== 11) FOM ======
mode_fom = ModeMatch(
monitor_name='fom',
mode_number='fundamental mode',
direction='Forward',
multi_freq_src=False,
target_T_fwd=lambda wl: np.ones(wl.size),
norm_p=1
)
# ====== 12) Create ONE Optimization per delta_w and initialize sim ======
dummy_optimizer = ScipyOptimizers(
max_iter=1,
method='L-BFGS-B',
scaling_factor=1e6,
pgtol=1e-5,
ftol=1e-5,
scale_initial_gradient_to=0.0
)
opts = []
for dw in DELTA_WS:
poly = FunctionDefinedPolygon(
func=(lambda p, dw=float(dw): cross(p, delta_w=dw)),
initial_params=initial_params,
bounds=bounds,
z=0.0,
depth=depth,
eps_out=eps_out,
eps_in=eps_in,
edge_precision=5,
dx=1e-9
)
opt_k = Optimization(
base_script=base_script_with_delta(crossing_base, float(dw)),
wavelengths=wavelengths,
fom=mode_fom,
geometry=poly,
optimizer=dummy_optimizer,
use_var_fdtd=True,
hide_fdtd_cad=False,
use_deps=True,
plot_history=False,
store_all_simulations=False
)
# Initialization trigger (your LumOpt requires dx)
opt_k.check_gradient(initial_params, dx=1e-3)
opts.append(opt_k)
def eval_one_case(opt_k: Optimization, params: np.ndarray):
"""
Evaluate single-case FOM and gradient for given params using initialized opt_k.
Your LumOpt: calculate_gradients() takes no args.
"""
fom_val = opt_k.callable_fom(params)
opt_k.calculate_gradients()
if hasattr(opt_k, "last_grad") and opt_k.last_grad is not None:
grad = np.array(opt_k.last_grad, dtype=float)
elif hasattr(opt_k, "gradients") and opt_k.gradients is not None:
grad = np.array(opt_k.gradients, dtype=float)
else:
raise RuntimeError("Gradient not available (opt.last_grad / opt.gradients).")
grad = grad.reshape(params.shape)
return float(fom_val), grad
def eval_mean_objective(params: np.ndarray):
"""
Expected objective: mean over fabrication cases.
Returns (J, gradJ) for maximization.
"""
J = 0.0
g = np.zeros_like(params, dtype=float)
for opt_k, w in zip(opts, WEIGHTS):
fk, gk = eval_one_case(opt_k, params)
J += float(w) * fk
g += float(w) * gk
return J, g
# ====== 13) Run external SciPy optimizer (maximize J -> minimize -J) ======
x0 = np.array(initial_params, dtype=float)
def fun_and_jac(x):
x = np.array(x, dtype=float)
J, g = eval_mean_objective(x)
return -J, -g
res = minimize(
fun=lambda x: fun_and_jac(x),
x0=x0,
method="L-BFGS-B",
jac=True,
bounds=bounds,
options={"maxiter": 35, "ftol": 1e-5, "gtol": 1e-5}
)
best_params = np.array(res.x, dtype=float)
np.savetxt(os.path.join(PROJECT_DIR, "2D_parameters_mean_robust.txt"), best_params)
# ====== 14) Save MODE project file (nominal: delta_w = 0) ======
EXPORT_DELTA_W = 0.0 # nominal mask geometry
# もし毎回ユニーク名にしたいなら timestamp を付ける
TAG = time.strftime("%Y%m%d_%H%M%S")
project_name = f"crossing_mean_robust_nominal_{TAG}.lms"
project_path = os.path.join(PROJECT_DIR, project_name)
# ※ .lms 保存には cell 名のユニーク化は基本不要（GDS向けの処理）
# ただし将来GDS化する可能性があるなら付けてもOK
# lsf_for_save = add_cell_tag_to_names(base_script_with_delta(crossing_base, EXPORT_DELTA_W), f"save_{TAG}")
lsf_for_save = base_script_with_delta(crossing_base, EXPORT_DELTA_W)
with lumapi.MODE(hide=False) as mode:
mode.cd(PROJECT_DIR)
mode.eval("switchtolayout; selectall; delete;")
mode.eval(lsf_for_save)
# 最適化形状（名目：delta_w=0）
mode.addpoly(vertices=cross(best_params, delta_w=EXPORT_DELTA_W))
mode.set('name', f'optimized_poly_{TAG}')
mode.set('x', 0.0)
mode.set('y', 0.0)
mode.set('z', 0.0)
mode.set('z span', depth)
mode.set('material', 'Si: non-dispersive')
# プロジェクト保存
mode.eval(f"save('{project_name}');")
print("Optimization finished.")
print("Saved params: 2D_parameters_mean_robust.txt")
print("Saved MODE project:", project_name)
print("SciPy status:", res.message)
このスクリプト群は、SOIクロス光導波路を対象に、varFDTD（2D）とLumOptの随伴法を用いたロバスト最適化を実装しています。LSFファイルでは、Si/SiO₂材料、4方向の入出力導波路、varFDTDソルバ、光源およびFOM用モニタを定義し、交差構造の電磁界解析の基盤を構築しています。導波路幅の加工ばらつきは、半幅バイアス delta_w として導入され、太る・細る誤差を連続的に表現しています。
Python側では、この delta_w を −20 nm、0、+20 nm の3ケースに展開し、それぞれについて独立にFOMと勾配を計算します。交差部形状は少数の設計変数からポリゴンとして生成され、対称性を保ったまま最適化されます。各加工ばらつきケースのFOMを等確率で平均化することで、名目条件だけでなく製造誤差に強い設計を目的関数として定義しています。最適化自体はSciPyのL-BFGS-Bを用いて外部から制御され、LumOptは随伴法による高速な勾配評価に専念します。最終的に得られた名目形状はMODEプロジェクト（.lms）として保存され、後段の3D検証やGDS化に利用できる形で出力されます。
出力結果
こちらが出力結果になります。
右側の従来と比べて、構造の揺らぎがつぶれて消えていることが確認できます。
ちょっとテーパがついてしまってますが。。40 nmも導波路幅を変えたのはやりすぎだったかも。
次回以降もう少しパラメータを調整して報告します。
※本記事は筆者個人の見解であり、所属組織の公式見解を示すものではありません。
問い合わせ
光学シミュレーションソフトの導入や技術相談、
設計解析委託をお考えの方はサイバネットシステムにお問合せください。
光学ソリューションサイトについては以下の公式サイトを参照：
👉 光学ソリューションサイト(サイバネット)
光学分野のエンジニアリングサービスについては以下の公式サイトを参照：
👉 光学エンジニアリングサービス(サイバネット)