統計数理（2005）
第 53 巻 第 2 号 297–315
c 2005 統計数理研究所


特集「計算推論―モデリング・数理・アルゴリズム―」［研究詳解］

2 次錐計画問題による磁気シールドのロバスト最適化1

2

土谷 隆 ・笹川 卓（受付 2005 年 4 月 27 日；改訂 2005 年 10 月 4 日）

要

旨

本論文では凸計画の新潮流の代表例の一つである 2 次錐計画を取り上げ，超電導磁気浮上式
列車の磁気シールド最適設計問題への応用を紹介する．超電導磁気浮上式列車は超電導磁石によって浮上・推進するため，車体を強磁性体のシールドで覆い，車内を強力な磁場から遮蔽する必要がある．一方，車体はできるだけ軽いことが望ましい．そこで，磁場を考慮してシールドの厚みを場所によって調整し，車内へ磁場が漏洩しないという条件下でシールド重量を最小化する問題を考える．適切なモデル化により，この問題は 2 次錐計画問題として定式化できる．
本論文では代表的な内点法の一つである主双対内点法を実装して 2 次錐計画問題を解くことにより最適設計を行う．そして，さらに解法の高速性と頑健性を生かして確率的シミュレーショ
ンによるロバスト最適化を行った結果，および最尤法によるロバスト性能の事後評価について報告する．
キーワード： 2 次錐計画，内点法，超電導磁気浮上式列車，最適設計，ロバスト最適化，磁気シールド．
1.

はじめに

最適化の世界では，長らく線形計画問題のみが実用レベルで確実に解ける問題であるとされ，線形計画法を中心として分野が発展してきた．しかしこの 10 年あまりの間に内点法の登場や計算機の高速化に伴い一般の凸計画問題の実用化が進み，制御や信号処理，機械学習，統計や組合せ最適化など諸分野への応用が着実に広がっている．その代表例が 2 次錐計画問題と半正定値計画問題である．
線形計画問題は，アフィン空間と第一象限という凸錐の交わり上で線形目的関数を最適化する問題として定式化できる．2 次錐計画問題（Ben-Tal and Nemirovski, 2001; Lobo et al., 1998;
Nesterov and Nemirovskii, 1994）と半正定値計画問題（Alizadeh, 1995; Ben-Tal and Nemirovski,
2001; Nesterov and Nemirovskii, 1994; Wolkowicz et al., 2000）は線形計画問題における “第一象限” をより複雑な “2 次錐や半正定値対称行列錐という特別なクラスの凸錐” に置き換えた凸計画問題である．これらの問題は，対称錐計画問題というより一般的な凸最適化問題のクラスの一例となっている．対称錐計画問題は 21 世紀の線形計画問題などとも呼ばれ，古典的線形計
画問題の素直な（しかし驚くべき）拡張となっており，線形計画問題の有する種々のよい性質を
12

統計数理研究所：〒106–8569 東京都港区南麻布 4–6–7
鉄道技術総合研究所：〒185–8540 東京都国分寺市光町 2–8–38

298

統計数理 第 53 巻

第 2 号 2005

引き継いでいる．特に日本で生まれ，線形計画問題を解く際には今や定番の一つとなっている
多項式時間解法である主双対内点法（Kojima et al., 1989; Tanabe, 1988; Wright, 1997）が，その主要な理論的結果も含めて拡張可能であることは，対称錐計画問題の特筆すべき性質の一つである．本論文では，対称錐計画問題の典型例である 2 次錐計画問題の具体的応用として，磁気浮上式列車の磁気シールドの最適設計問題を取り上げ，得られた研究成果を Sasakawa andTsuchiya（2003），笹川（2005）を元にして紹介する．
現在開発が進められている超電導磁気浮上式列車（鉄道技術総合研究所, 2005）は車上の超電導磁石と地上の常電導コイルの電磁相互作用によって浮上・推進する．その際，車体を強磁性体（たとえば純鉄）からなるシールドで覆い，車内を超電導磁石が作る強力な磁場から遮蔽する必要がある．一方，車体はできるだけ軽く作りたいという要請がある．そこで，磁場を考慮してシールドの厚みを場所によって調整し，磁場が車内に洩れないという条件の下でシールドの重さを最小化することを考える．適切なモデル化を行うとこの問題は 2 次錐計画問題に定式化できるので，Nesterov と Todd，筆者の一人によって開発された 2 次錐計画法に対する主双対内点法によって多項式時間で求解が可能である（Nesterov and Todd, 1997; Tsuchiya, 1999）（Monteiro and Tsuchiya, 2000; Schmieta and Alizadeh, 2001 も参照のこと）．そこで本研究では，主双対内点法を実装してこの問題を解いて最適設計を行い，さらに解法の高速性と頑健性を生かして確率的シミュレーションによるロバスト最適化を試みる．
ロバスト最適化（Ben-Tal and Nemirovski, 1998, 2001; Calaﬁore and Camp, 2005; El Ghaoui
and Lebret, 1997; El Ghaoui et al., 1999; Goldfarb and Iyengar, 2003）はモデルの不確定性まで考慮した上でロバストな最適化を行うより高次の最適化技術である．本論文においては数千変数の 2 次錐計画問題を主双対内点法により 10000 回繰り返して解くことによりロバスト最適化を行い，さらに設計されたシールドのロバスト性能を最尤法による統計的手法で評価する．
我々の実装ではこの計算を通常のパソコンで数時間で実行できる．これは，高速で安定した最適化アルゴリズムの登場によって，数千変数の凸最適化問題の求解の繰り返しを含むような新
しい種類のシミュレーションや高次元積分が可能となりうることを示している．
本論文の構成は以下の通りである．2 節で磁気シールド最適設計問題を導入する．3 節では，
2 次錐計画問題とそのアルゴリズムについて簡単に説明する．4 節では，最適化を行った結果について述べる．5 節では，ロバスト最適化について述べる．6 節はまとめである．
2.

磁気シールド最適設計問題

磁気シールド設計問題は磁気浮上式列車の開発や，MRI，極微小磁気測定などの分野で現われる．以下では本論文でとり扱う磁気浮上式列車における磁気シールド最適設計問題について述べる（鉄道技術総合研究所, 2005; 佐々木, 1993）．磁気シールド最適設計問題のためのモデル
を物理的立場から導入し，有限要素法による離散化を行って 2 次錐計画問題に定式化する．
2.1 磁気シールド最適設計問題の定式化と 2 次錐計画問題
現在開発が進められている超電導磁気浮上式列車では，車両間の廊下（貫通路ともいう）の左右両側に超電導磁石ユニット（SCM）が配置されている（図 1 参照）（鉄道技術総合研究所, 2005）．
各超電導磁石ユニットは，直列に配置された 4 つの超電導コイル（SCC）から構成されている．
各超電導コイルは長さ 1.07m，高さ 0.5m の競技場のトラックのような形をしており，これに700kA の電流が流される．車体の左右で相対する形で配置されている 2 つの超電導コイルは同じ方向（S-N-S-N）に励磁される．この超電導コイルが生成する強力な磁場から車内の乗客を磁気遮蔽するためには，車体をシールド材（例えば純鉄）で囲む必要がある．シールド材で囲むことにより超電導コイルが生成する磁束の大部分は車体内部の空間には到達せず，シールド内を

