2013年９月4日
平成２５年度電気学会C部門大会
先端技術セミナー

電磁界解析による最適設計技術五十嵐 一
北海道大学 情報科学研究科
システム複合情報学研究室
http://hbd.ist.hokudai.ac.jp

計算知能による設計とは人工知能の進展：日産ミーブの自動運転
コンピュータが将棋に勝った
IBMワトソンがアメリカのクイズ番組で優勝
顔認識(USJ)

設計に適した方法は？
エキスパートシステム：どのように知識を表現し，活用するか
計算知能：
計算(シミュレーション)により知能処理を実現する

電磁界解析と進化論的最適化
1970年代1st COMPUMAG 1976
Weiland, 有限積分法(FIT) 1977

Holland “Advances in Naturaland Artificial Systems” 19781980年代

中田・高橋 1982
「電気工学の有限要素法」
Silvester&Ferrari “Finite elementsfor electrical engineers” 1983
Bossavit, 辺有限要素法 1988

Goldberg “Genetic Algorithm in
Search, Optimization and Machine
Learning” 19891990年代３次元解析
GAによる最適化
Berenger, Perfect Matched Layer

Michalewicz， “ Genetic Algorithm+ Data Structure=Evolution Programs”

1994

IEEE Trans. Evolutional Computation

1994

1997

2000年代

商用ソフトの普及

Multi-objective, Particle Swarm…

Genetic Algorithm + Computational Electromagnetism = ?

Maxwell 方程式
D
rotH  J t
BrotE  
tdivD  
divB  0
D  EB  H

電磁界解析手法のまとめ

有限要素法

FDTD (FIT)モーメント法原理

重み付残差法Yeeアルゴリズム

積分方程式

主な用途

低周波Δt, jω
高周波 jω

高周波 Δt

長所

形状表現能力

時間陽解法

無限領域

短所

時間陰解法
無限領域

形状表現能力不均一媒質

高周波 jω

電磁界の離散化
電界 E(x) を節点ベースの補間関数で近似する方法
nE( x )   E i N i ( x )i 1EzEyEx

有限要素要素の界面で E の全成分が連続となる
E の法線成分の不連続性を表せない
後に述べる可換性も満足できない

電磁界の境界条件
 H  dl   H  dl

c1c2

c1

1 , 1c1

 2 , 2c2

 D  dS   D  dS

S1

 E  dl   E  dlc2

 B  dS   B  dS

S2S1

1 , 1S1

 2 , 2S2S2

直接計測可能な量を変数とすべし
ei   E  dl
起電力ci

起磁力hi   H  dlci

磁束bi   B  dS
Si

電束アンペアの法則 (静磁界)di   D  dSSi

h1  h2  h3  j1
ファラデーの法則

db1e1  e2  e3  dt

c3
S1
c1c2

勾配行列 G と 回転行列 Rn2(節点接続行列と閉路行列)1 ( 1)  Ci が nj に入る (出る)Gij   Ci が njを含まない
0+

c1+n1
c1  n1  n2

  G1 j n jjc2

1 (1)  S1 が cj と平行 (反平行)Rij  0 S1 が cj を含まない
RG  0  rot grad  0

c3

S1

c1c4
S1  c1  c2  c3  c4

  R1 j c jj

離散化により保存される性質
連続系

rot grad u  0

離散系

RGu  0c2c3境界の境界はゼロ

Si    Rij c j   Rij G jk nk  0j

 RG  0
同様に
連続系

div rot A  0

離散系

DRa  0j ,kS
c4c1

マクスウェル方程式の離散化
アンペア則

ファラデーの法則

D
rotH t

微分形

 H  dl   D  n dS

積分形

S

離散計

h3

S

rotE  Bt E  dl    B  n dS

SSRh  dRe  b

h2e2d1h4h1e3b1e4e1

FDTD法 (FIT)
sub-cellh2

h3

d3h4

h1

e3n 1 / 2n 1 / 2h
h[ R] e n  [  ]te2b1e4
primary celle1n 1nee[ R]t hn1/ 2  j n1 / 2  []t有限要素法の基底関数が満足すべき条件i

直交性

 N e  dl   eeeeae   A  dl

A   ae N ee

可換性ejrotN   R t [ M ]A   ae N eerot

