計算数理工学論文集 Vol. 20 (2020 年 12 月), 論文 No. 08-201219

JASCOME

共分散行列適応進化戦略と S 行列を用いた
フォノニック構造のパラメータ最適化
Parameter optimisation for phononic structure using CMA-ES and S-matrix
飯盛浩司 1) ，松島慶 2) ，鶴田健二 3) , 高橋

徹 4) , 松本敏郎 5)Hiroshi ISAKARI, Kei MATSUSHIMA, Kenji TSURUTA, Toru TAKAHASHI, Toshiro MATSUMOTO
1) 名古屋大学大学院工学研究科(〒 464-8603

名古屋市千種区不老町,

E-mail: isa@nagoya-u.jp)

2) 名古屋大学大学院工学研究科(〒 464-8603

名古屋市千種区不老町,

E-mail: k matusima@nuem.nagoya-u.ac.jp)

3) 岡山大学大学院自然科学研究科(〒 700-8530

岡山市北区津島中 3-1-1,

E-mail: tsuruta@okayama-u.ac.jp)

4) 名古屋大学大学院工学研究科(〒 464-8603

名古屋市千種区不老町,

E-mail: toru.takahashi@mae.nagoya-u.ac.jp)

5) 名古屋大学大学院工学研究科(〒 464-8603

名古屋市千種区不老町,

E-mail: t.matsumoto@nuem.nagoya-u.ac.jp)

This paper presents a parameter optimisation strategy for phononic structures. The
phononic structures here are assumed consisting of multiple scatterers. The rotation angles of each scatterer are used as the design parameters to optimise an objective function
related to the performance of the structure. To explore an optimal design, we here use
a metaheuristic approach called CMA-ES in which the S-matrix method is employed toevaluate the objective function eﬃciently. After checking the performance of the proposed method with a simple optimisation problem, we demonstrate that it can provide
an optimal design of a phononic waveguide.
Key Words : CMA-ES, S-matrix, Parameter Optimisation, Phononic structure
づらい幾何学形状を導く. 実際, 複雑な構造を造形できると

1. 諸言
近 年, 音 波 の 波 長 程 度 の 構 造 周 期 を 持 つ フォノ ニック 構

される積層造形を用いたとしても, 素朴なトポロジー最適化

造 (1) , 波長よりも短い周期をもつ音響メタ構造 (2) など、自による設計案はそのまま製造できないことが指摘されてい然界の材料にはない音響特性を示す人工構造体の研究開発とる (10, 11) . また, 物理的な考察に基づき大まかには形状が定応用が進められている.

まっているデバイスの形状を最適化する目的には使いづらい

フォノニック構造の設計開発においてトポロジー最適化をはじめとする構造最適化が有用であると考えられており, 特という側面もある.
一方で, フォノニック構造は必ずしも複雑な形状である必に, 広いバンドギャップを持つフォノニック構造のトポロジー

要はなく, 比較的単純な構造を周期的に配置することでバン

最適化 (3, 4) に関する研究が盛んに行われている. また, フォ

ドギャップを実現できることが知られている (1) . したがって,

ノニック導波路のトポロジー最適化 (5) に関する研究も見ら

製造の容易な単位構造を, 所望の機能を持つように配置する

れる. 著者らの研究グループにおいても, 波動問題の解析に方法論を構築することには意義がある (12) . このような最適(6, 7)に基づくフォ

化問題は寸法最適化に分類され, 設計変数に対する目的関数の開発を行い, 従来のの感度が評価できれば比較的簡単に勾配法に基づく最適化を

有限要素法に基づく方法と比較して高効率な設計が可能であ構築することができる. しかしながら, 設計変数の数が多いることを示してきた. その他, フォノニック構造のトポロジー

場合や目的関数が多峰性の場合には, 初期値やステップ幅な適した高速境界要素法とレベルセット法
ノニック構造のトポロジー最適化(8)

最適化に関する研究についてはレビュー論文(9)

どの最適化に関するパラメータを適切に選択することは難しに詳しい.

トポロジー最適化はその設計自由度の高さ故, 従来のデバ

い. そのような場合には, 所謂メタヒューリスティックな方法

イスとは全く異なるデバイス形状を創成し, 抜本的に性能をの利用もひとつの選択肢である. 中でも, 共分散行列適応進改善する, あるいはこれまでにない機能を持つデバイスを実化戦略 (CMA-ES)(13, 14) は最適化にあたってユーザが設定

現することを可能にする場合がある一方で, しばしば製造しすべきパラメータの数が極めて少ないという点において魅力的である. 開発当初から形状最適化問題への適用例が見ら2020 年 10 月 20 日受付，2020 年 11 月 14 日受理− 65 −