2 次錐計画問題による磁気シールドのロバスト最適化299

図 1. 磁気浮上式列車と超電導磁石の配置（単位：ミリ）．

図 2.

問題の状況．

通ってコイルに戻る（図 2）．一方，磁気浮上式列車は空中に浮上して走行するので，シールド
重量はできるだけ軽い方が望ましい．そこで，超電導磁石ユニットが生成する磁場を考慮した
上で，シールドの厚みを場所に応じて調整し，確実に遮蔽効果が得られる中でできるだけ軽いものを設計したい．これが本論文で考察する磁気シールド最適設計問題である．以下に紹介するモデルは Sasakawa and Tagawa（1997）において提案され，現実によく合うことが確かめられているものである．
まず下記の仮定を置く．
1. 車体内部 Ω は完全に強磁性体（シールド材）で囲まれている．
2. 強磁性体の初期透磁率 µ は無限大である．
3. 強磁性体は飽和磁束密度 Bs を持つ（単位：テスラ）．

仮定 2 と 3 は，現実の強磁性体の B-H カーブを図 3 に示す形状で近似することに相当する．このように仮定することで，我々は飽和磁束密度 Bs だけを問題として非線形性の詳細を無視することになるが，結果として得られる最適設計問題は十分に元の磁気シールド問題の本質を捉えた凸最適化問題となる．
シールド材は強磁性体なので，超電導磁石ユニットにより生成される外部磁場は，シールドの厚みを変えてもほとんど変化しない．そして特に，µ を無限大にした極限，すなわち仮定 2

300

統計数理 第 53 巻

図 3.

第 2 号 2005

仮定された B-H カーブ．

の下では外部磁場（車体外部 R3 − Ω）を決定する問題はシールド材の厚みを決める最適設計問題とは独立して最初に解くことができる．また，仮定 2 の下では車体内部には磁場が漏れることはない．ただし，仮定 3 よりシールド内を流れる磁束密度には上限である飽和磁束密度 Bs が存在するので，これらの簡略化は，シールド内の磁束密度が Bs 以下であることが前提となる．
µ = ∞ とした時のモデルのこれらの性質は，車体外部，シールド，車体内部（図 2 のそれぞれ R3 − Ω，∂Ω，Ω に相当）からなる全系の磁場を定める静磁場問題を変分法的に定式化して離散化し，結果として得られる重み付き最小二乗問題において µ を無限大にした極限をとった際の解の性質として理解することもできる（笹川, 2005）．この種の重み付き最小二乗法の極限を
層別最小二乗法という（土谷, 2005）．
上の考察に基づいてまず，超電導磁石ユニットが車体外部に生成する磁場を求める静磁場問
題を（数値的に）解いて，シールドに流入あるいはそこから流出する磁束の流れの法線成分 Bn（単位：テスラ）を定めることができる．ここで静磁場問題を解くにあたってはシールド表面 ∂Ω の接線方向成分が 0 であるという境界条件を置く．
においては B次に，シールド内の磁束について考察する．シールドは板状であることより，その厚み方向
 を厚み方向に積分には磁束密度は変化しないと仮定してよい．そこで，シールド内における B
したものとして定義される，2 次元ベクトル場 F = (F1 , F2 )（単位：ウェーバー／メートル）を導入する．すると，車体内部に磁束が漏れないことより，シールド内の磁束について流量保存則(2.1)

 = Bn (∂Ω 上で)divFが成立する．ここで，シールド表面上に（局所的）デカルト座標が定義されているとすると作用
素 div は，
=
divF

∂F1
∂F2
+
∂x1
∂x2

と書ける．
以下，
（各点での）シールドの厚みを T0 で表す．T0 が本問題の設計変数である．シールド内の磁束分布 F に対する磁束密度は F /T0 で与えられる．仮定 3 より，強磁性体の飽和磁束密度は Bs であるから，各点で
(2.2)


F
≤ Bs
T0

2 次錐計画問題による磁気シールドのロバスト最適化

301が成立する必要がある．本論文では典型的なシールド材である純鉄や鉄系合金の物理的性質をもとに Bs = 1.5 テスラ とした（Sasakawa and Tagawa, 1997）．
（2.1），（2.2）は内部への漏洩磁界が存在しない状況で，T0 ，F ，Bn が満たすべき条件である．
シールドの重みはシールドの体積に比例し，シールドの体積はシールドの各点での厚みを全表面に渡って積分して得られる．そこで，磁気シールド最適設計問題を，（2.1），（2.2）の下でシールドの体積を最小化する問題として下記のように定式化する：(2.3a)

minimize

1
Bs



∂Ω

F0 dS,

 = Bn (∂Ω 上で),subject to divF

(2.3b)  ≤ F0 (∂Ω 上で).
F

(2.3c)

ここで補助変数 F0 = Bs T0 を導入した．
この問題の最適解を F0 = F0∗ ，F = F ∗ とする．この時，最適解を元にその厚みが T0 = F0∗ /Bs
で与えられるシールドを作成しても，内部を実際に流れる磁束が F ∗ となるという保証はない
ように見える．物理的にはシールド内の磁束はシールドの厚みが決まると自然と定まる筈であり，それが（2.3）の最適解の F ∗ となるかどうかは分らないからである．内部を実際に流れる磁束が F ∗ となることが遮蔽の前提となるので，この観点からモデルの妥当性を検討しておくことは重要である．
実は，実際のシールド内の磁束の流れとして F ∗ が実現されることを，“シールド内の磁束の流れは磁場のエネルギーが最小となるように定められる” という原理を用いて大略以下の通り示すことができる（笹川, 2005）．
（2.3）は連続版ユークリッドノルム和最小化問題(2.4a)

minimize

1
Bs



∂Ω

  dS,
F

subject to divF = Bn (∂Ω 上で),(2.4b)と等価である．そして，その最適解 F ∗ は，F および（2.4b）に対応するラグランジュ乗数 λ を未知関数とする Euler-Lagrange 方程式
(2.5)

∂λ
F1
=
,
F  ∂x1∂λ
F2=,
F  ∂x2

divF = Bnの F の部分である．
一方，シールド上の各点での厚みが T で与えられている時，その透磁率が十分に大きく磁束密度に依存せずに一定であるとすると，シールド内の磁束の流れは，シールド内の磁場エネル
ギー最小化問題
1Bs F 2dS,
T∂Ω(2.6a)

minimize(2.6b)

 = Bn (∂Ω 上で)subject to divFの最適解 F ∗∗  で与えられる．但し，この結果が正しいためにはシールドに十分な厚みがあり各点で磁束が飽和しないこと，すなわち条件(2.7)