B   ae rotN eeB   bs M s   Rse ae M sss ,ee’

有限要素法の定式化重み付き残差法

 w ( x )  [rot( ( x )rotA( x ))  J ( x )] dV  0



残差

部分積分により弱形式を得る

 ( rotw  rotA  w  J ) dV  0a2

辺要素による近似A( x )   ai N i ( x )ガレルキン法

a1a6iw  N i , i  1,2,..., e

a3

a5a4

有限要素方程式
 a j  ( rotN i  rotN j ) dV   N i  J dV , i  1,2,.., ej

行列形式ではrot  rot A  J

R t NRa  j
特異行列

Ker (R )  G

Ker (rot )  grad

spanning treeの変数は冗長a2a1a6a3a5a4静磁界の有限要素解析
連続系

離散系

rot  rot A  J

R t NRa  j
特異行列

Ker (R )  G

Ker (rot )  grad

spanning treeの変数は冗長
a2a1

a6a3a5a4正則化は行列条件の悪化を招く
正則行列
特異行列





n



n

n

条件数

n

maxcond(K ) minH. Igarashi, "On the Property of the Curl-Curl matrix in Finite Element Analysis
with Edge Elements", IEEE Trans. Magn., vol. 37, no.5, pp.3129-3132 (2001)山登り法と集団的探索

遺伝子形と表現形
遺伝的アルゴリズム
x31011110111

x4x2f (x1 )x1 x2 x3 x4 上形状 下形状
x1f (x 2 )0100000000

遺伝子形

表現形

適合度

交 叉
x3

1011110111

x4x2x1 x2 x3 x4 上形状 下形状
x1f (x1 )1011110000

1011110111

f (x 2 )f (x1 )

子1

親1

親2

f (x2 )子20100000000

0100000111

電磁界解析を用いた最適化
手法

特徴

課題

パラメータ最適化

幾何パラメータや材料パラメータなど特徴量で表現

有限要素の変形や生成が必要

トポロジー最適化

ビットマップや関数で形状を表す

複雑形状

ICタグアンテナの最適化
目的：通信距離を最大化する．
ICとアンテナのインピーダンス整合を実現
アンテナゲインを最大化
基準線

263～

表現型 (Phenotype)2 6 3遺伝子型(Genotype)
H. Makimura, Y. Watanabe, K. Watanabe, H. Igarashi, Evolutional Design of Small Antennas for Passive UHF-Band RFID,IEEE Trans.Magn., vol.47, no.5, pp.1510-1513, 2011.5
牧村，渡辺，五十嵐，和木： 電波型パッシブRFIDによる電車線路設備のモニタリングシステム，
電気学会論文誌C，vol.132, no.5, pp.691-696, 2012

ICタグアンテナ最適化結果

免疫アルゴリズムによる最適化過程

実測結果

ICタグ用アンテナの最適化
half wavelength dipole @956Mhz
157

4120MLA

85,6

credit card size
(ISO/IEC 7810 ID-1)

12

PLA

53,98

54eMLA963(unit : mm)
本研究の一部はJR東日本(株)との共同研究による

インダクタ最適化

インダクタの断面写真
フェライトコアと巻き線

市場規模：約1兆円
日本シェア6割
小型化，高効率化，高性能化が求められている．

パラメータ最適化
最適化パラメータ：Y3x  ( X 1, X 3, Y 1, Y 2, Y 3) tY2
HY1LACL AC1  100 H
L AC 2  80 Hcoil
X1(X2) X3

WLAC1

LAC20.2

1.0

1/4モデル
最適化問題

f ( x ) | LAC1 106  100 | W 2 H  penalty
0 ( LAC2 > 80[μH]のとき)
penalty   80  L AC 2  106 (それ以外)

I [A]パラメータ最適化結果

最適化前

IAによる最適化結果

micro-GAによる最適化結果

K. Watanabe, F. Campelo, Y.Iijima, K. Kawano, T. Matsuo, T. Mifune, H. Igarashi,
Optimization of Inductors Using Evolutionary Algorithms and Its Experimental Validation,
IEEE Transactions on Magnetics,
vol. 46 no. 8, pp.3393 - 3396

パラメータ最適化の検証
120

インダクタンスＬ［Ｈ］