境界値問題に支配される.

れる (15, 16) 他, 近年ではトポロジー最適化 (17, 18, 19) への適用に関する研究も進められている. 特に, 藤井ら(18, 19)の一∇2 u(x) + k2 u(x) = 0

連の研究は顕著であり, 主にクローキングデバイスのトポロ

x∈Ω

(1)

本研究では, 多数の散乱体からなるフォノニック構造の最∂u(x)=0 x∈Γ∂n( sc)∂u (x)
lim− ikusc (x) = 0
|x|→∞
∂|x|

適配置決定問題を対象とし, 各散乱体の回転角を設計変数とここに, k := ω/c は波数, ω は角周波数, c は波速, p は入射波する最適設計問題を CMA-ES により解く方法を提案する.

の進む方向を表す単位ベクトル, n は Ω の外向き単位法線ベ

CMA-ES は多点探索に基づく方法であるから, 勾配法と比べ

クトル, usc := u − uin は散乱波である.

て目的関数の評価を多数回行う必要があり, したがって目的

2.2. 最適化問題の定義

ジー最適化を対象としつつ, 不等式制約の取り扱いなどに関する CMA-ES そのものの理論の改良を行っている.

関数の評価を効率的に行うことが鍵となる. 本研究では, 境界要素法と S 行列法を組み合わせた方法によりこれを実現∅ であると仮定する.
ci を中心として Ωsi を反時計回りに ϕi 回転させて得られるが, 近年, 境界要素法との組み合わせが提案されたことによりその適用可能範囲が拡大した(3)s以降, 各散乱体 Ωsi を包含する中心 ci の円 Bi を考え, ∩N
i=1 Bi =する. S 行列法は古くから知られる散乱問題の解法 (20) であ
(21)

(2)

s
s
る散乱体を Ωsi (ϕi ) と書く. また, Ω(ϕ) := R2 \ ∪N
i=1 Ωi (ϕi ) を. 当手法の利点の一つ

定義しておく. ここに, ϕ := (ϕ1 , · · · , ϕNs )t である. ここで,は, 多重散乱問題の解が与えられた際に, 個々の散乱体を回転以下の最適化問題を考える.

した場合の解を, 小さいサイズの行列に関する代数操作で構min J(ϕ)

成できることにある. したがって, 本研究では, CMA-ES と S

ϕ

行列法の組み合わせにより, 多数の散乱体の最適回転角を効(4)

ここに, J(ϕ) は Ω(ϕ) における多重散乱問題 (1)–(3) の解を率良く求める方法論を新規に構築することを目的とする.

用いて表される汎関数である.

以下, 2.1, 2.2 節において本研究で対象とする最適化問題

2.3. 共分散行列適応進化戦略 (CMA-ES)を定義し, 2.3 節で CMA-ES の概略を述べ, 2.4 節で目的関数

本論文では, 2.2 節で定義した最適化問題 (4) の解法としての評価に用いる S 行列法について解説する. その後, 3 節で共分散行列適応進化戦略 (CMA-ES)(13, 14) を用いる. CMA-

は, 提案手法の性能評価に関する数値計算結果と併せてフォES は多点探索に基づくメタヒューリスティックな最適化法でノニック導波路の最適設計案を示し, 提案手法の有効性を検ある. すなわち,

証する.

1. ある分布から乱択した多数の解候補に対する目的関数

2. 定式化

値を計算

2.1. 多重散乱問題

2. その結果を基に解候補の分布を更新

本 論 文 で は, Fig. 1 に 示 す よ う な Ns 個 の 剛 体 Ωsi (i =を繰り返すことで最適解を探索する. CMA-ES では, 解候補

1, · · · , Ns ) が配置された 2 次元領域における音波の散乱問を平均 m ∈ RNs , 共分散 σ 2 C ∈ RNs ×Ns (σ > 0) の多変量正題を考える. 平面音波 uin (x) = exp(ikp · x) が領域 Ω :=規分布から生成する. アルゴリズムの概略は以下のとおりである:

Step 0. 最適化のひとつのステップにおいて生成する解候補の数 (個体数) λ 及び m, σ, C の初期値を与える. 本研
究では, feasible set が li ≤ ϕi ≤ ui (i = 1, · · · , Ns ) と書けるとき, mi = li , σ = maxi (ui − li )/2, C = I =diag(1, · · · , 1) と与える. λ の設定については後述する.
Step 1. 標準正規分布 N (0, I) から λ 個の独立なサンプル z i ∈RNs を 生 成 し, ϕi = m + σy i = m + σBDz i (i =1 · · · λ) を解候補とする. ここに, B, D は, C の固有分
Fig. 1