が成立している必要がある．

T≥ ∗∗ 
F
Bs

302

統計数理 第 53 巻

第 2 号 2005

 ∗∗ は F
 および（2.6b）に対応するラグランジュ乗数 λ を変数とする Euler-Lagrange 方程式
F
(2.8)

F1
∂λ
,
=
T
∂x1

F2
∂λ
,
=
T
∂x2

divF = Bnの解である．
そこで，シールドの厚みを T = F ∗ /Bs ととった時の（2.6）の最適解について検討する．この時，（2.8）において F = F ∗ ，λ = Bs λ∗ と置くと，F ∗ ，λ∗ が（2.5）を満たすことより，これらが（2.8）の解となっていることが分る．すなわち F ∗∗ = F ∗ が成立する．さらに，条件（2.7）もT= ∗  F ∗∗ 
F=
Bs
Bs

より満たされていることが確認できる．したがって，（2.3）の最適解として定まる F ∗ がシール∗
ドの厚みを T = F /Bs ととった際のシールド内の磁束の流れとなる．
定式化（2.3）においては目的関数（2.3a）と制約（2.3b）は変数 F0 と F について線形である．2
番目の制約（2.3c）は線形ではないが F0 と F について凸である．この制約を 2 次錐制約といい，
本論文で重要な役割を果たす．このように，磁気シールド設計問題は凸最適化問題（2.3）として定式化できる．これは，後ほど導入する 2 次錐計画問題の連続版である．また，連続版最小費
用凸ネットワーク流問題と見なすこともできる（Taguchi and Iri, 1982）．
2.2 有限要素法による離散化
本節では前節で定式化された連続版 2 次錐計画問題を有限要素法（FEM）によって離散化す
る．まず，車体 1/4 の表面 ∂Ω を 1669 個の長方形有限要素に分割する．図 4（a）と（b）に分割を（離散化計算された）Bn とともに示す．図を見るとわかるように，車体本体と廊下部は大きさの異なる 2 つの直方体としてモデル化されている．進行方向に沿って見た時にメッシュが廊
下の 4ヶ所で細かく刻まれているところがあるが，これは超電導磁石ユニットの両端部分に対応して磁場の変化が大きくなると考えられるところである．
Gauss の定理とシールドの境界 ∂∂Ω に関する境界条件（車内を表す領域 Ω の境界 ∂Ω のさらに境界である ∂∂Ω は，境界の境界であるから列車全体を一まとめに扱うとすると本来空集合となるが，ここでは，対称性等を用いて計算対象の領域を図 4 のように縮小しているので空とはならない）を考慮すると（2.3b）の弱形式は