100

解析結果
試作結果

80604020

試作したインダクタの断面写真

0012
ＤＣバイアス電流［Ａ］

本研究は一部はJST産学協同シーズイノベーション化事業および太陽誘電(株)との共同研究による
K. Watanabe, F. Campelo, Y.Iijima, K. Kawano, T. Matsuo, T. Mifune, H. Igarashi,
Optimization of Inductors Using Evolutionary Algorithms and Its Experimental Validation,
IEEE Transactions on Magnetics,
vol. 46 no. 8, pp.3393 - 3396

3

形状変化の表現法
手法
メッシュ生成メッシュ変形*
ボクセルFEM

特徴
新たに有限要素を生成
有限要素を変形
ビットマップ状に形状を表現課題計算量 O(nlogn)
要素の裏返り材料境界が階段状

29

階層ボクセル有限要素法
材料境界近傍のボクセルを細分化することで，表現能力形状を向上

time (s)128voxel

nlogn

O(nlogn)O(n)4050000

100000 150000
number of nodes

200000T. Sato, H. Igarashi, Coupled Analysis of Electromagnetic Vibration Energy Harvesterwith Nonlinear Oscillation, to appear in IEEE Trans. Magn.

30

非適合接続要素間で磁束が連続になるように，スレーブ変数をマスター変数の線形結合により表現する．
a s   A  dl   a j  N j  dlesjes

 c1a1  c2 a2  c3a3  c4 a4a4a1a2
asa3
31

インダクタの３次元最適化目的：ジュール損失とヒステリシス損失の最小化

f 2  Wh ,Objectives: f1  R,
constraints:L  LT 2  0.03,
LT 2Ll  0.7 LT ,r2  w  1.5mm
L=1μH
Vin=3.3V
RIC
Cin=2.2 μF

Vout=5 V

Control logic

RFB
Cout
=1MΩ =4.7 μF
PWM signal,
3.5MHzBoost circuit.

RoR: 10.5mΩ, Wh: 9.08mW

hysretesis loss (mW)

パレート最適解

solutions (0.15)12solutions (0.1)

test model

R: 16.1mΩ, Wh: 4.8mW

108(A)

6(B)

(C)

420.1

0.12
0.14
0.16
winding resistance (Ω)

0.18R: 12.07mΩ, Wh: 6.67mW
T. Sato, et al., 3-D Optimization of Ferrite Inductor Considering Hysteresis Loss,
IEEE Trans. Magn., vol.49, no.5, pp.2129-2132, 2013

磁気シールドの３次元最適化
Fz

BtrWk min .
BT
WT
Coil,
1A, 100turn

1000

unit: (mm)zunit:[mm]
18.65
31.1
5.3

140

39.4

120
c4
c3
c2
c1

120

140
32.95

Iron

36.3
6.25

20
TargetRegion

x20

a1

a2 a3

a4

39.4

1000

xy
1000
120

140

140

120
a4
a3
a2
a1
20
TargetRegion

x20

a1a2

a3 a4

1000

(a):最適化形状：
Bbr:0.16, W:5908.99(b):等厚形状：
Bbr:0.337, W:5957.9

トポロジー最適化
On-Off法:セルの２値状態(材料あり，なし)を決める

f ( x)  +w

TTT
 w max min  min .
T0T+w
-v
-v

27.5mm+u+u8mm
56mm

IPM*のロータ構造の最適化*電気学会技術報告，no.776, 回転機のバーチャルエンジニアリングのための電磁界解析技術

空間フィルタの適用

トルク平均(Nm)

2.195

トルク平均(Nm)

3.728

リップル

0.136

リップル

0.103

目的関数

-0.167

目的関数-0.647

適応フィルタリング
 change in f filtering  1  exp  temperature 

トルク平均(Nm)

3.728→ 4.319

リップル

0.103 → 0.0767

目的関数

0.647 →0.870

日高，五十嵐，計算工学会2013にて発表予定

トルク[Nm]

適応フィルタ

固定フィルタフィルタなし回転角度[degrees]レベルセット法
レベルセット関数x
 d  x,  
0x   x    d  x ,   x  

材料中材料境界材料外d ( x )  min ( x , y )y









レベルセット法
静磁界の有限要素方程式：