解 C = BD2 Bt に現れる行列であり, したがって D は対Multiple scattering.

角行列である.

Step 2. 生成した全ての解候補 ϕi に対し, 目的関数 J(ϕi ) をss
R2 \ ∪N
i=1 Ωi に入射したとき, Ω における音場 u(x) は以下の

計 算 し, J(ϕi ) の 昇 順 に ϕi を 並 べ, 目 的 関 数 が j 番
目 に 小 さ い 解 候 補 を ϕj;λ と 書 く. こ れ に 対 応 し て,

y j;λ = (ϕj;λ − m)/σ も定義しておく.

− 66 −

Step 3. y j;λ の加重平均 < y >w =∑µi=j wj y j;λ を計算し, 以下

べる.

により m, σ, C を更新する.

m ← m + σ < y >wcσ||pσ ||σ ← σ exp  √ (dσ Ns 1 − 1 +4Ns

C ← C + c1 pc ptc + cµλ∑

以降, S 行列の成分を数値的に計算する方法について述
定式化に必要となる Graf の加法定理は以下で与えられる.


1
21Ns2

)

wj y j;λ y tj;λ

(5)
(6)

(7)

j=1

ここに, pσ , pc ∈ R は「進化パス」と呼ばれ, 零ベクト
ルを初期値として次式で更新する:√
1
pσ ← (1 − cσ )pσ + cσ (2 − cσ )µeﬀ C− 2 < y >w (8)√pc ← (1 − cc )pc + cc (2 − cc )µeﬀ < y >w
(9)
∑
1
2
ここに, C− 2 = BDBt , µeﬀ = 1/ µ
j=1 wj であり, 重み∑µwj は j=1 wj = 1 を満たすように選ぶ. また, µ < λ,c1 , cσ , dσ は定数である. なお, 本研究では λ 以外の全ての定数として, 原著 (14) に示された推奨値を用いる.

Hn(1) (k|x − y|)einarg(x−y)
 ∞
∑

(1)

Hm
(k|x|)Jn−m (k|y|)eimarg(x) e−i(m−n)arg(y)




m=−∞





if |x| > |y|
=
∞
∑

(1)


Hm
(k|x|)Jn−m (k|y|)eimarg(x) e−i(m−n)arg(y)




m=−∞