∂Ω



 , gradλ ) + λ Bn dS = 0(Fと書ける．ここで（,）は 2 次元ベクトルの内積，grad は ∂Ω を構成する 2 次元平面上（あるいはシールド表面）の作用素，テスト関数 λ は H 1 の任意の関数である．
本論文では，離散化後の λ としては H 1 に属する有限要素のうち一番簡単なものである長方形領域に対する双一次要素（Ciarlet, 1978）を用い，F = (F1 , F2 ) と F0 については定数要素を用∗いることにし，j 番目の有限要素上での F0 と F を各々F0j，Fj∗ と記すことにする．F0 と F の
自由度は，各要素の重心に存在するものとし，2 次錐制約をこれらの点で評価することにする．
また，Bn については，外部磁場を定める静磁場問題を解いて各要素上での関数値が与えられ
n∗ と記すことにする．
るので，これを B
そして，∂∂Ω における境界条件として，以下の通りディリクレ境界条件とノイマン境界条件を実現する（笹川, 2005）．
（1）ディリクレ境界条件 F × n = 0：これを実現するためには，離散化された λ を 0 と置く．
（2）ノイマン境界条件 F · n = 0：これを実現するためには，離散化された λ を自由にしておいて，関係する F には何の条件も付さない．

2 次錐計画問題による磁気シールドのロバスト最適化

図 4.

303

n∗ の法線成分．各長方形が有限要素でシールド表面 ∂Ω の離散化と（近似）外部磁場 B
n∗ に比例している．
ある．各要素の中央から伸びている矢印は B

我々は車体中央部断面ではノイマン境界条件を採用し，それ以外の境界ではディリクレ境界条件を採用した．ノイマン境界条件は標準的な有限要素法の境界条件の処理手続きで得られる．
一方，ディリクレ境界条件は常に満たされるわけではなく最適解（最適化されたシールド）においてのみ満たされる．その意味で “自然境界条件” に対応する．
以上より（2.3）は以下の通り離散化できる．
1 ∗(2.9a)
minimize
wj F0j,Bs j∈E(2.9b)

subject toj∈Ej∗ = (Bn∗ )i ,Dij F

i ∈ V,

304

統計数理 第 53 巻

第 2 号 2005∗
j∗  ≤ F0j
F,(2.9c)

j ∈ E,

ここで E と V は各々有限要素の添字集合およびメッシュの節点集合（ただし（1）のディリクレ
境界条件を課した境界上の節点については除く），wj は有限要素 j の面積，D は div 作用素を
n∗ )i は節点 i の試験関数と関連する Bn∗ の内積である．この問題の離散化したもの，そして (B∗
最適解における F0j が有限要素 j における離散化問題での最適なシールドの厚みを表す．
3.

2 次錐計画問題と主双対内点法

本節では 2 次錐計画問題と主双対内点法について説明する．
3.1 2 次錐計画問題
2 次錐 K(p) とは，p 次元ユークリッド空間において下記の通り定義される集合である．
K(p) = {x = (x0 , x1 ) ∈ R × Rp−1 | x20 − xT1 x1 ≥ 0, x0 ≥ 0}.

これは円錐を高次元に一般化したものであるが，自己双対性と等質性を満たす対称錐というクラスの凸錐の典型例である（Faraut and Korányi, 1994）．
x ∈ K(p) と x ∈ Int(K(p)) を x  0 と x  0 で簡単に記すことにする．K(p) の双対錐を K∗ (p)と記すことにする．K(p) は自己双対錐なので，
K∗ (p) = {s ∈ Rp | xT s ≥ 0， x ∈ K(p)}= {s = (s0 , s1 ) ∈ R × Rp−1 | s20 − sT1 s1 ≥ 0, s0 ≥ 0} = K(p)を満たす．
2 次錐計画問題は，線形関数をアフィン空間と 2 次錐の直積の交わり上で最適化する問題として定義され，その標準形は以下のように書ける：(P)

minimizenT

ci xi ,

i=1nsubject toAi xi = b, xi = (xi0 , xi1 )  0, i = 1, ..., n.

i=1

ここで Ai ∈ Rm×ki (i = 1, ..., n), b ∈ Rm , ci ∈ Rki (i = 1, ..., n). n は 2 次錐の数である．以下これを主問題と呼ぶ．線形計画や半正定値計画と同じように，2 次錐計画は多くの応用を有する．例
えば，多面体上で凸 2 次関数を最小化する凸 2 次計画問題やさらにその一般化である楕円体上で凸 2 次関数を最小化する問題も 2 次錐計画として表現できる（Lobo et al., 1998）．前節で定式化した（2.9）も 2 次錐計画問題である．
（P）の双対問題は(D)

maximize bT y,
subject to si = ci − ATi y, si = (si0 , si1 )  0, i = 1, ..., nと定義される．双対問題も 2 次錐計画問題である．
K= ni=1 ki とする．さらに，
A = (A1 A2 · · · An ) ∈ Rm×K , c = (c1 , ..., cn ) ∈ RK ,x = (x1 , ..., xn ) ∈ RK ,K = K1 × · · · × Kn ,s = (s1 , ..., sn ) ∈ RK ,

2 次錐計画問題による磁気シールドのロバスト最適化

305と記すことにすると,（P）と（D）は線形計画問題に類似の形で以下のように書ける：(P) minimize cT x,
subject to Ax = b, x  0,(D) maximize bT y,
subject to s = c − AT y, s  0.

ここで  の意味は自明であろう．半直線は 1 次元 2 次錐なので線形計画問題は 2 次錐計画問題の特殊な場合となる．主問題と双対問題の任意の許容解 x と (s, y) においてcT x − bT y = xT (c − AT y) = xT s ≥ 0,が成立する．ここで最後の不等式は x ∈ K，s ∈ K∗ = K から従う．このように，主問題の目的
関数値は双対問題の目的関数値よりも値が大きいか等しい．xT s(= cT x − bT y) を双対ギャップという．さらに，もし x と (s, y) が双対ギャップが 0 であるような許容解であれば，すなわちもし (x, s, y) が(PD)

sT x = 0, Ax = b, s = c − AT y, x  0, s  0,を満たすならば，x と (s, y) は（P）と（D）の最適解である．
（P）と（D）が内点許容解を持てば最適解は常に存在することが知られている．
なお，錐の形を変えない線形変換は群を成す．この群を錐の自己同型群という．2 次錐計画
問題は，錐の自己同型群の元による変換に関して問題の形が不変であるという大きな特徴がある．この変換の自由度は性質の良いアルゴリズムを作る上で重要な役割を果たす．
3.2 中心曲線と主双対内点法
主双対内点法は主問題（P）と双対問題（D）の許容領域の内部に定義され最適解に至る曲線を辿って問題を解く方法である．この曲線は中心曲線と呼ばれ，以下に述べるユークリッド的ジョ
ルダン代数を用いて定義される（Faraut and Korányi, 1994; Faybusovich, 1995, 1997; Tsuchiya,
1999）．Rki の元 xi = (xi0 , xi1 ) ∈ R × Rki −1 と si = (si0 , xi1 ) ∈ R × Rki −1 の積 xi ◦ si を以下のように定義する：xi ◦ si = (xTi si , xi0 si1 + si0 xi1 ).
Rki 上の 2 次錐 K(ki ) に関連するユークリッド的ジョルダン代数は，上に述べた積が定義された ki 次元実ベクトル空間のことである（Tsuchiya, 1999）．K(ki ) は，上の積を用いて K(ki ) ={v | v = w ◦ w, w ∈ Rki } と表せる．ei = (1, 0, ..., 0) は単位元である．ユークリッド的ジョルダン代数は下記のように積を定義すれば容易に全空間に拡張可能である：x ◦ s = (x1 ◦ s1 , ..., xn ◦ sn ).

この場合 e = (e1 , ..., en ) が単位元となる．
(3.1)

eT (x ◦ s) = T

ei (xi ◦ si ) =i Tx i si = x T siTであるから，x ◦ s = 0 ならば x s = 0 が成立する．したがって x ◦ s = 0 を満たす（P）と（D）の許容解を求めれば 2 次錐計画問題が解ける．
（P）と（D）の中心曲線は ν ∈ (0, ∞] をパラメータとする下記のような双一次方程式の解の成す集合として定義される：(3.2)

x ◦ s = νe, Ax = b,

AT y + s = c, x  0, s  0.

306

統計数理 第 53 巻

第 2 号 2005

主問題（P）と双対問題（D）に内点許容解が存在するという条件の下で，（3.2）は許容領域の内部に滑らかなパスを定義し，そのパスは ν → 0 で最適解に収束する（Faybusovich, 1997）．このパスを中心曲線という．
（3.1），（3.2），そして eT e = n が成立するので，中心曲線上では ν = xT s/nが成立する．
主双対内点法は，（P）と（D）の内点許容解の組を初期値とし，ν の値を 0 に近づけつつ（3.2）を近似的に解きながら，最適解を求めるアルゴリズムである．典型的には，(x, s, y) = (ν0 e, ν0 e, 0)(ν0 > 0) を初期点として用い，反復列は Int(K) × Int(K) × Rm の内部に生成される．
（3.2）に対するニュートン方向は AHO 方向と呼ばれる（Alizadeh et al., 1998）．多くの主双対内点法は，スケーリング付ニュートン方向といわれる方向を採用している．これは対称錐の有する対称性を生かし，錐の自己同型群の元による線形変換を行ってからニュートン方向を計算するもの
で，種々の好ましい性質を有する．スケーリング付ニュートン方向には，HRVW/KSH/M 方向（Helmberg et al., 1996; Kojima et al., 1997; Monteiro, 1997; Tsuchiya, 1999），Nesterov-Todd（NT）方向（Nesterov and Todd, 1997; Tsuchiya, 1999）などがある．
これらの方向を用いた主双対パス追跡法は一般に以下のように記述できる．
主双対内点法
精度パラメータ ε ∈ (0, 1)，ステップ幅パラメータ θ ∈ (0, 1) を設定し，(x0 , s0 , y 0 ) ∈ RK × RK ×m
R を初期点として設定する．ここで (x0 , s0 ) ∈ Int(K) × Int(K) である．µ0 := ((s0 )T x0 )/n，k := 0とする．
Repeat until µk ≤ ε do
1. (x, s, y) := (xk , sk , y k )，µ := µk とする．
2. σ ∈ (0, 1) を決める．
3. (x, s, y) におけるパラメータ ν := σµ の中心曲線上の点を求めるためのスケーリング付
きニュートン方向 (∆x, ∆s, ∆y) を計算する．
4. ステップ幅 α > 0 を，K × K の境界まで比率 θ だけ進むように選び，(xk+1 , sk+1 , y k+1 ) :=(x, s, y) + α(∆x, ∆s, ∆y) ∈ Int(K) × Int(K) × Rm とする．
5. µk+1 := ((xk+1 )T sk+1 )/n として k に 1 を加える．
End

なお，ステップ 3 において，探索方向は双対ギャップを比率 σ だけ減らすように選ばれている．
主双対内点法の実装においては，以下の 2 つが良く用いられる．
（1）基本アルゴリズム：この方法では，σ は一定値（たとえば σ = 0.1）に固定され，上述したスケーリング付きニュートン方向がそのまま用いられる．
（2） Mehrotra 予測子 修正子（MPC）アルゴリズム：この方法では σ を適応的に更新し，高次項を考慮にいれ修正した探索方向が用いられる．
MPC アルゴリズムは主双対内点法を加速するための標準的手法である（Mehrotra, 1992; Todd
et al., 1998; Wright, 1997）．
主双対内点法の一反復の主要な手間は，（i）探索方向を計算するための m 次元連立一次方程式の係数行列を計算する部分と（ii）この連立一次方程式を解くための部分に分けられる．A が密行列の場合，（i）には O(m2 K) 回，（ii）には O(m3 ) 回の四則演算が必要となる．しかし，本応用の場合は，問題の疎性を生かすことにより，ずっと高速に問題を解くことができる．なお，
探索方向の具体的な計算公式などについては Tsuchiya（1999）や Anderson et al.（2003）などを参照のこと．

2 次錐計画問題による磁気シールドのロバスト最適化

4.

307

計算結果と解析

本節では 2 節で説明した磁気シールド最適設計問題に対し，3 節で述べた主双対内点法を適用した結果を報告する．以下に報告する計算機実験はすべて Pentium III 700MHz デュアル CPU，
主記憶 1GB，OS: Windows NT 4.0 のデスクトップパソコンで行われた．
（ただし計算には 1 つの CPU しか使用していない．
）使用したプログラム言語は Fortran（Microsoft Visual FortranVer. 5 for Windows）で計算は倍精度で行われた．また，スパース行列に対するコレスキー分解のためのパッケージとしては，Visual Fortran に添付されている ISML ルーチンを用いた．
問題は既に 2 節と 3 節で説明した通りであるが，大きさについて補足する．この問題に現れる 2 次錐はすべて 3 次元のものであり，その数は 1669，そして各 2 次錐の次元は 3 なので，主問題と双対問題の変数 x，s の次元は 5007 である．また，y の次元は 1630 である．この問題は小さい同じ次元の 2 次錐を多数含み，さらに行列 A は微分作用素を離散化したものであるため疎性を有するという特徴がある．このような問題の場合，疎性を生かした実装のための技術は，大規模で疎かな線形計画問題に対する内点法のそれと類似したものとなる．なお，以下の実装
においては NT 方向を用いたが，HRVW/KSH/M 方向など他のスケーリング付きニュートン
方向を用いた場合もほぼ同様の結果が得られる（笹川, 2005; Sasakawa and Tsuchiya, 2000）．
4.1 アルゴリズムの性能
我々の実装で問題を解くのに要した時間と反復回数は，MPC アルゴリズムで 1.8 秒，21 回，
基本アルゴリズムで 3.3 秒，41 回であった．
（我々は反復を双対ギャップが 10−12 以下となったところで停止した．
）問題の大きさと計算時間の関係を調べるために，我々は，メッシュを k × k（k = 2, ..., 6）と細かくした問題を作成し，解くのに必要な計算時間を計測した．以下，最大の問題 k = 6 の場合について結果を報告する．2 次錐の数は 60084，主双対変数の次元は 180252（= 60084 × 3）である．また，y の次元は 59850 である．我々は MPC アルゴリズムでこの問題を 725 秒，34 反復で解くことができた．また，基本アルゴリズムでは 1940 秒，110 反復を要した．
（停止条件は上述のものと同じである．
）我々の実装は，特殊な問題に対するものではあるが，計算速度に関しては，MOSEK（Anderson et al., 2003）や SeDumi（Sturm, 1999）など，世界第一線で用いられているものよりも高速であるか，少なくとも遜色ないものである．
4.2 結果の解析
本節では最適化結果について物理的立場から解析する．最適値，すなわち車体の四分の一のシールドの体積は，2.525069751 × 10−2 m3 であった．したがって，最適化されたシールドの体積は，1.0100279004 × 10−1 m3 (= 4 × 2.525069751 × 10−2 m3 ) となる．結果を図 5 に示す．図中
x, y, z-軸と厚み（左上の棒）の単位はメートルである．
この図を説明する前に，超電導磁石ユニットによって生成される磁場について説明する．2.2
節で述べたように，全領域は，車体部と廊下部から構成される．
（図 1 および 4 を参照のこと．
）
車体は，底板，側板，天井板，そして接合板（車体と廊下の接合部；専門的には妻と呼ばれる）からなり，廊下は底板，側板，および天井板からなっている．
車体に流入あるいはそこから流出する磁束の流れについて説明する．図 4（a）を見ると分かるように，磁束は底板の廊下側，そして接合板の下部から流入する．一方，図 4（b）から見て取れるように，側板の接合板に近い部分から流出している．
廊下部については以下のようになる．図 4（a）を見ると分かるように，側板の車体側下半分と底板の車体側が磁束の流入部である．一方，図 4（b）からも見て取れるように，車体から遠い方の側板の下半分と底板からは磁束が流出する．これが，超電導磁石ユニットにより生成されるシールド表面における Bn の様子である．

308

統計数理 第 53 巻

図 5.

第 2 号 2005

最適化されたシールド．図中左上の長方形は厚みを表すグレースケール（単位：メートル）．視点は x 軸無限大方向が右，z 軸無限大方向が上である．
（a）下左側から見た図；（b）左側から見た図；（c）下側から見た図．

2 次錐計画問題による磁気シールドのロバスト最適化

309

さて，以下では最適化されたシールド内部における磁束の流れを図 5 に基づいて説明しよう．
廊下部の底板では，磁束は進行方向に垂直（x-方向）に流れる．これは 2.1 節で説明したように，
廊下の両側に配置された超電導磁石ユニットが同じ方向に励磁されるためである．シールドは，
進行方向に沿ってみた時に，y = −7.8 と y = −9.5 の 2ヶ所で厚くなっている．これらの厚い部
分は磁束の流れが強い部分に対応する．これらの 2ヶ所の厚い部分で磁束の流れは互いに反対となる．
次に，車体部においては，磁束は底板から側板に，そして接合板から側板に流れる．これは，底板と接合板が超電導磁石ユニットからの磁束の流入部であり，磁束は側板からまた超電導磁石ユニットに戻っていくからである．このような磁束の流れを反映して，側板においてシール
ドは (x, y, z) ∼ (0, −7, 0.3) および (x, y, z) ∼ (0, −7, 0) の 2ヶ所で厚くなっている．しかしながら，
この部分だけでは流入する磁束を完全にバランスさせることができないため，接合板の廊下下部に x 方向に車体の反対側に向けて大変強い磁束の流れが存在する．そして，この部分がシールドの一番厚い部分（厚み約 1.5cm）となる．ここで興味深い点は，この部分については外部から流入したり外部に流出する磁束の流れはあまり強くないことである．外部とのやりとりはあ
まりないが，全体の流れのバランスをとるために強い磁束の流れが存在し，したがって厚くしておく必要があるというわけである．これは自明なことではなく，モデルによる計算推論ではじめて認識，確認される点であることを強調しておきたい．
5.

ロバスト最適化への応用

本節では主双対内点法の磁気シールドのロバスト最適化への応用について紹介する．まず，我々が扱っている問題をモデルに従って完全に最適化することが妥当であるかどうかを今一度検討してみよう．モデルは近似的なものであるから，ある程度の誤差や不確定な点が存在する．
したがって理想的には，ある固定されたモデルに対して最適なシールドを設計するよりは，モデルの不確定要素を考慮した上で現実の状況がモデルと多少異なってもシールドがうまく働くという条件の下で “最適化” することが望ましい，とも考えられる．このような考えに基づいて最適化を行うことをロバスト最適化といい，近年積極的に研究が進められている（Ben-Tal and
Nemirovski, 1998, 2001; Calaﬁore and Camp, 2005; El Ghaoui and Lebret, 1997; El Ghaoui etal., 1999; Goldfarb and Iyengar, 2003）．従前においても実際の設計において，モデルと現実の差は意識されており，モデルを元にする場合であっても十分な安全係数などをかけた上で設計が行われてきたことはいうまでもないが，ロバスト最適化は，この部分をある程度モデル化に含めようという試みであるとも考えられる．
ロバスト最適化においては，想定される最悪の場合を考慮した上での最適化を行うことが必要となる．与えられた最適化問題をこのように再定式化した問題をロバスト化した最適化問題という．ロバスト化した最適化問題は元の問題よりも難しいが，想定される摂動集合や元の最適化問題の性質によっては 2 次錐計画問題や半正定値計画問題になることもある（Ben-Tal and
Nemirovski, 2001）．
本節では，ロバスト磁気シールド最適化問題を考え，開発した主双対内点法を用いてこれを近似的に解き，さらに解のロバスト性能を最尤法を用いて定量的に評価することを試みる．以
下では 2.2 節と同じ記法を用いる．我々のモデルの不確定性は主として “シールド自身を有限要素で離散化したこと” と “外部磁場が近似解であること” に由来する．ここでは特に，外部磁場の不確定性に対するロバスト最適化を行う．外部磁場は近似的に計算されたものであるから，それがある程度想定と異なる場合でも，磁場が遮蔽できるように設計することが必要となる．
今，外部磁場の計算値が有限要素毎に ±5% の相対誤差を含むものとしよう．すなわち，外部310

統計数理 第 53 巻

第 2 号 2005

磁場が以下のような集合 Γ 内に存在するとする：

n∗ )j + (∆Bn∗ )j , |(∆Bn∗ )j | ≤ 0.05|(Bn∗ )j |, j ∈ E}.
Γ ≡ {Bn∗ ∈ R|E| |(Bn∗ )j = (Bn∗ は前節で最適化を行った際に採用したものである．そして，Gi (Bn∗ ) を Bn∗ と i 番目
ここで，Bの節点に対応する試験関数の内積とし，G(Bn∗ ) = (G1 (Bn∗ ), ..., G|V | (Bn∗ )) と定義する．すると，
n∗ は Bn∗ = G(Bn∗ ) と書ける（（2.9b）参照）．Bn∗ が Γ の中にあるとして変化し
2.2 節で定義した Bたときの（2.9b）の右辺集合を G(Γ) と記す．G(Γ) は R|V | 内の多面体である．ロバスト化した最適化問題を以下書き下す．
∗∗∗
以下，F0iを i 番目の要素のシールドの厚みとし，シールドを，F0∗ ≡ (F01
, ..., F0|E|) として表∗現する．これをシールド F0 と呼ぶことにする．
（以下シールドについて同様の記法を用いる．
）
与えられた外部磁場 Bn∗ に対し，シールドがそれを遮蔽できるだけ十分に厚いこと，すなわち，


(5.1)

j∗ = Gi (Bn∗ ),Dij F

i ∈ V,∗j∗  ≤ F0j
F
,

j ∈E

j∈E
∗
∗
∗
を満たすような F ∗ ≡ (F1∗ , ..., F|V| ) が存在する場合に，“シールド F は外部磁場 Bn を遮蔽で∗∗∗きる” ということにする．そして，F0 が任意の Bn ∈ Γ を遮蔽できる場合に，F0 を “（Γ について）ロバストなシールド” と呼ぶことにする．言葉を替えていうならば，Γ についてロバストなシールドとは，Γ の範囲で外部磁場が変化しても，磁場が車内に漏洩することがないだけの十分な厚みがあるシールドのことである．
さて，以上のような設定の下で，ロバスト磁気シールド最適化問題は以下のように書くことができる：

minimize
(5.2)

1 
∗
wj F0j,Bs j∈E∗∗, ..., F0|E|) ∈ Θ,
subject to F0∗ ≡ (F01∗|E|
| F0∗ は Γ に関してロバストなシールドである．
}．
Θ ≡ {F0 ∈ R

この最適化問題は半無限凸計画問題であり，これを厳密に解くことは困難である．そこで，
以下ではこの問題の近似解を求めることを考える．Γ から N 個の標本点 (Bn∗ )1 , ..., (Bn∗ )N を一様分布に従う乱数を用いて抽出し，Γ を ΓN ≡ {(Bn∗ )1 , ..., (Bn∗ )N } で置き換えた問題を解く．この近似問題は，
minimize
(5.3)

1 
∗
wj F0j,Bs j∈Esubject to F0∗ ∈ ΘFinite ,}ΘFinite ≡ {F0∗ ∈ R|E| | F0∗ は ΓN に関してロバストなシールドである．

と記すことができる．ΓN は有限集合であるから，（5.3）は以下のような |E| × N 個の 2 次錐制約を有する 2 次錐計画問題として定式化できる：
minimize(5.4)

subject to1 ∗
wj F0j,Bs j∈Ej∗ )k = Gi ((Bn∗ )k ),Dij (Fi ∈ V,

k = 1, ..., N,j∈E∗j∗ )k  ≤ F0j
(F,j ∈ E,

k = 1, ..., N.

Γ をその有限部分集合で近似する接近法を取る限りにおいては Γ についてロバストなシール
ドを設計するのは原理的に不可能である．一方，上で述べた ΓN についてロバストになるよう

2 次錐計画問題による磁気シールドのロバスト最適化

311

なシールドは ΓN に属する高々有限個の外部磁場を遮蔽できることが保証されているにすぎな
いが，実際にはそれで Γ に含まれる外部磁場のかなりの部分を遮蔽できる可能性もある．したがって，“設計されたシールドが Γ の中でどの程度の割合の外部磁場を遮蔽できるか” を知ることが，得られるロバスト最適化の近似解の評価において重要である．以下，シールド S が遮蔽できる Γ 内の外部磁場の集合を，Γ(S)と記すことにする．そして，vol(Γ) − vol(Γ(S))
vol(Γ)を “非許容度” と呼ぶことにする．ここで vol(·) は領域の体積である．Calaﬁore and Campi（2005）では（5.4）の最適解の非許容度を理論的に評価し，1 − β の確率で非許容度を ε 程度にするためには，N = O(1/(εβ)) とする必要があることを示している．これは理論的に興味深い結果であるが，今の文脈でこの結果を適用しようとすると，現実的な値として，ε = 0.01, β = 0.01ととった場合，必要な N は 10000 程度となり，解くべき 2 次錐計画問題は，1669 × 10000 個の2 次錐制約を含むため，実際に解くことは不可能である．
そこで以下では，主双対内点法の高速性と頑健性を生かした以下のような反復近似解法によって（5.3）の近似解を求め，さらにその非許容度を最尤法で統計的に評価する．
 が想定された時に，シーシールド S = (S1 , ..., S|E| ) が与えられているとする．“外部磁場 B
 を遮蔽できる一番軽いシールドを求める問題” は，以下のような 2
ルド S よりも各点で厚く B
次錐計画問題となる．
1 ∗
minimize
wj F0j,Bs j∈E(5.5)

subject toj∗ = Gi (B),Dij Fi ∈ V,j∈E∗,
Sj ≤ F0jj ∈ E,∗j∗  ≤ F0j,
Fj ∈ E.

 を遮蔽できる場合は，F0∗ = S がこの問題の最適解ともし，シールド S がすでにそのままで B
なる．
そこで，4 節で得られた最適シールドをシールドの初期形状として与え，以下のような考え方でシールドを逐次変形する．すなわち，ΓN の各要素を順番に並べておき，想定される外部磁場として逐次与え，（5.5）を主双対内点法で解くことにより，“その時点でのシールドよりも厚く新たな外部磁場を遮蔽できる最軽量のシールド形状” を求め，それを新しいシールドの形状とすることを繰り返す．このような手続きでシールドの形状を反復毎に少しずつ厚くしていく
ことをすべての ΓN の要素について行い，最終的に得られたシールドをロバスト最適化問題の近似解として採用する．一反復では自由度が 3000 程度の線形不等式制約を含む 2 次錐計画問
題を解くことが必要となる．この反復の終了時に得られているシールドは，全ての ΓN の要素に対して十分な厚みを持つので，（5.4）の許容解となっており，したがって，その体積は，（5.4）の最適値の上界となっている．前節で求めた最適シールドの体積は（5.4）の最適値の下界となっているため，両者を比較することにより，近似解としての良さも評価できることに注意したい．
この近似解法を実装し，反復を 10000 回繰り返した結果，10000 個の想定外部磁場に対して磁場の漏洩が起こらないシールドを元のシールドの 1.1% の重量の増加で設計することができた．計算は 4 節で用いたのと同じパソコンで行われ，要した実行時間は 6 時間程度であった．
図 6 に反復回数とシールドの体積の増加の様子を示す．
さて，このシールドの非許容度は以下のようにして推定することができる（笹川, 2005）．S k
を上述の手続きの k 反復目で得られたシールドとする．上述の手続きの中で，シールドの更新

312

統計数理 第 53 巻

図 6.

第 2 号 2005

反復回数とシールドの体積 (単位：10−2 m3 ) の関係．

は K = 129 回起こっている．S 129 が最終的に得られたシールドである．j 回目のシールドの更新が起こった反復を kj と記す．
pj ≡

(5.6)

 kj ))