Ka  jレベルセット関数の修正

df
f f a  t a df : 目的関数
随伴変数法

~f ( x )  f ( x )  z t ( j  Ka )レベルセット関数 の等高線~t a t K
f f  f ta
    z K   z   a  fKz a２段階法：On-Off法とレベルセット法
Global Search with GA
Generations

01nLocal Search based on Level set approach
Steps01n41

２段階法による磁気シールド形状の最適化
f    WM

B average

105

SS design min .

42

On-Off法による大域探索

43

レベルセット法による局所探索

44

A.P. Wills, On the magnetic shielding effect of trilamellar
shperical and cylindrical shells, Phys. Review, vol.4, 193, 1899

46

最適化結果
On-Off 法

On-Off + レベルセット|B|average /10-50.124|B|average /10-50.0990

Volume of the core (cm2)2.736

Volume of the core (cm2)2.776

Objective function

0.0678

Objective function

0.0631

Y. Hidaka, T. Sato, H. Igarashi, Topology Optimization Method Based on On-Off Method and Level Set Approach,presented at COMPUMAG2013, to appear in IEEE Trans. Magn.

47

解の変化
F    WM

B average
10 5SS design

48

On-Off法と2段階法の比較
On-Off method

Torque average (Nm)

5.280

On-Off + Level Set method

Torque average (Nm)

5.309

Torque ripple

0.184

Torque ripple

0.112

Objective function-0.806

Objective function-1.018Y. Hidaka, T. Sato, H. Igarashi, Topology Optimization Method Based on On-Off Method and Level Set Approach,presented at COMPUMAG2013, to appear in IEEE Trans. Magn.

On-Off法と2段階法の比較

49

Y. Hidaka, T. Sato, H. Igarashi, Topology Optimization Method Based on On-Off Method and Level Set Approach,presented at COMPUMAG2013, to appear in IEEE Trans. Magn.

フラックスバリアの効果
Non Flux Barrier

Flux Barrier50２段階法による磁気コア形状の最適化
y[mm]

Coil
2.0×106[A/m2]

DesignRegion

点(90, 40)

60
TargetDomain30202040

90 100x[mm]On-Off法と2段階法の比較
On-Off + レベルセット

On-Off 法

目的関数67.480

65.198関数重ね合わせによるトポロジー最適化
NGnet(1)を用いる方法(2)Gi x Gˆ i x   N

 Gk x k 1

NGnet outputs

i 10.20-10-50input

Normalization

510-10-5-1

0510input

1.51

Weight sum

0.50-10Gi(x): ガウス関数

210Basis outputs

y x    wi Gˆ i x 30.4

Gaussian outputs

N-50input

510(1)J. Moody: "Fast Learning in Networks of Locally-Tuned Processing Units,"
Neural Computation, Vol. 1, No. 2, pp. 281-294, 1989.
(2)佐藤孝洋，渡辺浩太，五十嵐一，“正規化ガウス関数ネットワークを用いた電気機器のトポロジー最適化法”，静止器回転機合同研究科，SA-13-064,RM-13-078，2013

ガウス基底の配置

最適化の過程
-0.4
100

200

-0.6

-0.8

-1(10)

世代数
トルク

リップル

トルク (Nm)

3.5

0.5

32.5

0.3

21.5

0.1
0

100
世代数

200

リップル

目的関数値

0最適化結果

4optimized

original

トルク(Nm)3.5
32.5
Torque ( Nm)

2
0

10

2030

回転角(度)
トルクの変化

最適化されたフラックスバリアの形状

本研究の一部は(株)明電舎との共同研究による
佐藤孝洋，渡辺浩太，五十嵐一，“正規化ガウス関数ネットワークを用いた電気機器のトポロジー最適化法”，静止器回転機合同研究科，SA-13-064,RM-13-078，2013

コンピュータによる設計
設計原則： 設計者はコンピュータが設計したものについてよく理解していなければならない．

弱い設計原則：設計者はコンピュータどのような原理でものを設計するのかを理解していなければならない．

今後の課題(1) 現実の設計に適用するために，さまざまな要因を考慮して最適化を行える技術を構築する(2) トポロジー最適化に関する理論を構築する(3) 電磁界解析の高速化，最適化の効率化(4) ものづくりとの協調