if |x| < |y|(13)
さて, 単一の散乱体に入射する適当な平面音波に対する散乱波は以下の積分表現を持つ.
∫)i∂ ( (1)
usc (x) = −H0 (k|x − y|) u(y)dΓ(y) (14)4 ∂Ωs ∂n(y)(14) に現れる Hankel 関数に対して Graf の加法定理 (13) 上

以上を繰り返すことで, m は最適解 ϕ∗ に, σ は零に収束する

式を用い, (11) と比較すれば
∫
)
i
∂ (
bm = −Jm (k|y − c|)e−imarg(y−c) u(y)dΓ(y)4 ∂Ωs ∂n(y)ことが知られている (14) . 本研究では, CMA-ES の実装として(15)

Step 4. 1. から 3. を繰り返す.

libcmaes(22) を用い, 負の重みを許容する active CMA-ES(23)を選択した.

となることが直ちに分かる. S 行列の成分 Smn は (15) の u
を以下の入射波2.4. S 行列を用いた多重散乱解析
Jn (k|x − c|)einarg(x−c)

2.3 節で述べた CMA-ES を用いて最適化問題 (4) を解く場
合, 多重散乱問題 (1)–(3) の解で構成される目的関数 J(ϕ) を多数回計算する必要があり, 通常の偏微分方程式ソルバを用いることは計算効率の観点から現実的でない. そこで本研究では, S 行列を利用した高速算法により目的関数を評価する.

S 行列の定義を記述するため, 単一の散乱体 Ωs を考える.
Ωs の補領域 R2 \ Ωs で定義される外部境界値問題 (1)–(3) の解 u = uin + usc は, Ωs を包含する円 B の補領域において, 以
下のように B の中心 c における円筒波展開の形で書ける.
M∑

uin (x) = lim

M →∞

M →∞

に対する散乱問題の解 un に置き換えれば得られる. 本研究では, un (n = −M, · · · , M ) の計算に境界要素法を用いる.
したがって, 2M + 1 回の境界要素解析を実行する必要がある
が, 境界積分方程式を離散化して得られる 2M + 1 本の代数方程式の係数行列 A は全て共通であるから, A を予め LU 分
解しておき, その結果を再利用することで効率的に解析可能である.
次に, S 行列を用いた多重散乱問題の解法について述べる.

am Jm (k|x − c|)eimarg(x−c)

(10)

(1)
b m Hm(k|x − c|)eimarg(x−c)

(11)

m=−M
M
∑

usc (x) = limΩsi の S 行列を S(i) , Ωsi の中心 ci における入射平面波の円筒
波展開の係数を並べたベクトルを a(i) , 放射する波のそれを

b(i) と書く. j 番散乱体から放射された波

m=−M

lim

こ こ に, arg(x) = tan−1x2 /x1 で あ り, Jm ,

(16)

(1)
Hm は 各々m

M →∞

M
∑

(j)

(j) (1)bn
Hn (k|x − c(j) |)einarg(x−c

)

(17)

n=−M

次の Bessel 関数, 第一種 Hankel 関数である. また, am =は, i 番散乱体に入射するが, (17) に現れる Hankel 関数に(p2 + ip1 )m exp(ikp · c), bm ∈ C は円筒波展開の係数であGraf の加法定理 (13) 下式を用いれば, その c(i) における円り, bm は未知である. 今, 問題の線形性から, {ai }M
i=−M を筒波展開の m 次の係数は{bi }M
i=−M に変換する作用素もまた線形であり, したがって以下のように行列 S で表現できる.

lim

M →∞

b = Sa(12)

ここで, a = (a−M , · · · , aM )t , b = (b−M , · · · , bM )t ∈ R2M +1を定義した. (12) に現れる S ∈ C(2M +1)×(2M +1)

を S (=scat-tering) 行列と言う.

− 67 −

M
∑

(j)
T(ij)
mn bn

(18)

n=−M

と書けることが分かる. ここに, T(ij) は T (translation) 行列
と呼ばれ,
(j)

(j)
T(ij)
− x(i) |)e−i(m−n)arg(xmn = Hm−n (k|x(1)−x(i) )

(19)

2. O(Neval (2M + 1)Ns ): S 行列の回転 (25)と書ける.
さて, Ωsi には, Ns − 1 個の散乱体 Ωsj (j ̸= i) から放射された波に加え, 入射波も入射する. したがって, Ωsi に入射する全ての波の円筒波展開の係数はa(i) +∑

T(ij) b(j)(20)

る. 通常, M, Ns ≪ Ne , Neval であるから, S 行列を用いた方法が効率的であることが分かる. なお, Ns 個全ての散乱体の形状が異なる場合においては, 各々の散乱体に対応する S 行となるが, S 行列の定義 (12) より, これは()∑ (ij) (j)
(i)
(i)
(i)
b =S
a +
T b

列は全て異なる. この時, precomputation にかかる時間は上

(21)

j̸=i

記の estimate のおよそ Ns 倍となる一方で, 上記 2. 及び 3. にかかる計算時間は全ての散乱体が合同である場合と比べてほを満たす. (21) を i = 1 から Ns について連立して解けば, 未知のbである. 一方, 通常の境界要素法を用いる場合には, 反復法
を用いて代数方程式を解いた場合でも O(Neval (Ns Ne )2 ) であ

j̸=i

(i)

3. O(Neval (2M + 1)2 Ns2 ): 代数方程式 (21) の求解s
なお, (21) を {b(i) }N
i=1(i = 1, · · · , Ns ) が全て求まる.

について整理して得られる代数方程式の係数行列は I − Uと書けるが, この行列は良条件であることが知られており,||U|| < 1 を満たす (24) . したがって本研究では, (21) の解法として前処理なしの GMRES を採用する. 実際の数値計算においては, S 行列のサイズ 2M + 1 を有限の値に設定することんど変化しない. ただし, 各散乱体の要素数及び S 行列の
サイズ M は異なり得ることに注意が必要である.

3. 計算結果
本節では, 本論文で提案する CMA-ES と S 行列に基づくパラメータ最適化法の validation と性能評価 (3.1 節) とフォノ
ニック導波路の最適設計 ( 3.2 節) に関する数値例を示す. いとで近似解法を構成する。ここでは、所謂 Rokhlin の式 (25)ずれの例においても, Ω は波速 c = 1.0 の媒質で満たされていM = kd + ntol log(kd + π)

タ ntol は 5 とし, 代数方程式 (21) の求解に用いる GMRES の(22)を基に M を定める. ここに, d = maxi (diam Bi ), ntol は要求する精度に応じて定める定数である.
次に, Ωs の S 行列 S が既知である時に, これを半時計周りに ϕ だけ回転して得られる Ωs (ϕ) のそれ S̃ を求めることを考える. Ωs 周辺の波動場が (10), (11) と書けている時, Ωs (ϕ)
周辺のそれはuin (x) = lim

M →∞

m=−M

,··· ,eらの距離 r が r = 0.8 + 0.4 cos t (0 ≤ t ≤ 2π) と表される 2−iM ϕ) S diag(e−iM ϕ

測点 xobs = (1.0, 0.0)t に収束させる問題を考える. すなわち,

1
J(ϕ) = − |u(xobs )|22(26)

ここに, u は R2 \ ∪2i=1 Ωs (ϕ) で定義された境界値問題 (1)–(3)と書けるから, S 行列の定義 (12) からS̃ = diag(eFig. 2 に示すような, 中心 c1 = (−1.5, 0)t , c2 = (1.5, 0)t か

以下の目的関数 J(ϕ) を最小化する最適化問題を考える.

(1)b m Hm(k|x − c|)eim(arg(x−c)−ϕ)

(24)

iM ϕ

はじめに, 提案法の validation を行うため, 設計変数の数
が 2 つである簡単な最適化問題を考える.

せ, ω = 6.0 で x1 正の向きに進む単位振幅の平面音波を, 観

m=−M

M
∑

3.1. validation と性能評価

各々の散乱体を ci (i = 1, 2) 中心に反時計回りに ϕi 回転さ

am Jm (k|x − c|)eim(arg(x−c)−ϕ)(23)

usc (x) = lim

許容誤差は 10−5 とした.

つの閉曲線を境界とする散乱体 ∪2i=1 Ωs (0) を初期形状とし,

M
∑

M →∞

るとした. また, S 行列の精度コントロール (22) のパラメーの解である. u の計算には, 単一の散乱体の表面を 200 分割し,··· ,e

iM ϕ

)

(25)

た境界要素を用いた. この時, 波長と最も長い境界要素の長さの比は 22 である. また, 形状の対称性から制約条件として,となる. なお, 実際には (25) を用いて陽に S̃ を計算せず, 代数方程式 (21) を GMRES で解く際に, 「右」から順にベクト
ルとの乗算を計算する.

0◦ ≤ ϕ1 , ϕ2 ≤ 120◦(27)を課す. なお, libcmaes においては box-type の実行可能領域

さて, CMA-ES による最適化の過程において, 目的関数の評価を Neval 回行うとする. また, 簡単のため, Ns 個の散乱を設定することができる. 実装の詳細についてはソースコー
ド (22) を参照されたい.

体は全て合同であり, 各散乱体は Ne 個の境界要素に分割さ実行可能領域 (27) において目的関数 J を sweep した結れているとする. この時, 上記の S 行列を用いた方法の計算

果を, Fig. 3 に −J の等高線として示す. いくつかの局所の内, 主要な部分の計算量は,解が存在することが確認でき, (26), (27) が多峰性の最適化
問 題 で あ る こ と が 分 か る. な お, global な 最 適 解 は ϕ∗ =1. precomputation• O(Ne3 ): 単一の散乱体に対する境界積分方程式の求解にかかる LU 分解(58.7, 84.0)t と ϕ∗ = (36.0, 61.3)t 付近に存在し, この時の目的関数値は J(ϕ∗ ) = −1.19 であった. また, global な最適解を除いた局所最適解の中で最も目的関数値が小さいものは• O((2M + 1)Ne2 ): S 行列の生成ϕ∗∗ = (11.1, 109) であった (J(ϕ∗∗ ) = −1.06).

− 68 −

定した場合, ほとんど全ての試行において少なくとも ϕ∗∗ に3は収束することが確認できる. また, 参考として Table 2 に

2

S 行列を用いない素朴な境界要素法を用いた場合の同じ表を示す. λ = 6 の場合に ϕ∗ に収束する回数がやや増加するこ

1と, λ = 30 の場合に ϕ∗ に収束する回数がやや減少すること
0などの差異はあるものの, 概ね S 行列を用いた場合と同様の傾向を示すことが分かる. 以上の結果より, 先行研究 (14, 19)-1でも指摘されているように, 個体数 λ を推奨値よりも大きく取ることで大域的最適解に近い目的関数値を取る設計変数の-2

組み合わせが見つかる可能性は高くなると期待される. ただ
-3
-3

-2

-1

0

1

2

し, 一ステップあたりの計算時間は個体数に比例することに3

注意する必要がある. また, 適切な λ の設定値は問題依存でFig. 2あると予想される.

Settings of the validational optimisation
problem in which the angles of two star-

Table 1

shaped rigid scatterers are optimised to

The number of successful trials with Smatrix method.

focus a sound wave on xobs .

120

1.00
0.75
0.50
0.25

80

ϕ∗ に収束ϕ∗ または ϕ∗∗ に収束

6 (推奨値)

55

80

18

71

99

30

78

99

60

77

100

Table 2400λThe number of successful trials withBEM.

0Fig. 3λϕ∗ に収束ϕ∗ または ϕ∗∗ に収束

6 (推奨値)

67

84

18

78

98

Contour map of −J in (26) within the fea-

30

67

96

sible set (27) for the validational optimi-

60

83

100

40

80

120

sation.
Fig. 4, Fig. 5 に, ϕ∗ = (58.7, 84.0)t (大域的最適解のひと
以下, 提案法を用いて最適化問題 (26), (27) を解いた結果つ), ϕ∗∗ に対する散乱体周辺の音場強度を示す. ϕ∗ , ϕ∗∗ のを示す. CMA-ES には全ての戦略パラメータに推奨値が与いずれの場合にも, 想定した観測点付近に音場が集中してい

えられているが, 個体数 λ に関しては多峰性の最適化問題にる様子が見て取れる. ϕ∗∗ は上下対称な形状となっている一対しては推奨値 4 + [3 ln N s ] よりも大きな値をとる必要があ方で, 大域的な最適解 ϕ∗ に対応する設計案は上下非対称でるとされている (14, 19) . そこで, 適切な λ の設定値についてあることは興味深い.

検討するため, これを推奨値の 1, 3, 5, 10 倍 (対応する個体3.2. フォノニック導波路の最適設計

数は各々λ = 6, 18, 30, 60) として, 各々100 回最適化計算を次に, 提案法を用いてフォノニック導波路を最適設計した

実行した. Table 1 に, S 行列法により加速した境界要素法と例を示す. トポロジー最適化 (8) により得られた 2 次元フォ

CMA-ES に基づく提案法による 100 回の試行のうち global

ノニック構造 (周期 1 とし, 第一フルバンドギャップを最大化な最適解 ϕ∗ に収束した回数及び ϕ∗ あるいは ϕ∗∗ に収束したした構造, ω ∈ [2.5, 4.0] 付近にフルバンドギャップを持つ) の回数を示す. λ を推奨値に設定した (λ = 6) 場合には半数程度

周期単位を 7 × 11 個並べたものから 8 個取り除いたフォノの試行において局所解に陥ったことが分かる. λ を大きく設定

ニック導波路を考える (Fig. 6). 図に青で示す散乱体を各々することで global な最適解に収束する確率は改善し, λ = 30の重心まわりに ϕi [◦ ] (i = 1, 16) 回転させ, フルバンドギャッに設定することで 8 割程度の「成功」率となるが, λ = 60 とプの中心の角周波数 ω = 3.25 を持つ x1 方向正の向きに入射しても更なる改善は得られなかった. また, λ = 18 以上に設した単位振幅の平面波を欠陥部に伝播させるような構造を− 69 −

4

2

0

-2

-4
Fig. 4

Sound intensity around the optimised objects with ϕ

∗

-3

-2

-1

0

1

2

3

(The global optimum ob-tained by setting λ = 30). The number

Fig. 6of design variables in this case is 2.

Setting for the phononic waveguide design. Blue scatterers are rotated so thatthe sound intensity on the observation
points is maximised.

では, 各ステップにおいて λ = 61 種類の設計変数の組み合わせに対して目的関数を評価するが, Fig. 7 は各ステップにおいて最大であったものを plot した. 最適化の初期のステップ
においては目的関数が振動しているが, 徐々にその幅は小さくなり, 最終的に −J = 191 に収束した. Fig. 8, Fig. 9 に, 初期形状, 最適化後の散乱体配置とその周辺の音場強度を各々示す. 初期形状に対しては「曲がり角」付近までしか音波が到達していない一方で, 最適化後は欠陥部に音場が局在しつ

Fig. 5

つ, 「出口」付近まで音波が到達している様子が見て取れ, 提Sound intensity around the optimised ob-

案する手法でフォノニック導波路を設計することに成功した

jects with ϕ∗∗ obtained by settting λ =30. The number of design variables in thisと言える. ただし, 前節の結果からも明らかであるが, ここでcase is 2.

得られた設計案が大域的最適解である保証はない. なお, S
行列を用いた提案法および S 行列を用いない従来の境界要素
法 (代数方程式の解法には許容誤差を 10−5 とした GMRES

探索する. この目的のため, [2.5, 3.5] × [−1.0, 0.0] に等間隔に11 × 11 点配置した観測点 xobs
(i = 1, · · · , 121) での音場強度iを用いた) を用いた場合の最適化の第一ステップ (=目的関数を λ = 61 回計算) に要した時間を Table 3 に示す. 提案法はの −1 倍を目的関数に設定した. 目的関数の計算においては,従来法と比較して効率的であることが分かる.

単一の散乱体の表面を 121 分割した境界要素を用いた (波長4. 結言と最も長い要素の長さの比は 75). λ 以外の全ての戦略パラ

本研究では, 多点探索に基づく連続最適化法の一つである

メータに推奨値を用い, 前節の結果を踏まえ λ は推奨値の 5

CMA-ES と多重散乱問題の効率的な解法である S 行列法を倍 (=61) に設定した. また、単位構造の持つ対称性から, す組み合わせた, 多数の剛体からなるフォノニック構造のパラべての設計変数に制約 0◦ ≤ ϕi ≤ 90◦ を課した. なお、77 個メータ最適化法を提案した. 公開されている実装 (22) を用

全ての散乱体の角度を最適化する問題を CMA-ES で取り扱いることで, 容易に最適化システムを構築することができた.

うことも可能であるが、ここでは、フォノニック周期構造に応用例として, フォノニック導波路の最適設計を行い, 周期単乱れを導入することで導波路を実現するという物理的な背位の回転という簡単な操作のみで所望の周波数の音波に対景に基づき、導波路付近の散乱体の角度のみを設計変数とする導波路が設計できることを示した. 一方で, 個体数を多した。

く設定したとしても大域的な最適解に収束しない場合も見Fig. 7 に目的関数 (の −1/121 倍) の推移を示す. 提案手法

− 70 −られた. 今後, 収束性を改善するための CMA-ES の改善案と

2

objective function

1.5

1

0.5

0
0

Fig. 7

50

100

150
step

200

250

300

History of the averaged sound intensity oneach observation points. The maximum
value (among 61 oﬀsprings) at each stepis plotted against the iterative step of optimisation. The number of design variablein this case is 16.

Fig. 8Soundintensity

around/withinthephononic waveguide obtained by removTable 3

Computational time for one iteration

ing several scatteres from a phononicof the optimisation with the S-matrix-

structure desinged by a topology optimi-

accelerated and naive boundary element

sation(8) .

method. In the naive one, the algebraic
equations are solved by GMRES withtolerance 10−5 . The computations ran

elastic composites. Physical Review Letters, Vol. 71,

on a PC with Intel(R) Xeon(R) GoldNo. 13, p. 2022, 1993.

5118 using 24 cores.

(2) Z. Liu, X. Zhang, Y. Mao, Y. Y. Zhu, Z. Yang,

S 行列法

素朴な境界要素法

20.612302.27

C. T. Chan, and P. Sheng. Locally resonant sonic materials. Science, Vol. 289, No. 5485, pp. 1734–6, 2000.

計算時間 [秒](3) O. Sigmund and J. S. Jensen. Systematic design ofphononic band–gap materials and structures by topol-

して bi-population 型のアルゴリズム (26) を用いることが考えられる. また, 更なる計算効率の改善のため, CMA-ES を少ないステップ数で打ち切って得られた解を初期値として準ニュートン法を実行することも有効であると考えられる.

ogy optimization.

Philosophical Transactions of theRoyal Society of London. Series A: Mathematical, Physical and Engineering Sciences, Vol. 361, No. 1806, pp.
1001–1019, 2003.

本論文では設計変数の数が 2 または 16 の場合の結果を示(4) H. W. Dong, X. X. Su, Y. S. Wang, and C .Zhang.

すにとどまったが、CMA-ES ではより高次元の最適化問題をTopological optimization of two-dimensional phononic(19)

crystals based on the finite element method and genetic

取り扱うことが可能である。より大規模な問題へ本手法を適用することも課題の一つである。

algorithm. Structural and Multidisciplinary Optimiza-

謝辞

tion, Vol. 50, No. 4, pp. 593–604, 2014.

本研究は, 科学研究費補助金 (17K14146, 19H00740) の助
成のもとに行ったものである.

(5) J. S. Jensen and O. Sigmund. Topology optimization ofphotonic crystal structures: a high-bandwidth low-loss
T-junction waveguide. Journal of the Optical Society of

参考文献

America B, Vol. 22, No. 6, pp. 1191–1198, 2005.

(1) M. S. Kushwaha, P. Halevi, L. Dobrzynski, andB. Djafari-Rouhani. Acoustic band structure of periodic− 71 −(6) T. Yamada, K. Izui, S. Nishiwaki, and A. Takezawa.
A topology optimization method based on the level set
method incorporating a fictitious interface energy. Com-

rial structures. IEEE Journal on Multiscale and Multiphysics Computational Techniques, Vol. 4, pp. 234–244,2019.
(13) N. Hansen and A. Ostermeier. Adapting arbitrary normal mutation distributions in evolution strategies: The
covariance matrix adaptation. In Proceedings of IEEE
international conference on evolutionary computation,
pp. 312–317, 1996.
(14) N. Hansen. The cma evolution strategy: A tutorial.
arXiv preprint arXiv:1604.00772, 2016.
(15) T. Lutz and S. Wagner. Drag reduction and shape optimization of airship bodies. Journal of Aircraft, Vol. 35,
No. 3, pp. 345–351, 1998.
(16) A. L. Marsden, M. Wang, J. E. Dennis, and P. Moin.
Optimal aeroacoustic shape design using the surrogate
management framework. Optimization and Engineering, Vol. 5, No. 2, pp. 235–262, 2004.
(17) E. Raponi, M. Bujny, M. Olhofer, N. Aulig, S. Boria,
Fig. 9Soundintensity

around/within

theand F. Duddeck. Kriging-assisted topology optimization

phononic waveguide (optimised by using
61 oﬀsprings).

of crash structures. Computer Methods in Applied Me-

The number of design

chanics and Engineering, Vol. 348, pp. 730–752, 2019.

variable in this case is 16.
(18) G. Fujii, M. Takahashi, and Y. Akimoto. CMA-ESbased structural topology optimization using a level set
puter Methods in Applied Mechanics and Engineering,boundary expression―application to optical and carpet

Vol. 199, No. 45, pp. 2876–2891, 2010.

cloaks. Computer Methods in Applied Mechanics and(7) 飯盛浩司, 高橋徹, 松本敏郎. B スプライン曲面のレベルセットを用いたトポロジー最適化. 計算数理工学論文(19) 高橋正幸, 秋本洋平, 藤井雅留太. 共分散行列適応進化
戦略に基づいた音響クロークのトポロジー最適化. 日本

集, Vol. 17, pp. 125–130, 2017.

(8) 釜堀瑞生, 飯盛浩司, 高橋徹, 松本敏郎. SS 法と境界要素法を用いた周期構造のトポロジー最適化. 計算数理工学
論文集, Vol. 16, pp. 85–90, 2016.

機械学会論文集, pp. 17–00590, 2018.

(20) L. L. Foldy. The multiple scattering of waves: I. general
theory of isotropic scattering by randomly distributed(9) G. Yi and B. D. Youn. A comprehensive survey on
topology optimization of phononic crystals. Structural
and Multidisciplinary Optimization, Vol. 54, No. 5, pp.

scatterers. Physical review, Vol. 67, No. 3-4, p. 107,1945.
(21) Z. Gimbutas and L. Greengard.

Fast multi-particle

scattering: A hybrid solver for the Maxwell equations

1315–1344, 2016.
(10) G. Allaire, C. Dapogny, R. Estevez, A. Faure, and
G. Michailidis. Structural optimization under overhang
constraints imposed by additive manufacturing technologies. Journal of Computational Physics, Vol. 351,
pp. 295–328, 2017.

in microstructured materials. Journal of Computational
Physics, Vol. 232, No. 1, pp. 22–32, 2013.
(22) https://github.com/beniz/libcmaes.
(23) G. A. Jastrebski and D. V. Arnold. Improving evolution
strategies through active covariance matrix adaptation.

(11) 山田崇恭, 正宗淳, 寺本央, 長谷部高広, 黒田紘敏. 幾何学的特徴量に対する偏微分方程式系に基づく幾何学的特徴制約付きトポロジー最適化 (積層造形における幾何学的特異点を考慮したオーバーハング制約法). 日本機械学会論文集, Vol. 85, No. 877, pp. 19–00129, 2019.

(12) B. Blankrot and C. Heitzinger.

Engineering, Vol. 332, pp. 624–643, 2018.

In 2006 IEEE international conference on evolutionary
computation, pp. 2814–2821, 2006.
(24) F. A. Amirkulova and A. N. Norris. Acoustic multiple scattering using fast iterative techniques. In ASME
International Mechanical Engineering Congress and Ex-

Eﬃcient computa-

tional design and optimization of dielectric metamate-− 72 −

position, Vol. 58486, p. V013T01A005, 2017.

(25) R. Coifman, V. Rokhlin, and S. Wandzura. The fast(26) A. Auger and N. Hansen. A restart CMA evolution

multipole method for the wave equation: A pedestrian

strategy with increasing population size. In 2005 IEEE

prescription. IEEE Antennas and Propagation maga-

congress on evolutionary computation, Vol. 2, pp. 1769–

zine, Vol. 35, No. 3, pp. 7–12, 1993.

1776, 2005.

− 73 −