vol(Γ) − vol(Γ(Svol(Γ)

 kj ) は S kj によって遮蔽できる外場の集合なので，最終段階での pj の値が非許とすると，Γ(S
容度を与える．我々の目的は，pj を合理的手段で推定することである．この問題は以下のような，2 項確率を推定する問題に帰着することができる．
外部磁場 ΓN は Γ 内で一様にランダムに生成されていることを考えると，kj 回目以降 kj+1
回目の反復までの間での反復で更新が起こる確率は（5.6）で与えられる pj であり，実際に反復が起こったのが kj+1 回目であったという事実は，2 項確率 pj のベルヌーイ試行を行った結果，
kj+1 − kj 回目に更新が起こったと見なすことができる．また，p1 ≥ p2 ≥ · · · ≥ p129 である．
これらの条件より p1 , ..., p129 をパラメータとする統計モデルおよびその尤度を書き下すことができ，これらのパラメータを最尤推定することができる．この最尤推定問題は，内点法によって多項式時間で解ける凸計画問題となり，高速にしかも厳密に解ける．この手続きを実装し計算を行った結果，最終的に得られた非許容度の推定値は約 0.003 となった．これは，シールドが遮蔽しきれない “たちの悪い” 外部磁場は，Γ の中で 0.3%に過ぎないことを意味する．さらに，尤度関数の形状を観察することにより推定値の妥当性についても考察することができる．
6.

まとめ

本論文では超電導磁気浮上式列車に対する磁気シールド最適設計問題を連続版の 2 次錐計画
問題として定式化し，それを有限要素法による離散化と 2 次錐計画問題に対する主双対内点法を用いて効率的に解くことができることを示した．さらに，主双対内点法をロバスト最適化に応用し，数千変数の 2 次錐計画問題を 10000 回解くことにより，軽量でロバストなシールドを設計できることを示し，そのロバスト性能について，統計的評価を行った．ロバスト最適化に

2 次錐計画問題による磁気シールドのロバスト最適化

313

際しては外部磁場の不確定性のみを論じ，元の問題を有限要素近似することによる離散化の影
響については論じなかったが，本論文のものより精密な有限要素を用いて離散化問題を構成して双対理論を用いることにより，元の連続問題の最適値の上下限をある程度見積もることができる（笹川, 2005）．
謝辞本研究の初期に，いろいろとご助言をいただいた野末尚次氏（（株）数理モデリング）に感謝いたします．また，図の作成にご協力いただいた堀内信吾氏（NTT）にも感謝いたします．最後になりましたが，著者らの研究生活の駆け出しの時期より，内点法をはじめとして最適化，数値
線形代数などについて数々の指導をいただいた田邉國士統計数理研究所名誉教授に深く感謝いたします．

参 考 文 献Alizadeh, Farid（1995）. Interior point methods in semideﬁnite programming with applications to
combinatorial optimization, SIAM Journal on Optimization, 5, 13–51.
Alizadeh, F., Haeberly, J.-P. and Overton, M.（1998）. Primal-dual interior-point methods for semidefinite programming: Convergence rates, stability and numerical results, SIAM Journal on Optimization, 8, 746–768.
Andersen, E. D., Roos, C. and Terlaky, T.（2003）. On implementing a primal-dual interior-point
method for conic quadratic optimization, Mathematical Programming, 95, 249–277.
Ben-Tal, A. and Nemirovski, A.（1998）. Robust convex optimization, Mathematics of Operations
Research, 23, 769–805.
Ben-Tal, A. and Nemirovski, A.（2001）. Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications, SIAM Publications, Philadelphia, Pennsylvania.
Calaﬁore, G. and Campi, M. C.（2005）. Uncertain convex programs: Randomized solutions and con（1）, 25–46.
ﬁdence levels, Mathematical Programming, 102
Ciarlet, P. G.（1978）. The Finite Element Methods for Elliptic Problems, North-Holland PublishingCo., Amsterdam.
El Ghaoui, Laurent and Lebret, Herve（1997）. Robust solutions to least-squares problems with uncertain data, SIAM Journal on Matrix Analysis and Applications, 18, 1035–1064.
El Ghaoui, Laurent, Oustry, Francois and Lebret, Herve（1999）. Robust solutions to uncertain semidefinite programs, SIAM Journal on Optimization, 9, 33–52.
Faraut, J. and Korányi, A.（1994）. Analysis on Symmetric Cones, Oxford University Press, NewYork.
Faybusovich, L.（1995）. Jordan algebras, symmetric cones and interior-point methods (manuscript).
Faybusovich, L.（1997）. Linear systems in Jordan algebras and primal-dual interior-point algorithms,
Journal of Computational and Applied Mathematics, 86, 149–175.
Goldfarb, D. and Iyengar, G.（2003）. Robust portfolio selection problems, Mathematics of Operations（1），1–38.
Research, 28
Helmberg, C., Rendl, F., Vanderbei, R. J. and Wolkowicz, H.（1996）. An interior-point method forsemideﬁnite programming, SIAM Journal on Optimization, 6, 342–361.
Kojima, M., Mizuno, S. and Yoshise, A.（1989）. A primal-dual interior point algorithm for linear
programming, Progress in Mathematical Programming, Interior-Point and Related Methods（ed. N. Megiddo）, 29–47, Springer-Verlag, New York.

314

統計数理 第 53 巻

第 2 号 2005

Kojima, M., Shindoh, S. and Hara, S.（1997）. Interior-point methods for the monotone semideﬁnite
linear complementarity problem in symmetric matrices, SIAM Journal on Optimization, 7,
86–125.
Lobo, M. S., Vandenberghe, L., Boyd, S. and Lebret, H.（1998）. Applications of second-order cone
programming, Linear Algebra and Its Applications, 284, 193–228.
Mehrotra, S.（1992）. On the implementation of a primal-dual interior point method, SIAM Journal
on Optimization, 2, 575–601.
Monteiro, R. D. C.（1997）．Primal-dual path-following algorithms for semideﬁnite programming, SIAM
（3），663–678.
Journal on Optimization, 7Monteiro, R. D. C. and Tsuchiya, T.（2000）. Polynomial convergence of primal-dual algorithms for thesecond-order cone program based on the MZ-Family of directions, Mathematical Programming,
88, 61–83.
Nesterov, Yu and Nemirovskii, A.（1994）. Interior Point Polynomial Methods in Convex Programming, SIAM Publications, Philadelphia, Pennsylvania.
Nesterov, Y. E. and Todd, M.（1997）. Self-scaled barriers and interior-point methods for convex
programming, Mathematics of Operations Research, 22, 1–42.
笹川 卓（2005）. 2 次錐計画問題を用いた直流磁気シールドの最適化，博士論文，東京大学大学院情報理工学系研究科．
Sasakawa, T. and Tagawa, N.（1997）. Preliminary design of magnetic shielding by FEM, IEEE Transactions on Magnetics, 33（2）, 1951–1954.
Sasakawa, T. and Tsuchiya, T.（2000）. Optimal magnetic shield design with second-order cone programming, Research Memorandum, No. 775, The Institute of Statistical Mathematics, Tokyo.
Sasakawa, T. and Tsuchiya, T.（2003）. Optimal magnetic shield design with second-order cone programming, SIAM Journal on Scientific Computing, 24（6）, 1930–1950.
佐々木寛（1993）．磁気遮蔽，川崎製鉄技報，25（1）, 64–70.
Schmieta, S. and Alizadeh, F.（2001）. Associative algebras, symmetric cones and polynomial time（3）, 543–564.
interior point algorithms, Mathematics of Operations Research, 26
Sturm, J. F.（1999）. Using SeDuMi 1.02, A matlab toolbox for optimization over symmetric cones,
Optimization Method and Software, 11&12, 625–653.
Taguchi, A. and Iri, M.（1982）. Continuum approximation to dense networks and its application to
the analysis of urban road networks, Mathematical Programming Study, 20, 178–217.
Tanabe, K.（1988）. Centered Newton method for mathematical programming, System Modeling andOptimization（eds. M. Iri and K. Yajima）, 197–206, Springer-Verlag, Berlin.
鉄道技術総合研究所（2005）. Railway Technical Research Institute, http://www.rtri.or.jp/
Todd, M. J., Toh, K. C. and Tütüncü, R. H.（1998）. On the Nesterov-Todd direction in semideﬁnite
programming, SIAM Journal on Optimization, 8, 769–796.
Tsuchiya, T.（1999）. A convergence analysis of the scaling-invariant primal-dual path-following algorithms for second-order cone programming, Optimization Methods and Software, 11&12,
141–182.
土谷 隆（2005）. 層別最小二乗法―重み付き最小二乗法の極限―，統計数理，53（2）, 391–404.
（2000）. Handbook of Semidefinite ProgramWolkowicz, H., Saigal, R. and Vandenberghe, L.（eds.）
ming. Theory, Algorithms, and Applications, International Series in Operations Research &Management Science, 27, Kluwer Academic Publishers, Boston, Massachusetts.
Wright, S.（1997）. Primal–Dual Interior Point Algorithms, SIAM Publications, Philadelphia, Pennsylvania.

Proceedings of the Institute of Statistical Mathematics Vol. 53, No. 2, 297–315 (2005)

315

Robust Optimization of Magnetic Shielding with
Second-order Cone Programming
Takashi Tsuchiya1 and Takashi Sasakawa2
1 The Institute of Statistical Mathematics
2 Railway Technical Research InstituteThis paper deals with an application of second-order cone programming to an optimal magnetic shielding design problem of the MAGLEV (superconducting MAGnetically
LEVitated vehicle), a new bullet train under development in Japan. The MAGLEV is heldin the air by a strong magnetic ﬁeld and propelled by linear synchronous motors. Eachcar is equipped with several super-conducting magnet units which generate the magnetic
ﬁeld. Passengers inside the car need to be shielded from the magnetic ﬁeld outside. The
optimal design problem of magnetic shielding is to minimize the weight of the shielding
by adjusting its thickness taking into account the external magnetic ﬁeld.
After some appropriate simpliﬁcation, this optimization problem is formulated as asecond-order cone program and is solved with the primal-dual interior-point algorithm.
Taking advantage of its eﬃciency and stability, the algorithm is further applied to robust
design of magnetic shielding.

Key words: Second-order cone programming, interior-point methods, magnetic shielding, MAGLEV
train, optimal design, robust optimization.