計算応用科学分野 - 名古屋工業大学 創成シミュレーション工学専攻

計算応用科学分野
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
有限要素法による非一様断面チモシェンコ梁の動的解析
学籍番号 22413504
氏名 有本 博和
指導教員名
1 はじめに
杉山 勝
一方、断面に一様なせん断変形が生じると仮定す
機械や構造物を安全に設計するためには、構造部
材の振動特性を理解することが重要となる。本研究
るのがチモシェンコ梁理論である。チモシェンコ梁
の曲げ振動方程式は次のようになる。
では梁の曲げ振動を例にとり、そのモデル化の違い
によって振動特性がどのように変化するか観察する。
梁はいわゆる連続体であり、
空間に3次元的拡がり
を持っているが、ある1方向への拡がりが他の2方向
への拡がりよりもかなり大きいため、支配方程式を
簡略化できる。これが梁理論の出発点である。梁の
ここで はせん断弾性係数、
はチモシェンコ定数、
は断面の回転角である[3]。
モデルとして代表的なものに、ベルヌーイ・オイラ
ーの仮定を満足するベルヌーイ・オイラー梁と断面
3 両モデルの違い
に一様なせん断変形が発生すると仮定するチモシ
3.1 細長比による違い
ェンコ梁がある。梁の境界条件は左端固定の片持ち
図 1 は横軸に細長比
、縦軸にベル
ヌーイ・オイラー梁とチモシェンコ梁の各固有振動
梁とした。
一般に、
梁の曲げ振動方程式は4階の偏微分方程式
となり、特別な場合を除いて解析解を求めることは
困難である。そこで有限要素法（FEM）を用いてこれ
数の相 対差
をとった ものであ る。 ここで
、 は梁の全長さ、
を数値的に計算した。プログラム言語としては
はそれぞれ、ベルヌーイ・オイラー梁とチモシェン
MATLABを用いた[1-3]。
コ梁の第 次固有振動数である。
2 支配方程式
ベルヌーイ・オイラー梁理論では、変形の前後で
梁の断面形が変化せず、変形前に梁の中立軸に直交
していた任意断面は変形後も中立軸に直交すると仮
定する。
（これをベルヌーイ・オイラーの仮定と呼ぶ）
。
この理論では梁のせん断変形を無視し、梁は中立軸
と直角方向の並進運動のみを行うとしている。ベル
ヌーイ・オイラー梁の曲げ振動方程式は次のように
なる。
図1：細長比と固有振動数
ここで
は梁のたわみ量、 は単位体積の質量
密度、 は梁の断面積、 は梁のヤング率、 は梁の
断面 2 次モーメント、
図 1 を見ると、細長比が大きくなると固有振動数
は外部から加えられた
の相対差が小さくなることがわかる。すなわち、断
圧力荷重、 と は時間と梁の軸にそった空間の座標
面に比べて長い梁では、ベルヌーイ・オイラー梁と
である [3] 。
チモシェンコ梁の差が小さくなっていく。逆に短い
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
梁では両モデルの差が大きくなる。また、高次の固
計算応用科学分野
ったものである。図3を見ると、振動しながら移動す
有振動数は低次の固有振動数に比べ両モデルの差が
る集中外力荷重が自由端に接近すると、自由端の振
大きくなることもわかる。
動数が変化していることがわかる。すなわち自由端
の振動が自由振動から強制振動に移行していること
3.2 テーパー度による違い
テーパー度を
がわかる。また、ベルヌーイ・オイラー梁の方が強
と定義する。
制振動に早く移行していることもわかる。
はそれぞれ左端と右端の断面高さである。図 2 は横
軸に梁のテーパー度
、縦軸に振動数をとったも
のである。
図3：強制振動による梁の振動数の変化
4 まとめ
本研究では、梁の形状や外力条件を変更すること
で現れるベルヌーイ・オイラー梁とチモシェンコ梁
図2：テーパー度と固有振動数
の違いを有限要素法を用いて定量的に観察した。そ
の結果、形状（細長比、テーパー度）変化に伴う第1
図 2 を見るとテーパー度を変化させても第 1 次固
次固有振動数の変化は両モデル間で大きな違いが出
有振動数は両モデルでほとんど違いが現れないこと
ないことを示した。一方で高次の固有振動数は形状
がわかる。一方でベルヌーイ・オイラー梁の高次の
に大きく依存し、ベルヌーイ・オイラー梁の高次の
固有振動数はチモシェンコ梁に比べ非常に大きな値
固有振動数は、チモシェンコ梁のそれより常に大き
をとっている。また、両モデル共にテーパー度が大
くなることも確認した。
きくなると低次の固有振動数は上昇し、高次の固有
振動数は減少していることも観察できる。
また、固定端から自由端に向けて周期的に振動し
ながら移動する集中外力荷重に対して、ベルヌー
イ・オイラー梁では、集中荷重の接近に伴い、自由
3.3 外力の作用位置による違い
周期的に振動する外力が作用すると、梁は強制振
動を起こす。強制振動中の梁の振動数は強制外力の
端が急速に強制振動状態に移行するが、チモシェン
コ梁では比較的ゆっくりと強制振動状態に移行する
ことがわかった。
振動数に近づき、外力の作用がなくなると梁は自由
振動に移行する。
テーパー梁の上を固定端から自由端に向けて、周
5 参考文献
[1] Jacob Fish・Ted Belytschko:有限要素法 ABAQUS Student
期的に振動しながら、一定速度で移動する集中外力
Edition 付,丸善出版（2008）
荷重がある系を考える。ここで、この集中外力荷重
[2] 小松敬治：構造振動学 MATLAB による有限要素法と応
の振動数は両モデル共に共振を起こさない振動数を
答解析,森北出版株式会社（2008）
とる。
このときの自由端の振動数の変化を図3に示し
[3] Young W,Kwon Hyochoong Bang :The Finite Element
た。図3は、横軸に時間、縦軸に自由端の振動数をと
Method using MATLAB,CRC Press(2000)
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
サポートベクトルマシンを用いた多変量 2 標本検定とパス追跡による高速化
学籍番号 22413506 氏名 磯部 浩太
指導教員名 竹内一郎 准教授
1
まえがき
本稿では, サポートベクトルマシン (SVM) に基づく
多変量 2 標本検定について検討する.
多変量 2 標本検定とは, 多変量分布が 2 標本で異なっ
ているかどうかを統計的に評価することである. 2 標本
での分布の違いを評価するのに 2 クラス分類器から計
算される何らかの検定統計量を用い, その統計的信頼
性を評価することで多変量 2 標本検定を行うことがで
きる. 非線形やノンパラメトリックな 2 クラス分類器
を利用すれば, 複雑に分かれている 2 標本の違いを検
出できる利点がある. 本稿では, 2 クラス分類器として
SVM を用い, 学習後の SVM から計算される検定統計
量を用いて多変量 2 標本検定を行う問題を検討する.
統計的検定においては検定統計量の帰無分布を知る
必要がある. 帰無分布とは, 2 標本が同一の分布から生
じている場合の検定統計量の分布である. 2 クラス分類
器の検定統計量の帰無分布は一般に未知であり, ラベル
並べ替え演算などのランダムシミュレーションを用い
て推定しなくてはならない. ラベル並べ替え演算によ
り帰無分布を推定するアプローチでは, 分類器 (SVM)
の学習 (最適化) をラベル並べ替えの回数分だけ行わな
くてはならず, 計算コストが非常に大きくなってしま
う. 本稿では, これを効率的に計算するため, パス追跡
[1] を用いたラベル並べ替えサンプルの効率的計算方法
を提案する. 以下では, n 次元縦ベクトルを v ∈ Rn の
ように表記する. また, Nn は 1 から n までの自然数の
集合 {1, 2, . . . , n} を表すものとする.
2
ラベル並べ替え検定
統計的検定において 2 標本の差の有意性を定量化す
る (p 値などを求める) には, 検定統計量の帰無分布を
知る必要がある. 本節ではラベル並べ替え演算を利用
して帰無分布を推定するラベル並べ替え検定を説明す
る. y = {y1 , . . . , yn } を観測データのラベルの n 次元
ベクトルとする. ラベル並べ替え検定では, y をランダ
ムに並べ替え, 並べ替えたラベルに対して検定統計量を
求めることを繰り返すことで, 検定統計量の帰無分布を
推定する. ラベルをランダムに並べ替えることによっ
て仮想的な差のない 2 つの分布 (帰無仮説) を作成して
いるものと解釈することができる. 実用的には, ランダ
ムなラベル並べ替え演算を 1000∼10000 回繰り返して
検定統計量の帰無分布を推定し, 得られた帰無分布と正
しいラベルの基での統計量を比較することで統計的信
頼性を評価することができる.
SVM を用いた多変量 2 標本検定
3
SVM の 分 類 性 能 を 検 定 統 計 量 と し た 多 変 量 2
標 本 検 定 に つ い て の 説 明 を 行 う.
学習データを
{(xi , yi )}i∈Nn , xi ∈ Rp , yi ∈ {−1, +1} とすると, SVM
の学習は以下のような二次計画問題として定式化される:
max
{αi }i∈Nn
s.t.
∑
1 ∑ ∑
αi αj yi yj K(xi , xj ) +
αi ,
2
i∈Nn j∈Nn
i∈Nn
∑
αi yi = 0, 0 ≤ αi ≤ C.
−
i∈Nn
ここで, C ∈ [0, ∞) は正則化パラメータ, {αi }i∈Nn は
ラグランジュ未定乗数, K(xi , xj ) はカーネル関数を表
している. また, 分類境界は,
∑
f (x) =
αi yi K(xi , xj ) + b = 0,
(1)
i∈Nn
と与えられる. 2 標本が大きく離れている場合, 各データ
点と分類境界 (1) の距離が大きくなると思われる (デー
タ点 (xi , yi ) の分類境界までの距離は yi f (xi ) と計算さ
れる). 本稿では, 評価サンプルに対する SVM 分類境
界までの平均距離を検定統計量とする. Leave-one-out
交差確認により評価サンプルを作成するとし, データ i
を除いて学習された SVM 分類境界を f (−i) とすると,
本稿で採用する検定統計量は
∑
s := n−1
yi f (−i) (xi ),
(2)
i∈Nn
と定義される. 本稿では検定統計量 (2) を用いたラベル
並べ替え演算による多変量 2 標本検定を提案する.
3.1
計算機実験
実験により, 提案法の有効性について検証する. 実験
は, 実データからリサンプリングを 100 回行い, 有意水
準は 0.05 とした時の帰無仮説を棄却した割合を求め, 検
出力と頑強性を調べる. リサンプリングは 2 標本が同一
分布となるようにリサンプリングした場合 (Same) と,
2 標本が異なった分布となるようにリサンプリングした
場合 (Different) の 2 通りの方法を用いる. Same では棄
却した割合が有意水準に近ければ頑強性が高いことを示
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
し, Different では棄却した割合が高ければ検出力が高い
ことを示している. 用いる実データは UCI データベー
スから引用した BreastCancerPrognostic, Parkinsons,
を用いる. 比較する手法は, SVM 検定統計量 (2) を用い
た検定 (提案法), 最近傍検定 (NN), MMD[2], LSTT[3],
とする. 実験結果の例を表 1, 2 に示す. 提案法の検出
力は NN に劣っているが, LSTT, MMD より高いこと
がわかる. また, 頑強性については提案法の方が NN よ
り高い. このことから, 提案法が有効な手法であること
がわかる.
表 1: BreastCancerPrognostic
提案法
NN
MMD LSTT
3%
100%
24%
100%
5%
100%
表 2: Parkinsons
提案法 NN MMD
Same
Different
4
4%
51%
24%
78%
3%
28%
6%
100%
• 階層型クラスタリングに基づくスケジューリング,
を提案する.
計算機実験により, パス追跡による計算コスト削減を
検証する. 比較する手法は SMO アルゴリズム (SMO),
最小全域木に基づくスケジューリングを行ったパス追
跡 (MST), 階層型クラスタリングに基づくスケジュー
リングを行ったパス追跡 (HC) とする. 実験結果の例を
図 4.1 に示す. 結果から提案法により計算コストが大幅
に削減されていることが分かる.
LSTT
6%
24%
SVM ラベル並べ替え解の計算の
効率化
ラベル並べ替え検定を行う際に, ラベル並べ替え回数
は 1000∼10000 回と多い. 従って, これらのラベル群そ
れぞれに対して SVM の学習を行うと, 計算コストが膨
大になる. そこで, 本研究では以下のようなアプローチ
で計算コスト削減を図る.
1. パス追跡による並べ替えラベル群の最適解の追跡
2. パス追跡スケジューリング
パス追跡とは, 学習データ点の追加, 削除や正則化係
数の変化などが起こった場合, 最初から学習しなおすの
ではなく, 最適解の感度分析に基づいて変更後の最適解
を効率的に計算する方法である. ラベル並べ替え検定
において並べ替えラベル群には共通する部分が多いた
め, 最適解の変化が少ないと考えられる. そのため, パ
ス追跡を用いて効率的にラベル並べ替えデータの学習
を行うことができる.
多数の並べ替えラベル群に対してパス追跡をを行う
場合, パス追跡の経路が大きな影響を与える. 例えば,
ある SVM の解から別の SVM の解へパス追跡を行う時
には, 互いの解が近ければ高速に計算ができるが, 遠け
れば時間がかかってしまう. このため, すべての並べ替
えラベル群に対してパス追跡を行う際に, 最適解の総変
化量ができるだけ小さくなるようなスケジューリング
を行うことが重要である. 本稿では方針として,
計算機実験
4.1
cputime(sec)
Same
Different
• 最小全域木に基づくスケジューリング,
1000
900
800
700
600
500
400
300
200
100
0
C=1
C=10
C=100
図 1: 計算時間の比較 (秒) 赤:HC 緑:MST 青:SMO
5
まとめ
本研究では, SVM を用いた多変量 2 標本検定とその
高速化を行った. そして, 実験によりその有効性を示す
ことができた. 今後の課題として, 様々な実問題への応
用などが挙げられる.
参考文献
[1] M. Karasuyama, N. Harada, M. Sugiyama, and
I. Takeuchi. Multi-parametric Solution-path Algorithm for Instance-weghted Support Vector Machines. arXiv:1009.4791, 2010.
[2] A. Gretton, K. M. Borgwardt, and M. Rasch et.al.
A kernel method for the two-sample-problem. Advances in neural information processing systems,
19, 513-520, 2007.
[3] M. Sugiyama, T. Suzuki, Y. Itoh, T. Kanamori
and M. Kimura. Least-Squares Two-Sample Test.
Neural Networks, 2011.
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
Wolfram クラス II に属する
0- 静的 1 次元基本セルオートマトンのカオス性について
学籍番号
22413507 市川 貴紀
指導教員 大鋳 史男
1 はじめに
ス的振る舞いに落ち着くクラス III，そして規則的
セルオートマトン(以下 CA)とは，時間，空間，状
なパターンとカオス的パターンが混在するクラス
態は離散的な数理モデルであり，計算機との相性
IV と分類する．
が良いという特徴を持つ．現在では広範囲におけ
クラス IV は「カオスの縁」と呼ばれる重要なクラ
る複雑現象解明のための有効な手法の一つとなっ
スであるが，この分類は見た目による判断であり
ている．離散化された空間をセル，離散化された
厳密性にかける．クラス I や II であってもセルの
時間をステップと呼ぶ．すべての各セルの状態が
状態配置を特定のパターンに限定することで，複
与えられたものを configration と呼ぶ．
雑な時間発展をすることがある．ルール 56 はクラ
1 次元 CA とは以下の条件を満たす(S,r,g)である．
ス II とされていたが，ある種の interactive な複
(1) S={0,1,…,n}
雑な動きを見せる場合が存在する．
(2) r は正の整数
ルール 40 は，クラス I に属するが，特定の
(3) g:S 2r＋１ → S
configuration に対する力学系は，σL になり，従
本論で取り扱う 1 次元基本セルオートマトン(以下
って，Devaney カオスになることが示されている．
ECA と略す)とは，セルの配置が一次元であり，
[2]を参照.ルール 168 の基本セルオートマトンも
近傍数 1，セルの取り得る状態は 0 か 1 の二種類
またクラス I に属するが，これについても同様の
のみというもっとも基本的な CA である．
ことが観察され，σL，σR の部分力学系を持つこ
ECA はセル自体と隣り合ったセルを含めた三つ
とが示されている．[4]を参照．一つの CA は、た
のセルの状態によって次のセルの状態が決まる．
だ一つの力学特性を持つのではなく，いくつかの
すなわち，
特性を保持し，初期状態の違いによってその中の
g(1,1,1),g(1,1,0),g(1,0,1),g(1,0,0),
一つの特性が出現する．このような事態はより高
g(0,1,1),g(0,1,0),g(0,0,1),g(0,0,0)
いクラスではさらに複雑だと予想される．また CA
以上 8 つの対応関係(ローカルルール)によって決
全体において，より基盤的な CA が存在し，それ
まる．g(0,0,0)=0 であるとき 0-静的であるという．
らの組み合わせによって任意の CA を構成できる
ルール番号を∑ g(a, b, c) ∗ 2
*4a+2b+c+
と定める．左
と推測される．
シフトをσL，右シフトをσR とする.1 次元 CA の
本論ではこうした分解構成問題に見通しを立てる
分類にはランダムに与えた初期状態の時間発展の
道筋をつけるべく，従来クラス II とされたルール
様子によって 4 つに分類する Wolfram クラスとい
20 とルール 14 について複雑な時間発展を起こし
うものがある．[1]を参照．有限ステップですべて
うるパターンについて言及する．
のセルの状態が同じとなる死滅とも呼ばれるクラ
2 ルール 20 について
ス I，有限ステップで周期的，あるいは規則的な変
ルール 20 についてランダムな初期 configration
化に落ち着くクラス II，予測不可能な複雑なカオ
を与えたときの時間発展は究極的にg220 = σ2R とな
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
る．
やルール 14 の解析手法が適応できる可能性が示
しかし 10,1100 の二種類のブロックのみの組み合
唆される．
わせによる configration の集合 SI 上ではg220 = σ2L
今後の課題としてそのほかのクラス II の ECA で
となる．
独自の部分力学系がないかを解明し，それらがも
g220 = σ2L となる力学を構成しているg2 20 の対応関
しあればどの対応関係によって生じているかを見
σ2R となる力学を構成しているg2 20 の対
つけ出し，より基礎的な CA に分解し，共通して
応関係を調べ，
内包している近傍数 2 の一次元 CA
いる CA や独自の CA の発見がある．そして，そ
の対応を抽出した．
れらの基礎的 CA からクラス II の各 ECA を構成
このσL となるパターンを右側に，σR となるパタ
するための理論の構成がある.
ーンを左側に配置し時間発展させることにより，
クラスⅢについては configration の構成次第でク
その境界(FrontLine)が左や右に interactive に動
ラス II の規則的な時間発展をする可能性が残され
き，また FrontLine の移動量を定式化することが
ている.
できた．これによりカオスの条件である稠密な周
それらの部分力学系と対応関係の発見と,より基
期的 configration，推移的 configration を構成す
礎的な CA への分解と各ルールの構成理論の確立
ることができた．また 10,1000 の二種類のブロッ
がその先にある．
クの組み合わせ configration の集合 SII 上でもσ2L
ECA 以外の一次元 CA,さらには二次元 CA，多次
となり，周期的，推移的 configration を構成する
元 CA に対する構成理論の基礎を確立することが
ことができた．よってg20 は SI から SII への写像が
その先の課題として残されている．
全単射であることが確認でき，それによって
6 参考文献
(SI＼SII，g20)での Devaney カオス性を示せた．
[1] S. Wolfram: Statistical mechanics of cellular
3
automata, Review of Modern
係とg220
=
ルール 14 について
ランダムな初期 configration を与えたときの時間
Physics, 55 (1983), 601-644.
発展は究極的にσL となる．10,100,110,1100 ブロ
[2]
ックにより構成された configration に対しては反
Exponents and Spreading Rates for Rule 40, Complex
転右シフトとなることを確認できた．
Systems,Vol.19(2011),No.4, pp.323-341.
ルール 14 に関しては 1step で，interactive な力
[3]Fumio Ohi, Dynamical Properties of Rule 56
学特性を持っており，2,3steps を検証してみたが，
Elementary Cellular Automaton of Wolfram Class II, in
新規の力学を確認することはできなかった．
Automata-2010,
2steps での稠密となる周期的 configration と推移
Cellular Automata and Discrete Complex Systems, edt.
的 configration を構成することができた．
Nazim Fates, Jarkko Kari and Thomas Worsch, 2010,
4
pp.287-298.
まとめと今後の課題
Fumio
Ohi,
Exact
16th
Calculation
International
of
Lyapunov
Workshop
on
そのほかの Wolfram クラス II の ECA にもランダ
[4] Fumio Ohi, Chaotic properties of elementary cellular
ム初期 configration を与えたときの
automaton rule 168 of Wolfram class I,
部分的な時間発展に interactive な複雑な挙動が
in AUTOMATA-2008, Theory and Applications of
見られるものがあった．
Cullular Automata,
たとえばルール 58，74，134，142 にそうした
edt. A. Adamatzky, R. Alsonso-Sanz, A. Lawniczak, G.
interactive らしき動きを確認できた．
Juarez Matrinez, K. Morita, T. Worsch,
これらの ECA に対しても,本論で行ったルール 20
Luniver Press, 2008, pp.196-206.
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用分野
セルオートマトン法を用いた複雑流動のシミュレーション
学籍番号
22413509
氏名
伊藤
雄一
指導教員 大鋳
史男
第Ⅰ部 セルオートマトン法を用いた商用
することで歩行者の挙動を実現している．
高層建築物における避難歩行者流動のシミ
歩行者は一回の行動で中心セルの周り 8 セ
ュレーション ―最適な避難計画構築の方
ルのいずれかに移動することができる．歩
針について―
行者が次の移動先セル決定する際には以下
1 目的
の手順を踏む．
名古屋市内に実在する百貨店を
対象とした避難シミュレーションを行い，
（1）現ステップで移動を試みるかの判定
最適な避難方針を確定するための基本方針
（2）移動先候補セルである各セルについて
を探り出す事を目的とする．
いくつかの評価要素を計算し，各セルの評
2 セルオートマトン法 セルオートマト
価点を得る．ムーアー近傍セルのなかで最
ン法とは時間，空間ともに離散的なシステ
も評価点が高いセルを移動先セルとして選
ムである．空間はセルと呼ばれる格子に均
択する．
一に分割される．セルは内部状態を持ち，
（3）
（2）で選択した移動先候補セルに移動
近傍のセルの状態とローカルルールと呼ば
可能か判定する．
れるルールに従って新しい状態を決定する．
本研究では評価要素として目標方向，障害
3 モデル化
物の影響，慣性の 3 つを用いている．
3．1 建物のモデル化 建物をセルオート
4
マトン法を用いてモデル化する．建物を
が短縮される最適な避難指示を考えるため，
10cm×10cm のセルに分割し，分割したセル
4 種類の避難指示を与えた場合のシミュレ
に空間の情報を持たせる．
ーションを行った．その結果，各避難指示
ノードとは歩行者の最終的な目的地までの
を与えた場合を比較すると，退避において
中間目的地の役割を果たすもので，ノード
直近の階段に向かうことは必ずしも適切で
を出口までつなげることで歩行者の退避経
はないこと，全館退避完了時間が最も短縮
路が構成される．またノードから構成され
される避難指示は全ての階段の退避完了時
るグラフにダイクストラ法を適用させるこ
間が均一になる傾向があることが分かった．
とで出口までの最短経路が求まる．本研究
建物の構造が違う場合でも全館退避完了時
ではノードを出口，扉，標識，障害物の角
間が短縮される最適な避難指示を考えるた
などに設置をして建物空間のグラフ構造を
め建物の構造を変更した場合に4種類の避
構成する．
難指示を与えた場合のシミュレーションを
3．2 歩行者のモデル化 歩行者を 3×3 セ
行った．その結果，全館退避完了時間が最
ルと表現し，その中心セルを歩行者の中心
も短縮される避難指示は全ての階段の退避
セルとする．この歩行者の中心セルが移動
完了時間が均一になる傾向がある．このこ
シミュレーション 全館退避完了時間
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用分野
とから全ての階段の退避完了時間が均一に
層が現れたがそのままルールの場合は現れ
なるように指示すると全館退避完了時間が
ない．この複数の層を解析した結果，第一
短くなると考えられ，この指示が最適な避
車線の渋滞，第二車線の渋滞，両車線の渋
難指示であると思われる．
滞などの渋滞の種類によるものだと分かっ
5 まとめ
シミュレーション結果から，避
た．そのままルールは車線変更しやすいル
難計画を立案する際のポイントがいくつか
ールなので 1 車線だけが渋滞になっても車
わかった．
線変更することで渋滞が解消するか両車線
・各フロアでの退避においては直近の階段
が渋滞となるが，一旦ルールでは車線変更
に向かうことは必ずしも適切ではないこと．
しづらいルールなので第一車線の渋滞，第
・全ての階段の退避完了時間が均一になる
二車線の渋滞，両車線の渋滞などになり，
ようにする指示が最適な避難指示であると
それぞれ流量の値が異なるためメタ安定分
思われること．
岐で複数の層が現れることが分かった．
第Ⅱ部 セルオートマトン法を用いた交通
9
流シミュレーション
い交通流モデルを提案し，1 車線と 2 車線
6 目的
本研究ではセルオートマトン法
での車のパラメータの値を変化させたとき
を用いた交通流のモデルを提案し，1 車線
の交通流及ばす影響について調べた．その
と 2 車線それぞれに於いて，車のパラメー
結果，次のことが分かった．
タの値を変化させたときの交通流の変化を
• 1車線モデル，2車線モデルともに時間遅
考察した．また，シミュレーションによっ
れまたはクイックスタートがメタ安定状態
てメタ安定分岐現象を分析することで自由
の発生要因である．従来は時間遅れがメタ
走行相から渋滞相に至るメカニズムについ
安定状態の発生の要因であると考えられて
て解析することを目的とした．
きたが，時間遅れが無くても，クイックス
7 本研究の交通流モデル 本研究では最
タートを組み込むことによりメタ安定分岐
高速度，最大加速度，最大減速度，時間遅
が発生した．
れ，クイックスタートを組み込んだ決定論
・一旦モデルではメタ安定分岐領域で2次元
的な 1 車線交通流モデルと 2 車線交通流モ
的な拡がりを示した. そのままモデルでは
デルを提案する．2 車線交通流モデルでは
メタ安定分岐は現れるがメタ安定分岐領域
車線変更のルールはそのままルールと一旦
での2次元的な拡がりは再現できない．
ルールを提案する．
10
8 シミュレーション
まとめ
本研究では確率要素を含まな
参考文献
1 車線モデルと 2 車
[1] 伊藤義弥，セルオートマトン法を用いた高層
線モデルの車のパラメータの値を変化させ
建築物における避難時歩行者流動のシミュレーシ
たときの交通流の変化を観察した．メタ安
ョン―最適な避難計画構築のための指針について
定分岐は時間遅れまたはクイックスタート
―，名古屋工業大学修士論文，2011
を組み込むことにより発生した．
[2] 西成活裕，交通流のセルオートマトンモデル
2 車線モデルにおいて車線変更のルールが
について，応用数理，特集・交通流の数理，Vol.12，
一旦ルールの場合はメタ安定分岐に複数の
No2，pp.26-37，2002
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
拡張された熱力学に基づく流体力学的ゆらぎ理論に関する考察
学籍番号 22413515
氏名 岩月 聡
指導教員名
1 はじめに
杉山 勝
確な解析が可能な ET におけるゆらぎ理論は、
先にナ
近年、ミクロ・ナノスケールでの加工や制御技術
ノジェットなどで例を挙げた通り工学的にも必要性
といった工学技術への関心が高まっている。
例えば、
が高まっている。そして、ET における流体力学的ゆ
ナノスケールの口径を持つ液体ジェット(ナノジェ
らぎの理論として、
「拡張された熱力学に基づく流体
ット)は、
高密度小型回路のプリンティングなどに用
力学的ゆらぎ理論」[3]が最近発表された。
いることができると考えられている。このような微
本研究の目的は、ゆらぎに対する深い理解を得る
小なスケールでは、分子や原子のランダムな運動に
ことである。そのために、２つのゆらぎ理論のそれ
起因する「ゆらぎ」の影響が無視できなくなる。そ
ぞれのレビューを通して、ゆらぎの源泉であるラン
こで、ゆらぎに対する理解を深めることが大変重要
ダムフォースをどのように導入すればいいかを確認
になってきた。
した。また、２つの理論の間でゆらぎの影響や大き
流体におけるゆらぎを解析するための理論として
は1950 年代にL.D.Landau とE.M.Lifshitz により提
さの違いを評価し、ゆらぎの本質に迫るための方法
を考察した。
唱された「Landau-Lifshitz(LL)理論」[1]がある。
この理論はナノジェットの数値シミュレーション
2 研究内容
[2]等の流体の研究で注目を浴びている。しかし、LL
本論文ではまず、
ゆらぎ理論の準備として TIP と単
理論には適用限界がある。LL 理論は「不可逆過程の
原子希薄気体における ET について説明した。また、
熱力学(TIP)」
に基づいた流体力学的ゆらぎ理論であ
ET の前に「気体分子運動論」を述べている。希薄気
る。そのため、TIP において採用されている、
「系全
体においては ET と気体分子運動論の結果は一致す
体は非平衡にあるが局所的には平衡状態にある」と
るため、ET の理論構造を理解するためにも解説した。
いう「局所平衡の仮定」が破れるような領域では用
また、ゆらぎを解析するにあたり本論文では単原子
いることができない。特に、ナノスケールのような
希薄気体を扱っている。
とても小さい系においてはこの仮定は必ずしも成り
立つとは限らない。
次に、ゆらぎについて、平衡状態における安定性
条件と熱的ゆらぎを記述する Einsetein のゆらぎの
これらの問題点を克服し、適用限界を超えた領域
確率分布の説明から始めた。次に、ブラウン運動に
でも成り立つ理論として「拡張された熱力学(ET)」
おける Langevin モデルを議論し、
流体力学的ゆらぎ
[3]がある。
ETではTIPで扱う独立変数、
質量密度ρ、
の基礎的な考え方について述べた。
速度𝑣𝑖 、温度𝑇に加えて散逸的な流束である応力テ
ンソル𝑡𝑖𝑗 、熱流𝑞𝑖 も独立変数として扱い、これらの
こ の Langevin モ デ ル を 参 考 に し て
「Landau-Lifshitz(LL)理論」と「拡張された熱力学
量のバランス方程式を採用している。局所平衡の仮
に基づく流体力学的ゆらぎ理論」を説明する。これ
定が成り立たない小さいスケールの現象にも適用で
らの流体力学的ゆらぎ理論を導くためには、それぞ
きる ET の立場でゆらぎを解析すれば、LL 理論で見
れの基礎式にどのようにランダムフォースを導入す
ていたゆらぎについて新しい知見を得ることができ、
るかを考えることである。
ET の解像度のレベルに合ったより小さなゆらぎを
LL 理論では、エントロピー生成に着目してランダ
追うことができるようになると考えられる。このよ
ムフォースを導入した。エントロピー生成Σは熱力
うに、ET におけるゆらぎ理論は、ゆらぎの本質にせ
学的な力𝑋𝑎 と流束𝑥̇ 𝑎 で
まる意味で必要である。また、ナノスケールでも精
Σ = −𝑋𝑎 𝑥̇ 𝑎
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
と表される。ここから、𝑋𝑎 と𝑥̇ 𝑎 から構成式が
𝑠<𝑖𝑗> =
𝑥̇ 𝑎 = −𝛾𝑎𝑏 𝑋𝑏
と仮定できる。𝛾𝑎𝑏 は現象論的係数である。この構
𝑠𝑝𝑝𝑖 = −
成式は線形の緩和方程式になっているので、
「ゆらぎ
の平均的な振る舞いは現象論的な発展方程式に従う」
という「線形減衰の仮説」を採用するとランダムフ
ォース𝑦𝑎 は Langevin 型に
𝑥̇ 𝑎 = −𝛾𝑎𝑏 𝑋𝑏 + 𝑦𝑎
と導入できる。𝑦𝑎 は Gauss 型の白色雑音である。
ET のゆらぎ理論でもLL 理論と同様にエントロピー
生成に着目する。ET におけるエントロピー生成
Σ = λ<𝑖𝑗> 𝑠<𝑖𝑗> + λ𝑝𝑝𝑖 𝑠𝑞𝑞𝑖
1
𝑡
+ 𝑏<𝑖𝑗>
𝜏𝑡0 <𝑖𝑗>
ここで、𝑎 ≡
𝑘𝐵
𝑚
2
𝑞 + 𝑐𝑖
𝜏𝑞0 𝑖
であり、𝑘𝐵 はボルツマン定数、𝑚は
分子の質量、𝜏𝑡0 , 𝜏𝑞0 は緩和時間、下付き添え字の 0
は平衡状態における各々の物理量の値を表す。
𝑓<𝑖𝑗> , 𝑔𝑖 , 𝑏<𝑖𝑗> , 𝑐𝑖 はそれぞれランダムフォースであ
る。𝑡<𝑖𝑗>は𝑡𝑖𝑗 のトレースレス成分を表している。そ
して、
粗視化の近似を用いることで ET のゆらぎ理論
から LL 理論に帰着することも確認されている。[4]
から線形の構成式が仮定される。
𝑠<𝑖𝑗> , 𝑠𝑞𝑞𝑖 は生成項、
最後に、
２つの理論におけるゆらぎの大きさの違い
𝜆<𝑖𝑗>, 𝜆𝑝𝑝𝑖 はラグランジュ未定乗数である。よって、
を定量的に調べようとした。そこで、数値解析を提
この構成式に Langevin 型のランダムフォースを導
案した。２つの理論の構造の大きな違いとして、LL
入すればいいことがわかる。
理論ではランダムフォースが空間微分で入っている
この考えに基づいて得られたゆらぎを含んだ基礎
のに対し、ET では生成項に和の形で入っている。こ
式は以下の通りである。今は、緩和方程式が線形の
れは数値解析を行う上で ET の理論は扱いやすいと
範囲内で考えているので基礎式は線形化してある。
考えられる。LL 理論の数値計算に関しては先行研究
・LL 理論
として J.B.Bell らの論文[5]がある。この論文の手
∂ρ
∂𝑣𝑖
+ 𝜌0
=0
∂t
∂𝑥𝑖
ρ0
∂𝑣𝑖
∂ρ
∂𝑇 ∂𝑡<𝑖𝑗>
+ 𝑎𝑇0
+ 𝑎𝜌0
−
=0
∂t
∂𝑥𝑖
∂𝑥𝑗
∂𝑥𝑗
3
∂𝑇 ∂𝑞𝑘
∂𝑣𝑘
𝑎ρ
+
+ 𝑎ρ0 𝑇0
=0
2 0 ∂t ∂𝑥𝑘
∂𝑥𝑘
𝑡<𝑖𝑗>
∂𝑣<𝑖
= 2μ
+ 𝑓<𝑖𝑗>
∂𝑥𝑗>
𝑞𝑖 = −κ
∂𝑇
+ 𝑔𝑖
∂𝑥𝑖
順を参考にして ET においても同様の数値計算を行
い、
論文の LL 理論の結果と比較するという方針を提
案した。
3 まとめ
本研究では、LL 理論と ET に基づく流体力学的ゆ
らぎ理論のレビューを通じて、TIP と ET の２つの理
論へのランダムフォースの導入方法をレビューした。
また、二つの理論の間でゆらぎにどのような関係が
あるかを述べた。そして、この双方の理論における
ゆらぎの大きさの違いを知るために、数値計算の手
・ET に基づく流体力的ゆらぎ理論
∂ρ
∂𝑣𝑖
+ 𝜌0
=0
∂t
∂𝑥𝑖
ρ0
∂𝑣𝑖
∂ρ
∂𝑇 ∂𝑡<𝑖𝑗>
+ 𝑎𝑇0
+ 𝑎𝜌0
−
=0
∂t
∂𝑥𝑖
∂𝑥𝑗
∂𝑥𝑗
3
∂𝑇 ∂𝑞𝑘
∂𝑣𝑘
𝑎ρ0
+
+ 𝑎ρ0 𝑇0
=0
2
∂t ∂𝑥𝑘
∂𝑥𝑘
∂𝑡<𝑖𝑗> 4 ∂𝑞<𝑖
∂𝑣<𝑖
−
− 2𝑎𝜌0 𝑇0
= −𝑠<𝑖𝑗>
∂t
5 ∂𝑥𝑗>
∂𝑥𝑗>
∂𝑡<𝑖𝑗> 5 2
∂𝑞𝑖
∂𝑇
− 𝑎𝑇0
− 𝑎 𝜌0 𝑇0
= 𝑠𝑝𝑝𝑖
∂t
∂𝑥𝑗
2
∂𝑥𝑖
法を提案した。
4 参考文献
[1] ランダウ，リフシッツ：『流体力学２』
（東京図書、1971）
[2] M.Moseler and U.Landman : Science 289 (2000),1165.
[3] I. Muller and T. Ruggeri: Rational Extended
Thermodynamics, (Springer-Verlag, New York, 1998).
[4] A.Ikoma, T.Arima, S.Taniguchi, N.Zhao,and Masaru
Sugiyama : Physics Letters A 375 (2011), 2601.
[5] J.B. Bell, A.L. Garcia, and S.A. Williams: Phys. Rev.
E 76 (2007) 016708.
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
非飽和状態における IEEE802.16 無線通信ネットワークの性能解析
学籍番号
1 はじめに
IEEE802.16は，見通し距離50㎞で最大速度70Mbps
22413516
氏名
上田
指導教員名
馮
晋大
偉
る．ここで定数W0 はバックオフウィンドウの初期
の通信速度を実現する広域・高速の通信規格である．
値である．バックオフカウントは 1 コンテンション
現在，
その派生規格であるIEEE802.16-2004を用いた
スロット経過する毎に 1 ずつカウントダウンされ，
0
Wi-MAXという通信サービスが提供されている．
に到達すると帯域幅割当要求を送信する．
IEEE802.16の性能評価に関する研究は規格の策定
帯域幅割当要求を送信した端末は基地局からの帯
と前後して広く行われるようになっているが，多く
域幅割当を待つ状態となる．利用できる帯域幅があ
の研究は飽和状態（バッファに常に送信データがあ
れば，基地局は端末に帯域幅を割り当てる． M フ
る状態）の仮定の下で解析が行われており，バッフ
レーム以内に帯域幅割当が行われなければ，帯域幅
ァに保存されるパケット数の変化を無視していた
割当要求の再送のために再びバックオフを行う．
（例えば参考文献[1]，[2]）．しかし，バッファの
安定性，遅延やオーバーフローの確率は通信ネット
ワーク品質の重要な尺度である．
3 マルコフ連鎖によるモデル化
本研究では 1 基の基地局と n 台の端末から構成さ
そこで，本研究ではバッファ容量の性能評価を行
れる IEEE802.16 無線通信ネットワークについて考
うため，非飽和状態の仮定の下でIEEE802.16無線通
える．Fallah[1]のモデルを基に，ある特定端末（マ
信ネットワークの数理モデルを提案し，その性能解
ーク端末）の状態を 3 次元の離散時間マルコフ連鎖
{ X (t ), t = 0,1,2,L} によってモデル化する．時間推
析を行った．
移はコンテンションスロットを単位として起こるも
2 IEEE802.16 の通信方式
のとする．離散時刻 t におけるマーク端末のバッフ
IEEE802.16 標準規格では一対多接続方式が想定
ァ内のパケット数，再送回数，バックオフカウント
され，1 つの基地局がカバーエリア内のすべての端
の値もしくは帯域幅割当要求送信後の経過スロット
数を表す確率過程を x(t ), s (t ), b(t ) とする．これよ
末の帯域割当等の制御を行う．ネットワーク内のす
べての端末は，基地局から割り当てられた帯域幅を
利用してデータ送信を行うため，まず基地局に帯域
幅の割り当てを要求しなければならない。
帯域幅割当要求の衝突を回避するため，各端末は
り離散時刻 t におけるマーク端末の状態は
X (t ) = ( x (t ), s (t ), b(t )) と書き表すことができ，
X (t ) は 3 次元の離散時間マルコフ連鎖となる．
このマルコフ連鎖の推移確率行列 P を求めると
帯域幅割当要求の送信前にバックオフを行う．バッ
次のようになる．
クオフはコンテンション期間で行われる．コンテン
⎡C0 A1 A 2 A 3 A 4 L⎤
⎢C B B B B L⎥
0
1
2
3
⎢ 1
⎥
⎢
C2 B 0 B1 B 2 L⎥
P=⎢
⎥
C2 B 0 B1 L⎥
⎢
⎢
C2 B 0 L⎥
⎢
⎥
O O⎥⎦
⎢⎣
これは M / G / 1 型マルコフ連鎖の推移確率行列と
ション期間はコンテンションスロットと呼ばれる固
定時間単位の列として成り立っている．帯域幅割当
要求を送信したい端末は，まず，バックオフカウン
トと呼ばれるランダムな整数値を選ぶ． i 回目の再
送では 0 からWi − 1 の範囲から選ばれる．Wi はバッ
クオフウィンドウと呼ばれる値で，Wi = 2i W0 であ
なっている．
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
4
モデルの解析
バッファ 内の平均パケット数E[L]
端末の状態の定常分布を
5
π ( k ,i , j ) = lim P(( x(t ), s(t ), b(t ) = (k , i, j ))
4
π 0 = lim P( x(t ) = 0)
t →∞
パケット数
t →∞
3
2
1
と定義し，M / G / 1 型推移確率行列の定常分布を導
0
2
4
6
8
10
12
14
16
18
20
22
24
26
28
30
32
34
36
38
40
42
44
46
48
50
出する標準的なアルゴリズムを用いて定常分布を求
N
N
N
N
N
N
=
=
=
=
=
=
5
10
15
20
25
30
n
める．求めた定常分布からパケット数に注目した定
ネットワーク内の端末数
常確率ベクトル
図 1：N, n に対する平均パケット数 E[L]
πk = (π ( k ,0,0 ) , π ( k ,0,1), Lπ ( k , m, ( M − 2) )
バッファ内の平均パケット数E[L]
を求め，定常状態におけるバッファ内の平均パケッ
ト数 E[ L] を求める．
5
∞
E[ L] = ∑ kπk
パケット数
4
k =0
3
2
1
0
数値計算結果と考察
2
4
6
8
10
12
14
16
18
20
22
24
26
28
30
32
34
36
38
40
42
44
46
48
50
5
M=2
M=5
M=7
M=9
M = 11
n
3 次元のマルコフ連鎖でモデル化したところ，状
ネットワーク内の端末数
態数が膨大なものとなった．そこで、今回の数値計
図 2：M, n に対する平均パケット数 E[L]
算では計算量の抑制のため、バッファ容量を 5，
W0 = 8 とした．図 1 は 1 フレームあたりのコンテン
IEEE802.16無線通信ネットワークのモデルを提案し
ションスロット数 N とネットワーク内の端末数 n
に対する平均パケット数 E[L] を示したものである．
性能解析を行った．
モデルが3次元のマルコフ連鎖と
N ≥ 15 ではほぼ同じ値となったが N の値が小さ
制するために推移確率行列 P のサイズを制限した
いときには平均パケット数は増加する．特に
計算しか行えず一部の定性的な評価を行うのみとな
N < W0 では著しく増加する．図 2 は帯域幅割当を
った．今後の課題としては、近似による状態数の削
待つ最大フレーム数 M とネットワーク内の端末数
n に対する平均パケット数 E[L] を示したものであ
列 P のサイズ制限を広げ定量的な評価を行う必要
なったことで状態数が非常に多くなり，計算量を抑
減やより速い計算方法を用いるなどして推移確率行
がある．
る．図 2 を見ると M の値が小さい方が平均パケッ
ト数は少なくなるように考えられる．今回の計算で
7 参考文献
は，帯域幅割当要求の衝突確率が非常に高くなり，
[1] Y.P.Fallah, F.Agharebparast, M.R. Minhas, H.M.
複数回の衝突によって M の値が大きい方が処理時
Alnuweiri :“Analysis Modeling of Contention-Based
間のロスが大きくなったためこのような結果になっ
Bandwidth Request Mechanism in IEEE802.16 Wireless
たと考えられる．しかし，全体的に衝突率の高さと
Networks”, IEEE Trans. Veh. Technol., vol.57, no.5,
送信成功率の低さが著しかったため，今回の結果の
pp.2094-3107, (2008)
みではシステムパフォーマンスのパラメータへの依
[2]
存度は結論付けることができない．
Comprehensive Analysis of Bandwidth Request
D.Chuck,
K-Y.Chen,
J.M.Chang
:“A
Mechanisms in IEEE802.16 Networks”,IEEE Trans.
6 まとめ
本研究ではバッファの状態を考慮に入れた
Veh. Technol., vol.59, no.4 (2010)
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
最近傍分類のための適応的な距離学習アルゴリズム
学籍番号 22413517 氏名 植田覚
指導教員名 竹内一郎 准教授
1
はじめに
近年, 計測技術やネットワーク技術などの発達によ
り, 大規模なデータを取得出来るようになった. これら
のデータの背後に潜む規則を抽出する技術は機械学習
と呼ばれ, 様々な分野でその重要性が認識されている.
多くの機械学習アルゴリズムは距離情報を用いるため,
問題に応じて適切な距離尺度を導入することが重要と
なる. 本研究では, 学習データを用いて距離尺度の推定
を行う「距離学習」について考察する. 特に, クラス分
類問題の基本アルゴリズムである最近傍分類法に関す
る距離学習の問題を考察する.
最近傍分類法とは, クエリとして与えられた入力に対
して, その近傍点を見つけ, そのクラスラベルの多数決
によってクラス分類を行う手法である. 最近傍分類法
は距離尺度に依存するため, 最近傍分類法のための様々
な距離学習アルゴリズムが提案されてきた. これらの
従来法の多くは, 同じクラスの近傍点を近づけ, 異なる
クラスの近傍点を遠ざけるように設計されている. し
かし, 近傍点が学習前の距離尺度によって決定するた
め, 学習時に距離尺度が変化すると, 近傍点が離れてし
まったり, その逆が起こる問題点がある. そこで本研究
では, 学習途中で近傍点を適応的に更新する距離学習ア
ルゴリズムを導入し, その性能を実験的に検証する.
2
問題設定
サンプルサイズ n, 特徴数 p の G クラス分類問題を考
える. 学習データを D := {(xi , yi )}i∈Nn , xi ∈ Rp , yi ∈
{1, · · · , G} とする. ここで, xi は p 次元特徴ベクトル,
yi はクラスラベルを表し, Nn := {1, 2, · · · , n} である.
距離学習 [1, 2] とは, 最近傍分類法などの分類性能の
向上を意図した手法で, その多くはマハラノビス行列を
学習するものである. マハラノビス距離は, 任意の 2 つ
のインスタンス xi と xj の距離を
⊤
d(xi , xj |M ) := (xi − xj ) M (xi − xj ),
(1)
と表す. ここで, M ∈ Rp×p は距離尺度を定義するマ
ハラノビス行列で, 距離尺度の概念より, 半正定値性及
び対称性を持たなければならない.
多くの距離学習アルゴリズムはマハラノビス行列に
関する最適化問題として定式化される. 最近傍分類法
の性能向上のために, 分類誤差を小さくする M を求め
ることが目的となる.
2.1
target neighbors(TNs)
最 適 化 問 題 を 定 式 化 す る た め に, target neighbors(TNs) を導入する. ここでは, 2 クラス分類問題
を考える. 2 クラスの最近傍分類法では, インスタンス
xi の k 個の近傍点のクラスの多数決で分類を行う. そ
のため, 最近傍分類法の性能を左右するのは, 各クラス
における (k + 1)/2 番目の近傍点であることがわかる.
各インスタンス xi において, 同じクラスの (k + 1)/2
番目の近傍点を target hit と呼び, hi と表す. また, 異
なるクラスの (k + 1)/2 番目の近傍点を target miss と
呼び, mi と表す. 以下では, target hit と target miss を
合わせて, target neighbors (TNs) と呼ぶことにする.
また, 多クラス分類問題においても同様な定式化が可
能であるが, スペースの都合上省略する.
距離学習問題の定式化
2.2
最近傍分類法が高い汎化性能を持つために, target hit
までの距離よりも target miss までの距離を大きくした
い. そのため, 損失関数を全てのインスタンスにおいて,
target hit までの距離を小さくし, target miss までの距
離を大きくする距離尺度を得るように定式化する. さ
らに, 正則化項を導入し, M の制約条件を考慮すると,
本研究で考察する距離学習問題は,
min
M
θn−1
n
∑
{d(xi , xhi |M ) − d(xi , xmi |M )}
i=1
s.t.
+2−1 ||M − I||2F
(2a)
M ≽ 0,
(2b)
と定式化される. ここで, 目的関数 (2a) の第 1 項は損失
関数を表す. 第 2 項は正則化項を表し, ユークリッド距
離からの乖離の程度をフロベニウスノルムで定量化し
た. また, θ ∈ [0, ∞) は正則化パラメータで, 損失関数
と正則化項のトレードオフを制御している. θ = 0 のと
き最適化問題 (2) の解はユークリッド距離尺度となる.
3
提案法:適応的な距離学習
TNs を固定すると, 最適化問題 (2) が半正定値計画
問題となるため問題の扱いが容易となる. しかしなが
ら, 距離学習では学習の途上で行列 M が変化するため
距離尺度そのものが変更されてしまう. TNs は M に
より定義されるマハラノビス距離 (1) に基づいて定義
されるため, 距離尺度が変更することで学習前の距離尺
度での TNs が学習後の距離尺度でも TNs であるとは
限らなくなってしまう. したがって, 距離尺度に応じた
TNs の適応的な変更が必要である.
3.1
適応的な距離学習
提案法では距離尺度に応じて適応的に TNs を変更す
るアルゴリズムを考える. 提案法のアルゴリズムの概
要を説明する. まず, ユークリッド距離尺度 (M = I)
で初期化し, 初期の TNs を割り当てる. 学習が進むに
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
つれて距離尺度 M が変化すると, TNs の違反 (本来な
らば TNs でないものが TNs となる状況) が生じる. 違
反の程度がある基準を越えたとき, その時点の M にお
いて TNs の再割り当てを行い, さらに学習を続ける.
TNs の再割り当てを行う M を, 順に, M1 , M2 , . . .
とすると, t 回の再割り当てが行われた後では, 最適化
問題 (2) の単位行列 I を Mt に変更した問題:
min
M
θn−1
i=1
−1
+2
s.t.
n
∑
{d(xi , xhi |M ) − d(xi , xmi |M )}
||M − Mt ||2F
(3a)
M ≽ 0,
実験結果
4.2
図 1 は k = 5 の場合に選択された M をテストデー
タに用いた時の最近傍分類法の分類誤差の平均を示し
ている. 値が小さいほど良い結果であると言える.
図 1 より, 全てのデータにおいて提案法が最良の結
果となった. そのため, 提案法が有効に働くことがある
ことを示唆している. また, 実験において, TNs の違反
の基準を緩くしたほうがより良い結果を得ることがで
きた. これは, TNs の違反の基準を厳しくすると TNs
の更新が頻繁に起こるため, 学習が進まないのではない
かと考えられる.
(3b)
30
Eclid
NCA
を解く. θ = 0 のとき最適化問題 (3) の解は Mt となる.
t = 1, 2, . . . において最適化問題 (3) を解く際には,
TNs の違反が基準を越えないような最大の θ > 0 を見
つける必要がある. 提案法では, 簡単な二分探索を用い
ている. また, 過学習を防ぐため, 学習プロセスで得ら
れる全ての M に対して評価サンプルの最近傍分類誤
差を計算しておき, それが最小となる M を選ぶ.
4
25
error rate[%]
15
10
5
0
BreastCancerDiagnostic
実験概要
データは UCI Machine Learning Repository から入
手した. 表 1 に実験で用いたデータの概要を示す. 各
データをランダムに学習データ, 評価データ, テスト
データに分割する. データ数が 300 以下の場合には,
データを学習データと評価データとテストデータに 3
分割し, それ以上の場合には, 学習データ数と評価デー
タ数を 100 とし, 残りをテストデータとした. 学習デー
タで距離尺度を学習し, 評価データを用いて過学習を防
ぐ汎用性の高い解を選択する. そして, 選択された距離
尺度をテストデータに用いて最近傍分類法の分類誤差
を計算する. ランダムな影響を軽減するため, データの
分割を 10 通り行い実験した. 提案法の終了条件は, 損
失関数の変化が十分に小さくなったときとした.
表 1: データの概要
データ名
n
BreastCancerDiagnostic 569
Ionosphere
351
Parkinsons
195
SPECTFHeart
267
adaptive
20
計算機実験
提案法のベンチマークデータに対する数値実験結
果を報告する. 提案法 (adaptive) と他の距離学習
(NCA[1],LMNN[2]) とユークリッド距離 (Euclid) で得
られた距離尺度を用いて, 最近傍分類法の分類誤差を比
較する.
4.1
LMNN
p
30
33
22
44
G
2
2
2
2
Ionosphere
Parkinsons
SPECTFHeart
図 1: k = 5 における最近傍分類誤差の平均
5
まとめ
本研究では提案法として最近傍分類法のための適応
的な距離学習アルゴリズムを提案した. ベンチマーク
データに対する数値実験により提案法が状況によって
有効に働くことを検証した. 今後の課題として, 今回は
TNs の違反の基準を指定したが, 学習を行なう過程で
適応的に決定するアプローチを考えることなどが挙げ
られる.
参考文献
[1] J. Goldberger, S. Roweis, G. Hinton, and R.
Salakhutdinov. Neighbourhood components analysis. In L. K. Saul, Y. Weiss, and L. Bottou,
editors, Advances in Neural Information Processing Systems 17, pages 513-520. MIT Press, Cambridge, MA, 2005.
[2] K. Q. Weinberger and L. K. Saul. Distance Metric
Learning for Large Margin Nearest Neighbor Classification. Journal of Machine Learning Research
10. 2009.
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
集団送信 DCF を用いた IEEE802.11 無線通信ネットワークの性能解析
学籍番号
22413519
氏名
指導教員名
1 はじめに
内田
亮
馮
偉
非飽和状態を考えるために，バッファの情報を取り
IEEE802.11規格は，
主に無線LANで用いられる規格
入れた．これにより，ステーションの状態はバッフ
である．無線通信システムの複雑化に伴い，理論解
ァにあるパケットの数，バックオフステージの数，
析や数値計算による性能評価が重要視されるように
なってきた．これらに関する研究論文は，[1]-[3]
バックオフカウンタの値の 3 つとなる．時刻 t にお
ける，それぞれの確率過程を x(t ) ， s (t ) ， b(t ) と
がある．IEEE802.11規格の主なアクセス方法として
す る と ， あ る ス テ ー シ ョ ン の 状 態 は X (t ) =
DCF(分散コーディネーション機能)がある．Bianchi
( x(t ), s (t ), b(t )) と表すことができ，3 次元の離散時
氏の論文[1]は，DCFに関する研究の先駆けとなって
間マルコフ連鎖となる．このマルコフ連鎖の推移確
いる．また，Meerja氏らの論文[2]では連続送信の考
率行列は以下のようになる．(例:U = N = 3 )
えを中心に基本DCFに修正を加え，
通信速度の向上を
実現している．これらの研究は，バッファにパケッ
トが常にある(飽和状態)という仮定の下で解析を行
っている．しかし，実際の無線通信システムで常に
送信したいパケットが存在することは考えにくい．
これに対し，Bae氏らの論文[3]では，基本DCFに関し
て，パケットに到着率を与えることで非飽和状態を
作り出し，解析を行っている．
⎛ C0
⎜
⎜ C1
⎜ 0
P=⎜
⎜ D0
⎜ 0
⎜
⎜ M
⎝
A1
A2
A3
A4
B1
B2
B3
B4
C2
B1
B2
B3
K1
E0
E1
E2
K0
K1
E0
E1
M
M
M
M
L⎞
⎟
L⎟
L⎟
⎟
L⎟
L⎟⎟
O⎟⎠
3 モデルの解析
以上を踏まえ，本研究では，通信速度の向上のた
はじめに，推移確率行列をU ブロック毎に区切り
めに集団送信の考えを取り入れ，より現実の無線通
再構築する．再構築した行列は，M/G/1 型の行列と
信システムに近い非飽和状態の下で解析を行った．
なるので，標準的なアルゴリズムに沿って，定常分
布を求める．定常分布から処理量や HoL-delay など
2 モデルの設定
のシステムの性能解析を行う．
2.1 M-DCF
3.1 処理量
基本 DCF では，チャネルの状態に関わらず，必ず
バックオフ過程を経て送信を行っていたため，バッ
非飽和状態における処理量は，1 つのスロットの
間に送信できるペイロードの割合によって表される．
クオフ時間の無駄が多く存在した．無駄を減らすた
めに，本モデルでは，集団送信の考えを中心に基本
Γ=
DCF を修正した M-DCF を提案する．基本 DCF では，1
個のパケットを送信する度にバックオフを行ってい
( Ps1 + NPs N ) E[ P]
Pi Ti + Pc Tc + Ps1Ts1 + P TsN
(2)HoL-delay
たのに対し，M-DCF では，最初の送信が成功し，且
HoL-delay は，あるパケットの順番が回ってきた
つバッファにパケットがU 個以上ある場合， N 個
地点から，送信を終えるまでの時間を表す．本モデ
のパケットをまとめて送信することとした．これに
ルでは，各パケットの HoL-delay は，パケットが到
よってバックオフ時間の削減を期待する．
着した地点におけるチャネルまたはバッファの状態
2.2 3 次元マルコフ連鎖によるモデル化
によって異なる．最終的な式は次の式となる．
あるステーションの動きに注目し，M-DCF を元に
無線通信システムのモデル化を行う．
本モデルでは，
E[ S H ] = PN {Pη1 E[ S 0 ] + A} + PS A
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
は低下する．また，ステーションの数が35以上にな
ここで、 A は次式である。
(
~
A = Pη1 (1 − ξ )( P01Ts1 + P0N Ts1 ) + ξ (Tc + E[ S1 ])
)
ると，N が大きいほど処理量が高くなっている．こ
の理由として，集団送信によってバックオフ時間が
減ったため，全体に占める送信の割合が増えたこと
⎛ N −1 ~
⎞
+ (1 − ξ ) P ⎜ ∑ Pηl Ts1 + PηN Ts1 ⎟
⎝ l =2
⎠
N
0
が考えられる．
4.2 HoL-delay
ステーションの数が35 を超えたところから
4 数値計算結果と考察
バックオフステージの最大
初期CWサイズω 0 = 8 ，
HoL-delayが急激に大きくなっている．理由として，
このステーション数から衝突が増え，バックオフ時
数 m = 3 ，パケットの到着率 λ = 15 とし，ステー
間が長くなったことが考えられる．また，N が大き
ションの数を5から60まで5刻みに変化させたときの，
いほどHoL-delayが短くなっている．これは，処理量
集団送信するパケットの数 N と処理量及び
でも述べたように，集団送信によって，バックオフ
HoL-delayの関係をそれぞれ図1，2に示す．
時間が減ったためと考えられる．集団で送信される
とき，2番目以降の N − 1 個のパケットは，バックオ
処理量
フ時間が無い．これによって，1個のパケットあたり
0.6
の送信時間は短くなったと考えられる．
0.5
0.4
N=1
0.3
N=3
5 まとめ
N=5
0.2
本研究では，バッファの情報取り入れたことによ
0.1
り， N やU を調整することが可能となった．数値
0
5
10
15
20
25
30
35
40
45
50
55
60
計算結果では，N が大きいほど処理量が高くなるこ
ステ ー シ ョ ン の 数
とや，
HoL-delayが短くなるといった結果が得られた．
これは，集団送信を行うことによってシステムのパ
図1 : N と処理量
フォーマンスが向上したことを意味する．
HoL-delay
今後の課題として，行列解析法のアルゴリズムの
0.35
見直しが挙げられる．行列をU ブロック毎に区切り
0.3
再構築したため，U が大きいほど数値計算に時間を
0.25
N=1
0.2
N=3
0.15
N=5
要する．到着率やステーションの数に応じたU や
N の評価を行うために，アルゴリズムを改善して，
U を大きく設定できるようにする必要がある．
0.1
0.05
0
5
10
15
20
25
30
35
40
45
50
55
60
ステーションの数
6 参考文献
図2 : N とHoL-delay
4.1 処理量
処理量はステーションの数が30のところでピーク
[1] G.Bianchi :“Performance analysis of the IEEE802.11
Distributed Coodination Function”,IEEE J.Sel.Areas
Commun., vol.18, no.3, pp.535-547, 2000
を迎え，その後下がっている．これは，先に述べた
[2] K.M.Meerja and A.Shami:“Analysis of new
パラメータにおいては，ステーションの数が30のと
distributed-media access-control schemes for IEEE802.11
きに最もシステムのパフォーマンスが高いことを意
wireless local-area networks”,IEEE Trans.Veh.Technol.,
味する．ピーク後に下がる理由として，パケットの
vol56, no.4, pp.1797-1812, 2007
衝突率の増加が考えられる．競合するステーション
[3]Y.H.Bae, K.J.Kim, M.Moon, and B.D.Choi :“Analysis of
の数が増えるため，
送信率が上がり，
衝突が増える．
IEEE802.11 non-saturated DCF by matrix analytic
これにより，送信が成功する確率が下がり，処理量
methods”,Ann.Oper.Res., vol.162, pp.3-18, 2008
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
SVM パス追跡の GPGPU による並列処理
学籍番号 22413521 氏名 大坪拓也
指導教員名 竹内一郎 准教授
1
はじめに
の最適性条件は KKT 条件から次のように表せる:
本研究では 2 クラス分類問題を解く代表的な手法で
あるサポートベクトルマシン (SVM) に並列処理を適
用することを考える．SVM は未知の問題に対する分類
性能に優れることから，多くの 2 クラス分類問題に汎
用的に用いられてきた．しかし SVM の学習にはデー
タ数を n としたとき一般的に O(n2 )∼O(n3 ) の計算時
間を要することが知られている．そこで近年並列処理
の手法として注目されている GPGPU に着目する．
これまでに一般的な SVM 解法である SMO アルゴリズ
ムを GPGPU により高速化する研究などが発表されて
いる [1]．SMO アルゴリズムは解の一部を逐次的に更
新するため，GPGPU による並列化の利点を十分に享
受できない．一方，別の SVM 解法として正則化パス追
跡と呼ばれるアルゴリズムがある. このアルゴリズム
では SMO より計算コストが大きいが，カーネル関数行
列を用いた行列演算を大量に含むため，GPGPU によ
る効率化が期待できる．本研究では，文献 [2]，[3] で提
案されている正則化パス追跡アルゴリズムに GPGPU
による並列処理を適用し、その有用性の検証を行う．
2
SVM
SVM とは教師あり学習の 1 つであり，代表的な 2 ク
ラス分類問題の 1 つである．n 個のデータ点として特
徴ベクトル xi ∈ Rd ，クラスラベル yi ∈ {1, −1}, i =
1, · · · , n, が与えられているとする．ラグランジュ未定
乗数より得られる双対問題は次のようになる:
n
n
1 ∑∑
max −
α
2
s.t.
yi yj αi αj Kij +
i=1 j=1
0 ≤ αi ≤ C,
n
∑
n
∑
i=1
αi
(1)
yi αi = 0.
i=1
ここで α はパラメータベクトル，K はカーネル関数で
あり，Kij = K(xi , xj ) である．正則化係数 C は SVM
の過学習を回避するために調整するパラメータである．
双対問題において判別関数は，
f (x) =
n
∑
yi αi K(x, xi ) + w0 ,
(2)
i=1
となる．w0 はバイアスである．ここでマージン内のデー
タ点の集合を I ，マージン境界上のデータ点の集合を
M，マージン外のデータ点の集合を O とする．SVM
I:
M:
O:
2.1
yi f (xi ) ≤ 1 and αi = C,
(3)
yi f (xi ) = 1 and 0 < αi < C,
(4)
yi f (xi ) ≥ 1 and αi = 0.
(5)
SMO アルゴリズム
SVM 解法の従来法である SMO は SVM が凸計画問
題であることを活かしたアルゴリズムである．最適性
条件を最も満たしていないと思われるデータ点を 2 つ
選び，それらに対応するパラメータを暫定的な最適値
に更新する操作を繰り返すことで最適解に収束する．
正則化パス追跡アルゴリズム
3
本研究で扱う正則化パス追跡アルゴリズムは (1), (3),
(4), (5) の等式で表せられる最適化問題を利用する．最
適解 α と正則化係数 C の変化における連立方程式を
導くことで，C に対応した α の連続した変化を得る．
なおこれ以降重みつき SVM と同様に C を各データ点
に対応させたベクトル C として扱う．また正則化係数
の目標値までの差を ∆C とする．∀i ∈ M のとき (3)，
(4)，(5) の等式条件及び (2) より，
[
∆
]
w0
αM
[
=−
0
yM
y⊤
M
QMM
]−1 [
∗
y⊤
I
QMI
]
∆C I , (6)
として α の変化量 ∆α と C の変化量 ∆C の関係を連
立方程式で表すことができる．ただし Q = yy ⊤ ⊗ K ，
⊗ は要素ごとの積，添え字 O, M, I はそれぞれの集合
に対応した部分集合，部分ベクトルを表す．(6) から
目的の解を求める際に，O, M, I の集合が変化する地
点に遭遇することがある．これを breakpoint と呼び，
遭遇した地点で連立方程式を更新し直す必要がある．
breakpoint の検出には (3)，(4)，(5) より最も早く集
合を外れるデータ点を検出する．正則化パス追跡アル
ゴリズムの簡単な流れを以下に示す．
4
GPGPU
GPGPU とは，PC において CG 処理を行う GPU と
いうハードを汎用計算に用いる技術である．CPU と比
較して，小規模な計算装置を数百搭載し大量かつ単純
な並列処理が得意という特徴がある．このため行列演算
などについて非常に有効である．本研究では NVIDIA
社が提供する開発環境 CUDA[4] で実装を行っている．
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
Algorithm 1 SVM 正則化パス追跡
1: while C が ∆C 更新するまで do
2:
∆α を (6) より計算
3:
∆yf を (2) に ∆α を代入し計算
4:
breakpoint を (3)，(4)，(5) より検出
5:
α, C, yf を breakpoint まで更新
6:
集合 O, M, I を更新
7:
(6) 右辺の係数行列及びベクトルを更新
8: end while
4.1
正則化パス追跡と GPGPU
正則化パス追跡では 1 ループのコストが高い．大きな
計算コストとなるのは，カーネル関数行列がメモリに
収まらない場合の再計算と ∆yf の計算である．これら
の計算は行列演算で表現でき，カーネル関数行列では n
行 n 列の行列生成を，∆yf では n 行 n 列行列× n 次元
ベクトルを計算する．正則化パス追跡は breakpoint 数
の予測が困難なため，最終的な計算コストの予測も困
難である．しかしアムダールの法則より，並列処理を行
える部分が大きい程最大限に並列処理を行った場合の
性能向上率も高くなる．1 ループに O(n2 ) の纏まった行
列演算を含む正則化パス追跡は並列処理を行える処理
の割合が高いため，大量の並列処理を行える GPGPU
と相性が良い．また正則化パス追跡は既存の解から目
的の解まで移動する量が少ない場合，ループ数が少なく
済む傾向がある．この場合全体の計算量を抑えながらも
大量のデータを用いた計算を少ない回数行う，GPGPU
に適した処理になる．
実験
5
今回は大規模データにデータを追加した際，提案法及
び SMO を CPU 及び GPU で計算した結果を紹介する．
5.1
実験概要
• 使用データ
– ”Web” データ数 49249 (初期解分) + 500 (追
加分)，次元数 300
– C = 20, γ = 1/300，tolerance=5 ∗ 10−5
SMO
パス追跡
CPU
ループ数
885.739 s
66086
52277.5 s
76
GPU
ループ数
215.187 s
72973
179.291 s
47
×4.54
×180.32
1 ループ当りの性能向上率
SMO アルゴリズムと GPGPU
SMO は 1 ループが軽量だが，1 ループで更新される
のは α の内 2 値と判別関数ベクトルであり，更新する
2 値の選択もベクトル操作で行う．加えてループ数が
大きくなりやすく GPU は小規模な計算を何度も行う
ことになる．そのため GPGPU は性能を伸ばし切れず，
性能を享受しきれないと考えられる．
4.2
表 1: Web：計算時間と性能向上率
• 実行環境
– CPU：Intel Core i7 860
– GPU：GeForce GTX480 (VRAM 1536MB)
5.2
実験結果と考察
結果は表 1 のようになった．カーネル関数行列のキャッ
シュに 2412 行，単精度浮動小数点で保存する．従来法
である SMO では高速化されてはいるものの大きな値
ではない. 一方正則化パス追跡ではカーネル関数行列
の計算が高速化されたことで 180 倍の性能向上が見ら
れた．
6
まとめ
本研究で，パス追跡を利用した SVM の適応的更新
法に GPGPU を適用することを提案した．また実験に
より，GPGPU が有効に働いていることを示した．
参考文献
[1] Austin Carpenter，CUSVM:A cuda implementation of support vector classification and regression,
2009,
(http://patternsonascreen.net/cuSVM.html).
[2] T. Hastie, S. Rosset, R. Tibshirani, and J. Zhu.
The entire regularization path for the support vector machine. Journal of Machine Learning Research, 5:1391-1415, 2004.
[3] M.Karasuyama, I. Takeuchi, Multiple incremental
decremental learning of support vector machine，
NIPS2009，2009．
[4] NVIDIA Corporation, CUDA
(http://www.nvidia.co.jp/object/cuda home new
jp.html)
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
テイラークエット流れにおける鎖状高分子の影響
学籍番号 22413531 氏名 金 相佑
指導教員名
後藤俊幸
である．ここで𝑾1 と𝑾2 は両端の球でのブラウン運
動によるランダム力である．ランラム力は次の条件
1 はじめに
を満足する．
テイラークエット流れは同軸 2 重円筒の間に流体
𝑊𝑖 𝑡
を入れて円筒を回転させたときの流れである．テイ
(𝑚 )
𝑊𝛼,𝑖
ラークエット流の数値シミュレーションを行うこと
で，乱流への遷移の様子や詳細な場の構造を知るこ
台の間の流体の種類による抵抗の関係を調べること
(𝑛)
𝑡 𝑊𝛽,𝑗 (𝑠) = 𝛿𝑚𝑛 𝛿𝛼𝛽 𝛿𝑖𝑗 𝛿(𝑡 − 𝑠)
p
またダンベルのストレステンソルTij は
とができる．テイラークエット流の安定性の研究は
動力を伝達する軸の回転速度とその軸を支える支持
=0
𝑝
𝑇𝑖𝑗
𝜈𝜂 𝐿3𝑏𝑜𝑥
=
𝜏𝑠 𝑁𝑡
𝑁𝑡
𝑛=1
𝑅𝑖 𝑅𝑗
𝑹
𝑓
2
𝑟𝑒𝑞
𝐿𝑚𝑎𝑥
− 𝛿𝑖𝑗 𝛿 𝑥 − 𝑟𝑔 (4)
によって工業上の応用も考えられる．さらに，この
で計算できる[1]．ここで 𝜂 ≡(3req/4a)2𝛷𝑉 ，𝛷𝑉 は全
流れに高分子を添加することで抵抗減少によるエネ
流体と全高分子の堆積比率である．
ルギー損失の減少が期待できる．
そこで本研究では，
テイラークエット流の Reynolds 数と撹乱によって
流れが不安定になることを確認し，高分子の濃度と
Weissenberg 数による流れの不安定性の変化と内側
の円筒に働く力との関係を解析する．
2 高分子モデルと基礎方程式
流体に関しては次式の非圧縮の Navier-Stokes 方
図1：ダンベルモデル．二つの球がバネでつながって
程式を円筒座標系に変換した式を用いる．そして高
いるモデルで，reqはダンベルの平衡長さ，Lmaxは最大
分子が流体に及ぶす力を外力として加えた．
長さ，rgは重心ベクトル，Rは末端間ベクトル，そし
𝜕𝒖
1
+ 𝒖 ∙ 𝛻 𝒖 = − 𝑝 + 𝛻2𝒖 + 𝒇
𝜕𝑡
𝜌
(1)
𝑝
𝑓𝑖 =
てa両端の球の半径を意味する．
3 並列計算
𝜕𝑇𝑖𝑗
本研究では多量とダンベルを計算するので計算時
𝜕𝑥𝑖
実際の高分子を計算機上で再現することは難しく，
間の短縮するために MPI を用いて並列計算を行う．
本研究では大量の高分子モデルを計算するので各高
全体の並列の概念図を図 2 に示す．本研究では流体
分子を図 1 で示す様なダンベルモデルを用いて表す． 場の計算とダンベル運動の計算を分離した．流体に
そして高分子の運動はこのダンベルの重心ベクトル
𝒓𝑔 と末端間ベクトル𝑹の方程式はそれぞれ
𝑑𝒓𝑔 1
1 𝑟𝑒𝑞
= 𝒖1 + 𝒖2 +
𝑾1 + 𝑾2 (2)
𝑑𝑡
2
2 2𝜏𝑠
𝑑𝑹
1
𝑹
= 𝒖1 − 𝒖2 −
𝑓
𝑹
𝑑𝑡
2𝜏𝑠
𝐿𝑚𝑎𝑥
𝑟𝑒𝑞
+
𝑾1 − 𝑾2
2𝜏𝑠
𝑓 𝑧 = 1−
𝑧 2 −1
関しては z 軸方向に 1 次元並列を行う．並列計算を
するためには常にプロセッサー間の通信が必要とな
る．1 次元並列は一番単純な方法で最大並列数が比
較的に小さいが，プロセッサー間の通信にかかる時
間が一番少ない．本研究で扱う流体場は大きくない
ので 1 次元並列を行う．ダンベル計算に間しては
(3)
流体場の領域で並列化を行い，各プロセッサーが担
当する領域間のダンベルの移動の情報を交換するよ
うにした．
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
図3：高分子の濃度とRe数による内側円筒表面での
トルク
図2：並列化の概念
5 まとめ
本研究ではテイラークエット流れにおける高分子
4 計算結果
の影響を調べた．流れに高分子を添加すると同じRe
数値シミュレーションの結果，流体に高分子を添
数でより安定な流れになる結果が得られた．そして
加することでより高いRe数で不安定が起こる結果が
添加する高分子濃度が高くなることでより安定な流
得られた．そして臨界Re数は添加する高分子の濃度
れになった．内側円筒表面でのトルクに関しては高
の増加と共に高くなる．そしてより高いRe数で流れ
分子を添加し，その濃度が高くなることでトルク減
はより不安定になりӨ方向依存性の流れに遷移する．
少効果が強く表れた．この結果はYi & Kim[3]の実験
この流れを波打つ流れという．同じ波打つ流れの状
での結果と同じ傾向を示す．しかしトルク減少効果
態でもRe数によってその波の波数と振幅が異なる．
はどのRe数でも起こるものではなく，高いRe数では
Re数が高くなることで確実な波の波数の増加が表れ
その効果がなかった．そして流体の運動エネルギー
た．高分子を添加した流れでは同じRe数で波の波数
は高分子を添加し，添加濃度を高くすることでより
または振幅がより小さくなる結果が得られた．
大きい値が得られた．今後の課題としては内側の円
添加する高分子の濃度による内側円筒表面でのト
ルクを図3に示す．図3から添加する高分子の濃度が
筒だけではなく外側の円筒も回転させる計算とより
高いRe数と高分子濃度での計算が考えられる．
高くなることで内側円筒表面でのトルクの減少率が
高くなることがわかる．しかしトルク減少はどのRe
6 参考文献
数でも起こるものではなく，Re=250では100ppmの高
[1] T. Watanabe and T. Gotoh, preprint (2012)
分子溶液でもトルク減少効果が見られなかった．
[2] Kazuyasu Sugiyama, Enrico Calzavarini, Detlef Lohse :
𝒖2
流体の運動エネルギーの平均
を次式で表す．
1
𝒖2 =
𝑢𝑟2 + 𝑢𝜃2 + 𝑢𝑧2
(5)
𝑁𝑟 𝑁𝜃 𝑁𝑧
𝑖,𝑗 ,𝑘
ここでi,j,kはそれぞれの方向の格子点番号を表す．
計算結果
𝒖2
の値はRe数が高くなることことによ
Microbubbly drag reduction in Taylor-Couette flow in the
wavy vortex regime, J.Fluid Mech. (2008)
[3] Mu-Kweon Yi, Chongyoup Kim : Experimental studies on
the Taylor instability of dilute polymer solutions,
J.Non-Newtonian Fluid Mech. (1997)
って減少した．そして添加する高分子の濃度が高く
[4] Ch Hoffmann, S Altmeyer, A Pinter, M Lucke :
なると同じRe数で増加する結果が得られた．しかし
Transitions between Taylor vortices and spirals via wavy
この値は無次元化した値であるためRe数ごとの比較
Taylor vortices and wavy spirals, New Journal of Physics
2
はスケールを合わせる必要がある．そのためにRe
をかけるとRe数増加と添加する高分子濃度の増加と
共に流体の運動エネルギーは増加した．
(2009)
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用化学分野
ナノクラスターの衝突による微細加工に関する分子動力学シミュレーション
学籍番号 22413543
指導教員名
1
はじめに
下村拓也
尾形修司教授
上述したクラスター特有の照射効果は，クラスタ
近年，多数の原子の集合体であるクラスターを入
ーを構成する多数の原子が時間的，空間的双方のス
射粒子として用いるクラスターイオンビーム技術が
ケールにおいてきわめて局所的に標的材料原子と相
新しい表面ナノ加工技術として注目を集めている．
互作用をなす(多体衝突)ことによるが，クラスター
クラスターとは，原子または分子が多数集まって
サイズ(クラスター1 個に含まれる原子・分子数)や
塊となったものである．元の原子や分子が気体であ
衝突エネルギーに対する関係は不明な部分が多い．
る場合はガスクラスターと呼び，ファン・デル・ワ
このような多体衝突現象を計算機により解析するた
ールス力により結合してクラスターを構成する．ガ
めには，系内に存在するすべての原子の運動を追跡
スクラスターは，高圧の材料ガスをノズルと呼ばれ
する分子動力学法の適用が必須である．
る細い管を通して真空中に放出することで生成され
る．高圧のガスを真空中に噴出すると，断熱膨張に
ナノクラスターの衝突を用いた微細加工として,
よりガスが凝縮温度以下まで冷却され，お互いに結
クラスター原子としてAr原子が用いられることがあ
合しクラスターが生成される．生成したクラスター
るが,照射の入射条件について、
Ar原子の分子間力に
に高速の電子を衝突させることでクラスターの電子
用いられるLennard-Jonesポテンシャルを用いた場
が弾き飛ばされることを利用して，正の電荷を持た
合と,よりクラスターの結合力の大きい, Siの共有
せてイオン化する．イオン化されたクラスターは電
結合に用いられるStillinger-weberポテンシャルを
界をかけることで簡単に加速することができ、これ
用いた場合でシミュレーションを行い,使用するポ
を標的材料に照射することで表面の加工や改質に利
テンシャルの違いによって衝突効果にどのような違
用される．
いをもたらすのかを比較することで,ナノスケール
微細加工への知見を得ることを目的とした．
クラスターイオンビームが固体表面に照射される
と，標的材料の表面は１度に何千個の原子集団の衝
2
シミュレーション方法
突を受ける。これは高密度の照射であり，極めて大
クラスターを衝突させる標的材料として,364,000
きな質量電荷比による低エネルギー照射効果，標的
原子からなる Ar 結晶と 384,000 原子からなる
材料の基板に並行な方向へのスパッタリング（ラテ
Si(100)結晶を用意した.それぞれの大きさは,前者に
ラルスパッタリング効果）など，従来の単原子イオ
ついては衝突面が 217Å 217Å,深さが 217Å,後者に
ンビームで実現が困難であった新しい照射効果が得
ついては衝突面が 244Å 244Å,深さが 244Åとなる.
られる．低エネルギー照射効果を利用すると，極浅
ただし,標的材料粒子は衝突面に並行な方向の境界
イオン注入が可能であり，
次世代LSIデバイスの接合
と隣接しており,周期境界条件により無限に広がっ
形成，無損傷半導体表面プロセスが可能となる．ラ
ていると考える.
テラルスパッタリング効果を利用すると，ナノレベ
衝突させるクラスターとしてＡｒ原子の標的材料
ルの超表面平坦化プロセス，難加工材の高速エッチ
には2000 原子からなるAr 原子クラスターを用意し,
ングが可能となる．高密度照射効果を利用すると、
入射速度に換算して 30m/s,50m/s,100m/s の運動エ
スパッタリングによる高品位薄膜の形成を可能とす
ネルギーを与えて衝突面に対して垂直に衝突させ
る．
た.Si 原子の標的材料には 2000 原子からなる Si 原
子クラスターを用意し,入射速度に換算して
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用化学分野
100m/s,500m/s,1000m/s の運動エネルギーを与え
4
て衝突面に対して垂直に衝突させた.
微細加工において，アルゴンクラスターは、壁面の
アルゴンガスの分子間力のポテンシャルにはレナ
まとめ
破壊には大きな入射エネルギーを必要とする一方で、
ード・ジョンズ(Lennard・Jones，LJ)ポテンシャ
入射するクラスター構成粒子の壁面への侵入が浅く
ルを用いた．Si 原子の共有結合には stillnger weber
水平方向に広がりやすいため，浅いイオン注入や表
ポテンシャルを用いた．シリコンに関して現在使わ
面の平坦化に適しているといえる．
れいる中では最も古い古典ポテンシャルモデルであ
る．一般的な二体間ポテンシャルに三体間の相関項
一方で、Si 原子クラスターは,Ar 原子クラスターに
を付加した形で表現されている．
比べて固いことから,クラスターの形を保ったまま
壁面内部に侵入することで，Ar 原子クラスターより
3
結果
も垂直方向に対して大きく運動エネルギーを与える
Ar クラスターの Ar 標的への衝突では,初速 100m/s
での衝突（図 1 左側）ではクラスターは壁面原子を
押しのけるように水平方向に広がり，クラスターの
衝突点から放射状に原子が飛散する様子が確認でき
た．飛散する原子速度ベクトルは、多くが衝突壁面
対して水平方向の成分を大きく持っていた．衝突後
の壁面の配位数を観測すると，衝突の影響は広く浅
い分布であった．
Si クラスターの Si 標的への初速 1000m/s での衝突
（図 4.1 右側の列）については，衝突後 2ps ほどで
壁面内部にクラスタが入り込み，その後差壁面の表
層の一部を引き剥がすようにクレーターを形成した．
飛散した原子はわずかであった.衝突後の壁面の配
位数を観測すると，衝突の影響は垂直方向に深く分
布していた．
Ar (100m/s)
0.5ps
10ps
40ps
Si (1000m/s)
2ps
10ps
40ps
図1：クラスター衝突の様子
（左列：Ar 初速100m/s , 右列：Si初速1km/s）
ことがわかった.また Si 原子はほとんど蒸発するこ
とはなかった.このことから、微細加工において,Si
原子クラスターは壁面を深く掘ることには適してい
るが，壁面内部に Si 原子が残留する可能性があると
わかった.
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
分位点回帰分析における同時学習
学籍番号 22413545 氏名 杉浦 徹
指導教員名 竹内一郎 准教授
1
はじめに
した分位点回帰分析の一般式は次式となる:
回帰分析では，入力変数と出力変数の関係を定量的に
表現した回帰モデルと呼ばれる関数を学習する．分位
点回帰分析 [1] は主に異常値検出の分野で用いられる回
帰分析の一種であり，入力変数が与えられたもとで出
力変数の条件付分位点を入力変数の関数として学習す
る．同一のデータに対して複数の分位点を推定する場
合，分位点が近い程，互いの相関が大きいと考えられ
る．しかし，従来法ではこのような相関を無視して個
別に学習するため，推定精度が劣化してしまう現象が
みられる．そこで，本研究では複数の分位点の関係を
考慮した学習を行うことで，従来法よりも精度の高い
推定を目指す．複数の類似した問題を学習する枠組み
はマルチタスク学習，もしくは同時学習と呼ばれてお
り，様々なアプローチが提案されている．本研究では，
分位点回帰分析における同時学習法として同時特徴選
択 [2, 3] と同時特徴抽出 [3] を導入し，その有用性を数
値実験にて例証する．
2
分位点回帰分析
学習データを {(xi , yi ) ∈ RD × R}N
i=1 と表す．分位点
回帰分析は分位 τ ∈ [0, 1] における条件付分位点関数
fτ : RD 7→ R を推定するアプローチであり，以下の最
小化問題として定式化される:
min
fτ
N
∑
Lτ (yi , fτ (xi )) :=
i=1
{
ψτ (r) :=
N
∑
ここで，Lt (·, ·) は τt 分位点回帰分析の損失関数 (1)，
P(·) は正則化項を表し，λ は損失項と正則化項のバラ
ンスを制御するもので正則化パラメータと呼ばれる．wt
と w0t はモデルパラメータであり，W ∈ RT ×D は t 行
目が wt となる行列，w0 ∈ RT は t 番目の要素が w0t
となるベクトルである．
3.1
同時特徴選択
学習の過程で不要な特徴を省き，少数の重要な特徴
のみで回帰モデルを構築することが重要である．この
アプローチは特徴選択と呼ばれ，ℓ1 正則化が有力であ
ると認識されつつある．
同時特徴選択 [2] は ℓ1 正則化を拡張したアプローチ
であり，全ての分位点回帰モデル fτ1 , · · · , fτT におい
て不要な特徴を省き，関連性が強い特徴のみを選択す
るアプローチである．これは W の行に対して ℓ1 正則
化，列に対して ℓ∞ 正則化，もしくは ℓ2 正則化を用い
ることで以下のように定式化される:
ℓ1,∞ 正則化 : P(W ) =
D
∑
∥wj ∥∞ ,
(2)
∥wj ∥2 .
(3)
j=1
ψτ {yi − fτ (xi )} (1)
ℓ1,2 正則化 : P(W ) =
i=1
(1 − τ )|r|, if r ≤ 0
τ |r|,
if 0 < r.
分位点回帰分析は最小二乗法と異なり，絶対誤差を用
いているため外れ値に影響されにくい．また，分位 τ
を変化させることで様々な条件付分位点関数を推定で
きる．特に，τ = 0.1, 0.9 のように極端な場合は異常値
検出などに応用することができる．
3
T ∑
N
∑
Lt (yi , w⊤
min
t xi + w0t ) + λP(W )．
W ,w0 t=1 i=1
分位点回帰分析における同時学習
本研究では，T 個の分位点 0 < τ1 < · · · < τT < 1 に関
する線形分位点回帰モデル:
fτt (x) = w⊤
t x + w0t ,
を同時に学習する問題を考える．同時学習のアプロー
チとして正則化と呼ばれる方法 [2, 3] を用いる．正則
化とは，モデルパラメータ w に制約を加えることで平
滑化を行うアプローチである．この同時学習法を導入
D
∑
j=1
ここで，wj は W の j 列目を示している．ℓ1,∞ 正則
化 (2) を用いた線形分位点回帰分析は線形計画問題と
なり，ℓ1,2 正則化 (3) を用いたものは 2 次計画問題と
なる．
3.2
同時特徴抽出
特徴抽出は特徴選択と同様に回帰モデルを少数の特
徴で構築するためのアプローチである．特徴抽出では，
もとの入力変数をそのまま使うのではなく，入出力関
係を記述するのに有用な低次元空間を抽出して新たな
入力変数とするものである．抽出されたデータが低次
元になる程，元の入力変数から不要な情報を省き，重
要な情報のみを表現したデータとなる．抽出後のデー
タが元の入力変数の線形結合で表現されると仮定する
と，抽出後のデータに対する回帰モデルは直行行列 U
を用いて次式で定義される:
ft (x) =
D
∑
j=1
atj ⟨uj , x⟩ = ⟨at , U ⊤ x⟩, ∀t.
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
T ×D
の (t, j) 要
ここで，atj は重みであり，行列 A ∈ R
素とする．同時特徴抽出では，行列 A に対して ℓ1,2 正
則化 (3) を行うことで抽出後のデータを低次元化する:
P(A) =
D
∑
∥aj ∥2 .
表 1: データセット: Size は学習/評価/テストに使用し
たデータ数，Features は真/偽の特徴数を表す．
ID
Data Name
Size
Features
D1
Normal
100/100/10000
10/10
D2
Normal
100/100/10000
10/50
D3
Laplace
100/100/10000
10/10
D4 BankSimuration 100/100/7992
32/32
D5
BirthWeight
100/100/198177 15/15
D6
RedWine
100/100/4698
11/11
(4)
j=1
ℓ1,2 正則化 (4) の影響により，A は列毎にスパースな
行列となる．つまり，正則化を強めれば抽出後の特徴
の次元数が少数になることを意味している．
4
表 2: 汎化誤差の平均値
Method
STL1
MTFS1Inf
MTFS12
MTFE12
計算機実験
各アプローチの比較実験を行う．従来法には ℓ1 正則化を
用いており (STL1)，提案法には (2)，(3)，(4) を用いた
(MTFS1inf，MTFS12，MTFE12)．データには，人工
データと UCI Machine Learning Repository より取得
した実データを用いており (表 1)，平均 0，分散 1 の正
規乱数に従うダミーの特徴を追加する．T = 5 個の分位
点 τ1 = 0.05, τ2 = 0.25, τ3 = 0.5, τ4 = 0.75, τ5 = 0.95
を同時に学習する問題を考え，正則化パラメータは λ ∈
{0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100} の
中から評価データを用いて選択する．また，人工デー
タは平均ベクトルが 0 ベクトル，分散共分散行列の
(i, j) 要素が 0.5|i−j| となる多次元正規分布から生成し，
ノイズは平均 0，分散 1 の正規分布とラプラス分布を
使用した．
実験結果を表 2 に示す．表 2 は 5 つの分位における
汎化誤差の平均値を示しており，同時学習アプローチ
は D4 の MTFE12 を除いた全ての項目において従来法
より汎化誤差が向上している．D2 において汎化性能が
向上していることから，ダミー変数を多く含むデータ
においても高い性能を示すことが分かる．また，D3 に
おいて汎化性能が向上していることから，分布の裾が
広いノイズを持ったデータに対して高い性能を示すこ
とが期待できる．各実験における分位毎の汎化誤差を
図 1 に示す．人工データ (D1-D3) において MTFS12 と
MTFE12 は分位が極端であるほど性能が向上する傾向
にあるが，実データ (D4-6) ではこの傾向は見られない．
これは，データ密度が高い分位に関する回帰モデルに
他の回帰モデルが引き寄せられていることが原因であ
ると考えられる．
次に，各データにおいて選択や抽出された特徴数を
表 3 に示す．同時特徴選択は全データにおいて STL1
より少数の特徴を選択していることから，高い特徴選
択性能を持つことが分かる．また，同時特徴抽出で抽
出された特徴数に着目すると，D4 を覗いた全データに
おいて平均 2 個以下の特徴が抽出されている．これは，
回帰モデルを構築するために必要な情報は高々2 次元
の特徴で表現できることを意味しており，データの性
質を表現するための重要な特徴が抽出できたことが分
かる．
D1
0.167
0.164
0.164
0.154
D2
0.171
0.165
0.164
0.158
D3
0.182
0.179
0.178
0.169
D4
0.079
0.075
0.077
0.148
D5
0.241
0.238
0.238
0.236
D6
0.216
0.212
0.213
0.215
表 3: 選択 (抽出) された特徴数: 上段は選択された特
徴数，下段は抽出された特徴数を表す．
Method
D1 D2 D3 D4 D5 D6
STL1
18 38.4 14.6 40.3 16.3 19.4
MTFS1Inf 12 20.8 11.3 8.6 3.8 15.5
MTFS12 12.6 22.8 12.2 17.2 7.8 16.5
MTFE12 20 60 20 64 29.7 22
MTFE12 1.6 1.5 1.5 4.1
1
1.9
図 1: τ 毎の汎化誤差: 黒点線は STL1，青線は MTFSInf，緑線は MTFS12，赤線は MTFE12 であり，縦軸
は汎化誤差，横軸は各分位 (τ1 , · · · , τ5 ) を表す．
5
まとめ
本研究では，分位点回帰分析における同時学習の導入
を提案した．今後の課題として，極端な分位における
汎化性能の向上，大規模データや非線形回帰モデルへ
の応用が挙げられる．
参考文献
[1] R Koenker., “Quantile Regression”, Cambridge
University Press, England, 2005.
[2] BA Turlach., WN Venables. and SJ Wright.,
“Simultaeous Variable Selection”, Technometrics,
2005.
[3] A Argyriou, “Learning to Integrate Data from Different Sources and Tasks”，Citeseer，2007，pp.6599.
平 成 23 年 度 創 成 シ ミ ュ レ ー シ ョ ン 工 学 専 攻 修 士 論 文 梗 概 集
計算応用科学分野
一様定常乱流の統計性における高分子の影響
学 籍 番 号 22413546 氏 名 杉 本 大 輝
指導教員 渡邊威准教授
1
乱流中に高分子を加える事で，乱流と固体面と
の摩擦が減少する摩擦低減効果はよく知られてお
り，乱流中の乱れを抑える作用がある．この効果
は，流体中の物体輸送や，配管内の流体輸送の効
率向上のために実用化されている．高分子を添加
する事で，乱流の特性にどのような影響を及ぼす
かを十分に理解する事ができれば，それらの効果
をより有効に活用できる事に繋がる．そこで，本
研究は，壁から十分に離れた空間内で高分子モデ
ルが，乱流の基本統計量及びラグランジュ統計量
に及ぼす影響を理解する事を目的とする．
解析方法は，高分子モデルが流れ場に影響を及
ぼさない場合 (1-way) と及ぶ場合 (2-way) の二通り
の 結 果 を 比 較 す る 事 で ，高 分 子 か ら の 影 響 を 考
察 す る ．高 分 子 モ デ ル は ，２ つ の ブ ラ ウ ン 粒 子
が有限伸長性を持つ非線形弾性バネで繋がれた
FENE(Finitely Extensible Nonlinear Elastic) ダンベル
モデルを扱う．計算手法は，乱流場の直接数値計
算 及 び ，大 多 数 の 高 分 子 モ デ ル の ブ ラ ウ ン 動 力
学 計 算 を 結 合 し た 計 算 を 行 う ．ま た ，乱 流 と 高
分 子 モ デ ル 集 団 の 運 動 を ，MPI(Massage Passing
Interface) 通信を用いた並列計算で計算し，その並
列化効率についても議論する．
2
2.1
fs (z) = (1 − z2 )−1 .
はじめに
数値解析法
一様定常乱流の数値計算
扱う乱流場は，一様等方性を持つ定常な乱流で
ある．一辺が 2π の立方体中で周期的に広がる空間
において，一様等方性を持つ初期速度場を作り，ポ
リマーストレステンソルによる項を加えた Navierstokes 方程式 (1) を，非圧縮条件と共に解いていく．
空間微分はスペクトル法で計算し，時間微分に関し
ては２次の Runge-Kutta 法を用いる．外力は統計的
に定常な流れをつくるため，エネルギー注入率 εin が
一定となるように a(t) を設定し，式 (2) で求める．
∇ · u = 0,
1
1
∂u
+ u · ∇u = − ∇p + ν∇2 u +
f + ∇ · T P . (1)
∂t
ρf
ρf
{
a(t)u(x，t) for kl ≤ |k| < kh ,
f (x，t) =
(2)
0
otherwise.


( (n) )
Nt  (n) (n)
3
∑

 Ri R j
νη Lbox
|R |
p
 2
T i j (x，t) =
fs
− δi j  δ(x − r(n)
g ),
τs Nt n=1
req
Lmax
(3)
(
)2
ここで，η = ΦV 3req /4a ，ΦV = 2Nt (4π/3)(a/Lbox )3 , req =
1
(kB T/k) 2 で あ る ．ま た ，rg は ダ ン ベ ル モ デ ル の 重
心 ベ ク ト ル ， R(n) は 末 端 間 ベ ク ト ル で あ る ．
2.2 FENE ダ ン ベ ル モ デ ル
図1
高分子鎖のモデル（ダンベルモデル）．高分子鎖の末端間
の弾性力をバネの弾性力で表現する．2 つのブラウン粒子は独立
してブラウン運動する．
FENE ダンベルモデルの２つの球の運動はラン
)
(
(n)
(n)
ジュバン方程式で表され， r(n)
g = x1 (t) + x2 (t) /2 と
(n)
(n)
R(n) (t) ≡ x1 (t) − x2 (t) を用いて，以下の運動方程式に
従う．
))
(
)
dr(n)
req ( (n)
1 ( ( (n) )
g
u x1 , t + u x(n)
=
W 1 + W (n)
2 ,t + √
2 (, 4)
dt
2
8τs
( (n) )
(n)
(
)
(
)
dR
1
|R | (n)
(n)
= u x(n)
f
R
1 , t − u x2 , t −
dt
2τs
Lmax
)
req ( (n)
+√
W 1 − W (n)
. (5)
2
2τs
⟩
⟨ (n) ⟩
⟨
(n)
(s)
= δmn δαβ δi j δ(t − s),
W i (t) = 0， Wα(m)
(t)W
β， j
，i
(
)
α, β = 1, 2， i, j = 1, 2, 3, n, m = 1, 2 · · · Nt ．
ここで，球の質量が十分に小さいとして慣性項を
0 と近似している．また，希薄な高分子溶液を想
定しているために高分子間の相互作用を無視し，
内部間の衝突は無いものとする．球の位置での流
(n)
体場の速度 u(xi (t)，t) は，球の周りの格子点から
3 次精度の TS(Taylor Series)13 スキーム (3) を用いて
補間して求める．時間発展は上記の式 (4)，(5) を 2
次精度の Adams-Bashforth 法を用いて解いていく．
このとき，数値積分においてバネの最大伸び長さ
Lmax を超えるのを防ぐため，|R(n) | > Lmax /2 の場合，
1 − |R(n) |/Lmax に つ い て 陰 解 法 解 き ，e(n) = R(n) /|R(n) |
は陽解法で解く．また，|R(n) | > 0.9Lmax の場合には，
ダンベルの時間刻み幅 ∆t p = ∆t/10 として，より精
度を高くする．
平 成 23 年 度 創 成 シ ミ ュ レ ー シ ョ ン 工 学 専 攻 修 士 論 文 梗 概 集
計算応用科学分野
0.6
計算条件
0.5
表 1 流体場及び高分子ダンベルモデルの計算条件
0.4
_
ε(t)
ダンベルモデルの各パラメータはポリアクリ
ル ア ミ ド を 想 定 し ，慣 性 半 径 Rg = 0.5µm，バ ネ の
最 大 伸 び 長 さ Lmax = 77µm を 用 い た ．高 分 子 モ デ
ル の 重 要 な 無 次 元 パ ラ メ ー タ で あ る Weissenberg
数 Wi = τs /τK は ，1-way で Wi = 5.0 程 度 と な る よ う
に バ ネ の 緩 和 時 間 τs を 定 め た ．そ し て ，高 分 子
と 流 体 場 の 体 積 比 ΦV = 2Nt (4π/3)(a/Lbox )3 を 設 定 す
る た め ，実 際 に 数 値 計 算 で 解 く ダ ン ベ ル モ デ ル
の 数 Ñt を b 倍 し た 数 Nt = bÑt を 本 計 算 で 扱 う 高 分
子 の 総 数 と み な す ．ま た ，ダ ン ベ ル モ デ ル の 計
算 は ，乱 流 が 統 計 的 に 定 常 な 状 態 と な っ た 時 刻
において計算領域に分散させ開始する．
ダ ン ベ ル モ デ ル の パ ラ メ ー タ 及 び ，流 体 場 の
計 算 条 件 を ，表 1 に 示 す ．本 計 算 で は ，体 積 比
ΦV = 5ppm(×10−6 )，25ppm，50ppm の三通りの計算
を行った．
0.3
0.2
1-way
0.1
動粘性： ν f
1.0 × 10−2
エネルギー散逸率： εin
エネルギー注入範囲
1 ≤ |k| < 2
Weissenberg 数： Wi
5.0
ダンベルモデルの数：Ñt
108
体積分率： ΦV
4
0.5
5ppm,25ppm,50ppm
計算結果
表 1 の計算条件において，乱流場のエネルギー
散 逸 率 ε̄(t) = 2ν⟨(∇u)2 ⟩ 及 び 運 動 エ ネ ル ギ ー ス ペ ク
ト ル E(k) を 1-way と 2-way で 比 較 し た 結 果 を ，図
1，図 2 に 示 す ．図 1 よ り ，エ ネ ル ギ ー 散 逸 率 の 低
減 が 見 ら れ る ．こ れ は ，流 体 場 の 粘 性 に よ り 散
逸 す る 運 動 エ ネ ル ギ ー の 一 部 が ，流 れ に よ り 引
き延ばされたダンベルの弾性エネルギーとなる
か ら で あ る ．ま た ，図 2 よ り ，ス ペ ク ト ル の 低 波
数 に お い て ，1-way と 比 較 し て 2-way の 場 合 は 値
が 低 く な り ，高 波 数 に お い て 値 が 高 く な っ て い
る ．従 っ て ，乱 流 に 注 入 さ れ た エ ネ ル ギ ー は 高
分 子 に よ り 輸 送 す る た め ，高 分 子 添 加 に よ り 低
波 数 領 域 で の エ ネ ル ギ ー が 減 少 し ，高 波 数 領 域
でのエネルギーが増加される事が分かる．
120
130
t
140
150
160
170
180
102
101
100
10-1
10-2
10-3
10-5
128
4.0 × 10−3
110
103
3
時間刻み幅： ∆t
100
図 1 エネルギー散逸率 ε̄(t) = 2ν f ⟨(∇u)2 ⟩ の時間推移．t =
100 で Ñt = 108 個の FENE ダンベルモデルを計算領域に一様に
分散させ，ΦV =5ppm,25ppm,50ppm と変化させた時の様子．
10
格子点数： N
2-way(ΦV=25ppm)
2-way(ΦV=50ppm)
-4
3
2-way(ΦV= 5ppm)
0
90
____ _
ln (E(k,t)t ε-2/3 lK-5/3)
3
1-way
2-way(ΦV= 5ppm)
2-way(ΦV=25ppm)
2-way(ΦV=50ppm)
0.1
ln (k lK)
1
図 2 エ ネ ル ギ ー ス ペ ク ト ル E(k, t) の 時 間 平 均 値（ 縦 軸：
t
t
t
(ε̄(t) )2/3 (lK (t) )5/3 ，横軸：lK (t) で無次元化）．高波数側は，エ
√
イリアシング誤差除去のため，2 2N/3kmax ≤ k で E(k, t) = 0 で
ある．
5 まとめ
本 計 算 に よ り ，一 様 定 常 乱 流 中 に 高 分 子 モ デ
ル を 加 え る 事 で ，乱 流 の エ ネ ル ギ ー 散 逸 低 減 効
果や加速度のラグランジュ時間相関が減少する
事 を 確 認 し た ．ま た ，実 験 結 果 と 比 較 す る 事 で
結 果 の 妥 当 性 を 示 し ，本 計 算 に 用 い た プ ロ グ ラ
ム の 並 列 計 算 の 効 率 を 検 討 す る 事 が で き た ．今
後 は ，計 算 す る 高 分 子 モ デ ル の 数 の 増 加 や ダ ン
ベルモデルを鎖モデルに拡張するなど物理的効
果が本研究結果にどのような影響を及ぼすか検
討 し て い く.
参考文献
1) T. Watanabe and T. Gotoh， Phys. Rev. E 81，
066301(2010)
2) A. M. Crawford， N. Mordant， H. Xu and E. Bodenschatz，New J. Phys. 10，123015(2008)
3) P. K. Yeung and S. B. Pope， J. Comp. Phys. 79，373-
416(1988)
4) T. Watanabe and T. Gotoh, preprint(2012)
平成 23 年度 創成シミュレーション工学専攻 修士論文梗概集
計算応用科学分野
乱流中における雲粒子の運動
学籍番号：22413548 氏名：鈴木 祐有紀
指導教員名：後藤 俊幸 教授
1. はじめに
微分にはフーリエスペクトル法を，時間積分には連
地球環境問題や気象予測において雲は大変重要な
続体部分，
粒子部分共に Runge-Kutta-Gill 法を用い
要素であるため，雲の構造や発生，成長過程の解明
た．1 辺 25.6cm の立方体を計算領域とし，3 軸方向
は重要かつ大きな課題である．雲の動きには数10km
に周期境界条件を課して計算を実行した．
から数μmスケールまで渡る様々な物理過程が存在
し，観測や実験上の制約からこれまでに行われた実
3. 3 次元 FFT の改良
験，数値計算で得られたデータだけでは決して充分
NS 方程式の非線形項の評価にFFT を用いると計算
でない．本研究では直接数値計算(DNS)を用い，一様
のほとんどが FFT に費やされるため，FFT の高速化
等方乱流により輸送される半径10μm程度の雲粒子
は非常に重要である．我々は今まで図 1 左のように
の挙動を解析する．さらに，Navier-Stokes(NS)方程
計算領域を 1 つの軸方向に分割して FFT を計算して
式の非線形項の評価に用いる高速フーリエ変換
いたが，本計算では図 1 右のように計算領域を 2 軸
(FFT)を効率良く計算する方法についても触れる．
方向に分割し，それぞれのブロックに 1 つのプロセ
スを割り当てた FFT コードを開発した．並列化して
2. 構成方程式と計算手法
FFT を行う場合，座標の転置は不可避である．2 次元
本研究では地上温度6℃，高度1500mに存在する雲
領域分割の場合，4 回の転置が必要であるが，それ
を温度，浮力，水蒸気混合比などの基本的な物理過
を 2 回に減少させることにより，3 割程度の高速化
程を考慮して数値解析を行う．乱流場は非圧縮で一
につながった．テストとして Fourier 空間内の全格
様等方かつ定常状態であると仮定し，以下の方程式
子点でランダムに与えられた乱流のフーリエ係数に
を考える．乱流場，温度場，水蒸気混合比はEuler
対して逆 FFT と FFT を 1 回ずつ行い，変換に要した
的に，粒子はLagrange的に方程式を解いていく．
時間を計測した．格子点数 N
3
=20483 の場合におけ
る，ランク 0 での 1 回当りの FFT または逆 FFT に要
した時間の平均を図 2 に示す．
図1：1次元領域分割（左）と2次元領域分割（右）
また，格子点座標を使って粒子位置における物理量
を補間する際には，Taylor Series 13 (TS13)スキー
ム[1] を用いた．そして，連続体方程式における空間
図 2：1 回あたりの FFT の計算時間(N 3=20483)
平成 23 年度 創成シミュレーション工学専攻 修士論文梗概集
計算応用科学分野
図 2 中の M, Npx はそれぞれ全プロセス数，第 1 軸
方向の分割数を表す．図 2 から，第 1 軸方向の要素
数=N/Npx=512 にすると変換に要する時間が急激に
落ちていることが分かる．この傾向はほかの格子点
数の時にも確認できた．最適な分割数を考えるため
の端緒になるといえよう．また，開発したスキーム
を用いて N
3
=20483 の高解像度パッシヴスカラー乱
流の計算を行ったところ，1 次元領域分割の時に比
べ，4 割程度の計算時間を減少させることができた．
図4：雲粒子半径に関するPDFの時間変化
4. 計算結果
ここでは結果の一例として，雲粒子の分布と雲粒
成は起こらないようにしているため，雲粒子と未飽
和空気が触れ合うことにより，粒子半径は減少して
子半径の確率密度関数(PDF)の時間推移について述
いくのみである．また，雲粒子周りの空気の状態は
べる．粒子分布の変化に関して，雲粒子の初期位置
時々刻々変化するためにPDFの分布が非対称となり，
を図3上に，
1.0秒後の粒子位置を図3下にそれぞれ示
時間経過に伴ってその裾野を広げていく．
す．この結果から．雲粒子を入れてから1.0秒後には
空間全体に一様に広がっているため，乱流により雲
5. おわりに
粒子が急速に混合されている様子が分かる．また，
本研究では，一様等方乱流中に発生させた雲粒子
重力による効果は乱流速度揺らぎが重力による粒子
の挙動を解析した．雲粒子が乱流によって急速に混
の終端速度よりもはるかに大きいために見られない．
合され，空間全体に拡散する様子や，乾燥空気によ
次に，雲粒子半径のPDFの時間変化を図4に示す．図4
って雲粒子の半径が減少することを確認することが
から，粒子半径は減少傾向にあることが分かる．本
でき，大方予想通りの結果が得られた．そのため，
計算では，水蒸気が凝結することによる雲粒子の生
近い将来に行う，大規模数値計算をするための準備
ができた．さらに，2次元領域分割によるFFTコード
の開発により，格子点数を増大させた時に効率良く
計算する方法を見出せた．今後は解像度と雲粒子数
を上昇させた計算の実行と，粒子衝突や併合などを
考慮したより正確なモデルを構築する方針である．
参考文献
[1] P. K. Yeung, S. B. Pope，“An Algorithm for Tracking
Fluid Particles in Numerical Simulations of Homogeneous
Turbulence”，Journal of computational physics 79,
373-416 (1988)
[2] Bipin Kumar, Jörg Schumacher, Raymond A.Shaw
“Cloud microphysical effects of turbulent mixing and
entrainment”(To be submitted)
図3：雲粒子の位置（上：初期状態，下：1s後）
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
乱流スカラー輸送ハイブリッドコードの開発
学籍番号
22413561 氏名
畑中
祥吾
指導教員
後藤
俊幸
a1 ( f i ′+1 − f i ′−1 ) + h f i ′′+ b1 h ( f i ′+′1 − f i ′−′1 )
1 緒言
乱流場によって混合されるパッシブスカラー輸送
問題は，伝熱や物質の拡散をはじめ多くの輸送現象
=
c0
f + f i −1
f + f i−2
+ c2 i+2
f i + c1 i +1
h
h
h
に適用可能な基礎的問題であり，その理論の発展は
を連立させて解く事により微分係数を評価する陰的
広範な分野への貢献が期待される重要な問題である．
スキームである[2][3]．微分演算がフーリエ空間で
パッシブスカラーDNSの代表的手法であるスペクト
波数 ik の乗算として表わされるのに対応して，
差分
ル法は高い計算精度を持つ一方，
FFTに必要なノード
スキームを用いた場合の ik に対応する波数は
間多対多通信により大規模並列計算ではその実効性
Modified Wavenumber と呼ばれ，スキーム精度評価
能が大きく制限されている．
の重要な指標である(図 1)．
本研究では，パッシブスカラー輸送問題のDNS に
対する新たな手法の提案を目的とし，速度場の計算
をスペクトル法を用いて行う一方で方程式が空間的
に局所的であるスカラー場の計算には結合コンパク
ト(Combined CompactDifference:CCD) スキームに
よる高精度差分法を用いるというハイブリッドシミ
ュレーションを新たに開発し，その有効性の検証を
行った．
図 1：1 階微分に対する Modified Wavenumber
2 流体力学の基礎理論
非圧縮粘性流体の基礎方程式である連続の式，NS
方程式及びパッシブスカラー輸送方程式は以下で記
述される
r
∇ ⋅u = 0
r
r
r
r
∂u
+ (u ⋅ ∇)u = −∇p + υΔu
∂t
r
∂θ
+ (u ⋅ ∇)θ = Δκθ
∂t
ハイブリッドシミュレーションでは速度場をスペ
クトル法，
スカラー場を差分法により計算しており，
時間積分は全て4次Runge-Kutta-Gill法により行っ
た．差分法には 8 次精度 CCD スキーム(CCD8)，6 次
精度 CCD スキーム(CCD6)，4 次精度中心差分(FD4)を
用いたシミュレーションをそれぞれ実行し，計算結
果の比較を行った．
ここで，νは動粘性係数，κは拡散係数である．
4 CCD スキームの MPI 並列化コード開発
ハイブリッドシミュレーションを行う1 つの理由
3 数値計算法
諸言にも述べた通り，本研究ではスペクトル法と
は，
スペクトル法に用いるFFTは通信量が多く大規模
CCD スキームによるハイブリッドシミュレーション
並列化が難しい事であった．その点を検証するため
を行っている．CCD スキームとは陽的スキームであ
CCD スキームとFFT のMPI 並列化効率を調べた(図
る中心差分などとは異なり，各格子点上で方程式
2)．図2より，空間を3次元的に分割し並列化したCCD
a1 ( f i′+1 − f i′−1 ) + f i ′ + b1 h ( f i ′+′1 − f i′−′1 )
= c1
f i +1 − f i −1
f − f i−2
+ c2 i+ 2
h
h
スキームが最も良いスケーラビリティを示しており，
並列化効率100%に近い性能が出ている事が分かる．
一方で，
2 次分割のCCDスキームとFFT は並列化のス
ケーラビリティに関しては同程度であった．
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
5 結言
本研究ではパッシブスカラー輸送問題のDNS に対
する新たな手法として，速度場をスペクトル法で解
きスカラー場を結合コンパクトスキームで解くとい
うハイブリッドシミュレーションを試みた．
CCDスキ
ームのプログラム開発においては，高い並列化効率
を持つ3 次元分割のCCD スキームの開発を行うこと
が出来た．それにより，ハイブリッドシミュレーシ
図2：CCDスキームとFFT のMPI 並列化特性. 太線(右軸) は並
ョンでは計算時間をスペクトル法の7割程度に抑え
列化のスケーラビリティで黒線が並列化効率100%の曲線, 細
る事ができた．計算精度の検証では，スペクトル法
線(左軸) は演算性能.
とCCDスキームはほぼ同等のシミュレーションが行
える事が確認できた．特に，8 次精度CCDスキームよ
5 ハイブリッドシミュレーション
り精度の劣る6 次精度CCDスキームでもスペクトル
格子数1024のスペクトル法及びハイブリッド法で
法と同等の結果が得られた事は，
8 次精度CCDスキー
一様等方減衰BOX乱流のDNSを実行し計算時間及び計
ムを用いることで十分な精度のDNS が行える事を強
算結果を比較した．ハイブリッド法では計算時間が
く示唆しているものと考えられる．
スペクトル法の7割程度に減少しており，
スカラー場
今後の展望としては，まずハイブリッドシミュレ
の計算に限れば計算時間は半分程度であった．次に
ーションをスカラー場の記述がより複雑になる実際
格子数2048のスペクトル法による結果を基準場とし
的な物理系に適用すれば，方程式を実空間で平易に
て各スキームでの計算結果を比較した(図3,4)．図3
扱いながら，スペクトル法と同等の精度の高効率シ
のスペクトルでは各スキームとも奇麗に一致してい
ミュレーションを行える事が期待できるだろう．ま
る．一方で図4のスカラー微分量のPDFを見るとFD4
た，今後のペタスケールコンピューティングへの発
でのみPDFの劣化が確認されたが，スペクトル法と
展性を考えると，NS 方程式まで差分法で解いた，
CCDスキームでは同等の結果が得られた．
より高効率のシミュレーションを行う事が考えられ
る．そのようなシミュレーションをスペクトル法
DNS の代替とするためには，ポアソン方程式のみ
FFT で解くなど質量保存則を高精度で満たす特別
の取り扱いが必要になるだろう．
6 参考文献
[1] S.K.Lele : Compact finite difference schemes with
spectral-like resolution,
図 3：スカラー分散スペクトル
J.Comput.Phys. 103,(1992)16-42
[2] P.C.Chu, & C.Fan : A three-point combined compact
difference scheme,J.Comput.Phys. 140,(1992)370-399
[3] K.Mahesh : A family of high order finite difference
schemes with good spectral resolution,J.Comput.Phys.
145,(1998)332-358
[4] Nihei T., & Ishii K.: Parallelization of a highly
accurate finite difference
scheme for fluid flow calculations,Theor.Appl.Mech.Jap.
図 4：スカラー1 階微分の PDF
52(2003)7181
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
第一原理計算による水分子がシリカガラスの破壊に及ぼす影響の考察
学籍番号 22413572
1
はじめに
氏名
水越智彦
指導教員名
尾形修司
影響を，局所構造と機械特性の相関から明らかにす
ガラスは生活に密接した材料である．特にシリカ
にすることを目的とする．
(SiO 2 )ガラスは光ファイバーなどに用いられてお
り，現代産業における重要な材料である．
2
理論的には完全結晶と同程度の，Si-O結合強度を
€
シミュレーション手法
本研究では，密度汎関数理論(Density Function
持つと考えられる．しかし，実際のシリカガラス材
Theory: DFT)に基づいて，第一原理引張試験を行う．
料は脆く，
予想される理論強度の1/100程度の強度で
計算には第一原理 PAW(Projector Agumented-wave)
破壊してしまう．その原因は，ガラス内部や表面に
法計算ソフト QMAS[2]を用いた．
存在するき裂や空孔などの欠陥の先端に外力が集中
本研究で用いたスーパーセルを Fig.1 に示す．こ
の系は SiO 2 ユニットを 32 個，合計 96 原子を含んだ
し，理論強度に達するからと考えられている．
ガラスは結晶と異なり，原子配列に周期性がない
一辺 11.32Å の立方体である各境界面には周期境界
ため，
X線解析などの実験からは原子レベルの局所構
条件を課している．原子配置については，アモルフ
造解析は困難である．そのため，原子構造を再現で
ァス状態を実現するために，分子動力学(Moluculer
€
きるシミュレーションが有効的である．Si−O結合の
Dynamics: MD)計算で高温融解させたシリカを急冷
切断など，電子が支配的な現象には量子論に基づく
したものを採用している[3]．また，この系には三員
第一原理計算によるアプローチが必要である．
環が 2 つ含まれている．
シリカガラスはSiO 4 四面体を基本構造とし，頂点
を共有してランダムなネットワークを形成している．
特徴的な構造に，SiO 4 四面体3つでリングを作る，
三員環構造がある．三員環は自由度が少なく歪んだ
€
構造であり，機械・光学特性に影響する．
また不純物を含むことでその性質を変化させる．
€
特にH 2 O分子は強度低下の原因として知られている．
最近の計算によると，H 2 O分子は三員環と活性が高
く[1]，その歪みを解放すると予測されている．
€
シリカガラスとH 2 O分子の反応は次のようである．
≡Si-O-Si≡
+ H 2 O → 2≡Si-OH
€
この反応に従い，シリカガラスの破壊は進展する．
Figure 1: 引張試験に用いたSiO 2 ガラスの
しかしながら，実験からはH
2 O分子がどのように局
€
スーパーセル(黒球:Si，白球：O)
所構造に作用するかはわからず，現在も解明されて
€
いない．
引張試験では，
スーパーセルのxyz各軸の一方向の
みセルサイズを伸ばしていく．各歪み値で構造緩和
€
本研究では第一原理計算によりシリカガラスの引
€
を行うので，非常にゆっくりとした引張速度の計算
張試験を行う．第一原理引張試験はガラスなどのア
に相当する．なお本研究では，不純物を含まない場
モルファスを扱った例はまだ少ない．
合とH 2 O分子，OH 分子をそれぞれを含んだ場合の3
−
本研究は，シリカガラスの引張に対する応答や，
パターンについて行った．
−
H 2 O分子やOH 分子がシリカガラスの破壊に及ぼす
€
€
€
€
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
3
第一原理引張試験の結果
−
・ OH 分子を含む場合
−
・ 不純物を含まない場合
OH 分子を三員環の近傍，遠方に配置した場合
の7パターンについて引張試験を行った．Figure 4
各軸の引張試験の歪みに対する応力の変化を
Fig.2 に示す．応力が低下しているとこではボン
€
にその内の3パターンの歪みに対する応力の変化
ドが破断している．
x,z 軸とy 軸では傾きが違い， €
を，不純物を含まない場合と比較する．三員環近
その挙動に違いが見られる．
傍にOH 分子を配置するとボンドが早く破断した．
−
−
SiO 4 四面体同士のねじれ角の平均値をFig.3に
三員環とOH 分子が影響し，歪みを解放しやすく
なったと考えられる．
示す．y 軸引張の方がよりねじれている．原子構
造を観察すると，x,z 軸引張では三員環やその周
25
辺のボンドが，y 軸引張では三員環とは離れたボ
€
Stress (GPa)
€
€
ンドが破断していた．SiO 4 四面体のねじれにより
歪みを解放しているが，三員環はそのねじれを阻
害する傾向が見られた．
€
20
15
10
5
pure
near 3-mem-ring
far 3-mem-ring
near 3-mem-ring
0
-5
0
5
10
15
20
Strain (%)
−
Figure 4: OH 分子を含む場合の
歪みに対する応力の変化
4
まとめ €
シリカガラスを対象にした第一原理計算は数多
くあるが，本計算のように破壊のダイナミクスを扱
−
った例はまだ少ない．OH 分子と三員環は他の構造
に比べ反応しやすい傾向が見られた．三員環は歪ん
だ構造であり，シリカガラスの機械・光学特性に影
Figure 2: 歪に対する応力の変化
響することが分かっている．その構造を破壊する傾
€
−
向があるOH 分子はシリカガラスの特性を変化させ
る要因となると考えられる．
しかし，より現実的なモデルに近づけるためには
まだ課題がある．特に化学反応が再現でき，さらに
€
原子のダイナミクスを同時に扱えることが出来れば，
計算機シミュレーションはシリカガラスの破壊過程
の解明の手助けになると予想できる．
Figure 3: SiO 4 四面体同士のねじれ
5
参考文献
[1] G. Batyray, et al., Phys. Rev. Lett.,
100, 105503
・ H 2 O 分子を含む場合
H 2 O 分子を含んだ場合では，挙動に変化が見ら
€
€
€
€
(2008).
[2] http://qmas.jp
れなかった．原因として０K でのシミュレーショ
[3] T. Tamura, G. ‒H. Lu, R. Yamamoto and M. Kohyama, Phys.
ンだったため，活性化エネルギーを超えられず，
Rev. B, 69, 195204 2004)
(
H 2 O 分子と SiO 2 の化学反応を再現出来なかった．
€
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
乱流により輸送されるパッシヴスカラーの確率密度関数
学籍番号 22413578 氏名 和田 恭志
指導教員名
1 はじめに
後藤
俊幸
ると
現在我々は，大気汚染や水質汚染など地球規模で

2
  k ui k , t   Pij k N j k 
 t

環境問題を抱えている．乱流により輸送される温度
や濃度などのスカラー量の振る舞いを解明すること
は，これらを解決のための大きな手助けとなる．乱
流場の性質を知る上では，統計的な平均に注目する
のが有効であり，本研究では，乱流の統計的な性質
を表す量として代表的な確率密度関数(Probability

2
  k  k , t   ik  F  x, t ux, t 
 t

と書ける．この式を数値的に解いていく．
Density Function，以下 PDF)に注目する．この PDF
に注目した研究は過去にも多く行われてきたが，乱
3 計算方法
流中のスカラー揺らぎの PDF は速度の PDF ほど理解
本研究では，一辺 2 の立方体の一様等方性乱流
されてはおらず，PDF がガウス型になるという報告
を考え，周期境界条件を適用した．そして，非線形
[1]もあれば，PDF が指数型になるという報告[2]も
項を高速フーリエ変換(FFT)を用いたスペクトル法
ある．このように，多くの文献でパッシヴスカラー
で解き，
時間積分に４次の Runge-Kutta-Gill 法を用
の PDF の関数形が異なったタイプになる条件を説明，
いて計算した．
分類しようと試みているが，基礎的な法則に基づい
初期条件は，以下のエネルギースペクト，スカラ
た説明はまだなされていない．そこで本研究では，
ー分散スペクトルに従う正規乱数によるランダムな
速度場の積分長とシステムサイズの大小関係が PDF
場を与えた．
の違いを説明する一つの要因ではないかと推測し，
乱流輸送とそのスカラー場の DNS によって，スカラ
ーの PDF の振る舞いを解明することを試みた．
2 乱流場を支配する基礎方程式
非圧縮流体の速度場は連続の式と Navier-Stokes
方程式，スカラー場は移流拡散方程式

E (k ,0)  16 2  u02 k 05 k 4 exp  2k kv 0 
2


3 2
u0   E (k ,0)dk　，u0  1
0
2
2
E (k ,0)  16 3 2  k 05 k 4 exp  2k k 0 



1 2
 0   E (k ,0)dk
0
2
u  0
ここで， kv0,kθ0 はそれぞれ速度とスカラーの初期
u
1
 u   u   P   2u
t

スペクトルの最大値の波数を表している．エネルギ

 u     2
t
によって記述される．乱流の変動速度を波数空間で
ースペクトルは各波数に対応するスケールにどの程
度エネルギーが分配されているか表すので，この
kv0,kθ0 値を変え，支配的なスケールの揺らぎと領域
サイズの影響が PDF にどのような影響があるか検討
した．
取り扱うようにすれば，乱れの運動エネルギーのカ
スケードの仕組みや，
渦相互間のエネルギー輸送や，
4 計算結果
粘性によるエネルギー散逸などの詳細を理解するの
ゆらぎのスケールと領域サイズがPDFにどのような
に都合がよい．物理空間 x での NS 方程式と移流拡
影響を与えるか，他のパラメータは統一し，初期ス
散方程式をフーリエ変換して波数空間 k で表現す
ペクトルの最大値の波数kv0=4,16，
kθ0=4,16で組み合
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
わせを変え４つの条件で計算した．その結果，速度
場に関わらずスカラー分散スペクトルをkθ0=4とし
たときはガウス型(図1)，kθ0=16としたときは指数型
(図2)のPDFが確認された．kθ0=16の条件のスカラー
分散スペクトル(図3)が示すように，
高波数側(k>10)
でスカラー分散が大きい値をとるとき，すなわち，
領域全体に対して小さなスケールの揺らぎが支配的
であるときにスカラーのPDFは指数型に近くなると
いうことが言える．ScやReを変化させた計算も行っ
たが，
結局のところスカラー揺らぎのスケールがPDF
に影響しているということがわかった．
スカラーのPDFについてどのような特性なり表現
が可能であるか探ってみたとき， を与えたときの
Δθ,(∇θ)2の条件付き平均をそれぞれ<Δθ|θ>,
<(∇θ)2|θ>とすれば
図4：条件付き平均の指数型の例
 | 
 2 | 
 | 
 2 | 

1

2
2 

2 
 2 

 ならば P( )  exp  
 const. ならば P( )  exp  |  |
ということが解析的に導けた．図4に示すように，実
際にこの条件付き平均値を計算したところ上記の理
論が正当であるということを確かめることができた．
5 まとめ
過去の研究において，
スカラー揺らぎのPDFは速度の
PDFほど解明はされておらず，PDFがガウス型になる
図1：ガウス型のPDF(kv0=4,kθ0=4,kvou0t=3.2)
という報告と指数型になるという報告があった．
本研究では，積分長とシステムサイズの大小がその
違いの一つの要因ではないかと推測し，三次元等方
乱流のDNSによって，スカラー揺らぎのPDFの振る舞
いを解明することを試みた．その結果，スカラー場
の高波数域$k>10$にスカラー分散が多いとき，すな
わち，空間に対して小さなスケールの揺らぎが支配
的であるときにスカラーのPDFは指数型になるとい
う結果が得られた．しかし，この結果について十分
図2：指数型のPDF(kv0=4,kθ0=16,kvou0t=4.0)
な物理的な解析を付け加えることができず，それを
今後の課題としていきたい．
参考文献
[1] Watanabe, T., and Gotoh, T. 2004. Statistics of a
passive scalar in homogeneous turbulence. New J. Phys.,
6, 40.
[2] Emily, S.C.Ching., and Yuhai, Tu. 1994. Passive scalar
fluctuations with and without a mean gradient: A numerical
図3：スカラー分散スペクトル(kv0=4,kθ0=16)
study. Phys. Rev. E., 49, 1278-1282
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用分野
乱流中の固体粒子群の沈降速度に関する数値的研究
学籍番号
21413515
氏名
指導教員名
宇野聡志
渡邊威
1 はじめに
ような周囲 8 点の格子点での流体速度から線形補間
本研究の目的は,この乱流中の粒子の振る舞い,相互
して求める. 粒子の初期位置は計算領域において一
作用について調べるものである.乱流中の粒子の分
様に分布するように乱数を用いて配置し,それぞれ
散や,粒子の運動による乱流への影響だけでなく,特
の位置での流体の速度を初期速度とした.本研究で
に重力下で乱流が粒子の沈降速度に及ぼす影響を調
は乱流場と固体粒子の相互作用を考える時,1-way
べる事に焦点をあてる.静止流体中の自由落下する
と 2-way の二つの方法を考える.1-way では,粒子は
粒子の沈降速度に比べて乱流中の粒子の沈降速度は
周りの流体の影響を受けて運動するが,流体には影
大きくなる傾向がある事が,これまでの数値実験に
響を及ぼさないと仮定した計算を行う.つまり,粒子
おいてわかっている.しかしこれまでの研究では粒
の運動を計算する際に流体の速度場のデータを用い
子の運動の計算において,計算コストの問題からす
るが,流体場を計算する際には粒子からの反作用力
べての粒子について計算するのではなく,複数の粒
を0 とする.一方,2-way では,粒子は周りの流体の影
子を一つの粒子に代表させて計算結果を得てきた.
響を受けて運動するだけでなく,その運動によって
そこで本研究では,並列計算を用いる事ですべての
流体への影響を及ぼすものとして扱う.流体の運動
粒子に関して個別に運動方程式を解く事を目指した.
を計算する際は,粒子との運動量の授受を考慮する.
また,乱流と粒子の相互作用を計算するにあたって,
2-way で計算する際に,粒子から乱流場への影響を
両者をつなぐ方法についてOne-way coupling(以下
ストークス抵抗力の反作用力によって与える.
1-way)とTwo-way coupling(以下2-way)の2つを用い,
その違いによる結果の類似点と相違点について調べ
る.
一様定常乱流を生成する為に,フーリエ変換された
2 乱流場と粒子の基礎理論
速度場の低波数側に負の拡散項に類似した決定論的
流れは非圧縮流体の Navier-Stokes 方程式により支
な外力を付加する.流れの計算では,この外力と粒子
配されると仮定する.
からの反作用力を Navier-Stokes 方程式に加えて,4
次のルンゲクッタジル法を用いて微分方程式を解く.
また各時間ステップ毎に各方向の平均速度が 0 とな
るように速度場を修正する. また,乱流中の多数の
また, すべての固体粒子は粒子間の相互作用が無視
粒子を計算する為に Message Passing Interface
できるほど濃度が低く,流体運動のスケールよりも
(MPI)を用いた並列計算を行う. 乱流場の計算に 1
ずっと小さいようなサイズを考える.即ち粒子を質
つのノードを用い,粒子の運動の計算には複数のノ
点として扱う.また,固体粒子に働く力はストークス
ードを用いる.乱流場計算ノードからは速度場デー
の抵抗則からくる抵抗力と重力のみを考える.
タを,各粒子計算ノードからは反作用力を渡す. こ
れらの並列計算を行う際に,各粒子計算ノードと乱
流場計算ノードの計算量に偏りが出ないように粒子
数,ノード数を調整する事が望ましい.
3 数値計算手法
4 計算結果
乱流中の固体粒子の時間発展を求める時,その粒子
2-way で計算を行った場合,乱流場の統計量のうち
が存在する点上での流体の速度を,その粒子を囲む
運動エネルギー,エネルギー散逸率が 1-way に比べ
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
て減少し,マイクロスケールのレイノルズ数,
計算応用分野
にし反作用力を 10 倍にして返すような計算も行い,
Skewness は 1-way でも 2-way でも大きな変化は表れ
すべての粒子を計算した結果との比較を行ったとこ
なかった.固体粒子群の速度のアンサンブル平均を
ろ,乱流場の統計量,固体粒子群の統計量共に大きな
とり,静止流体中を自由落下する粒子の終端速度で
変化は見られなかった.
割ったものの時間発展を 1-way と 2-way それぞれ図
1,図 2 に示す.
図3：計算粒子を用いた粒子の平均速度（2-way）
5 まとめ
図1：粒子の平均速度（1-way）
本研究では,重力下で乱流中に多数の粒子を分散さ
せた時,どのような影響を及ぼすのかcoupling方法
を変えて調べた.その結果,2-wayにcoupling方法を
変えた時,固体粒子群から運動エネルギーを奪い,沈
降速度がより増加する事がわかった.また,計算コス
トを抑える為に1/10程度の計算粒子を用いる事は計
算結果を損なわずに計算コストを抑えるのに妥当で
あるといえる.
6 参考文献
図2：粒子の平均速度（2-way）
[1] 後藤俊幸“乱流理論の基礎”朝倉書店
どちらの計算結果も沈降速度が終端速度に比べて増
[2] 木田重雄 柳瀬眞一郎“乱流力学”朝倉書店
加しており,更に 2-way で計算された沈降速度は 1
[3] Thorsten Bosse and Leonhard Kleiser: “Small
-way で計算されたものに比べてより増加している
particles in homogeneous turbulence: Settling velocity
事がわかった. 重力下での乱流中の粒子は,渦構造
enhancement by two-way coupling”, Phys. Fluids, Vol. 18,
にぶつかる時に渦の周囲を下向きの流れに沿って移
227102, (2006)
動していく傾向があることが知られている.この選
[4] L-P Wang and M.R. Maxey: “Settling velocity and
択的掃過(preferential sweeping)と呼ばれる現象
によって乱流中の粒子は乱流場の速度ܷ３が負の値
concentration distribution of heavy particles in
となるような領域に集中すれば,沈降速度の大きな
256, 27, (1993)
増加をもたらす.2-way の場合では,こうした現象か
[5] 鈴木敏夫“固体粒子の乱流拡散とその乱流への影響” 平
ら沈降速度の増加した固体粒子群が流体を重力方向
成 19 年 修士論文
に押し出すような影響を与え,その結果として,更に
選択的掃過による固体粒子の沈降速度の増加が起こ
るのではないかと考えられる.
また,10 個の粒子を代表するような1 個の粒子(計算
粒子)を用いた計算,つまり計算する粒子数を 1/10
homogeneous isotropic turbulence”,J . Fluid Mech, Vol.
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
ランキングを考慮した推薦システムのための
行列補完アルゴリズム
学籍番号 21413546 氏名 棚橋 優樹
指導教員名 竹内一郎 准教授
1
はじめに
近年の Web の発展に伴い生まれた，amazon.com に
代表されるネット通販ビジネスには，ユーザーが求め
るコンテンツを，それまでのユーザーの履歴等と関連
の深いものから選択し提示する推薦システムというも
のが必要不可欠なツールとなっている．推薦システムの
構築は履歴データに基づいて行われるため，機械学習
を用いたアプローチが有効とされている．本研究では
機械学習による推薦システム構築法のうち，協調フィル
タリングを用いた行列補完の問題を考察する．協調フィ
ルタリングとは前もってユーザーからコンテンツに対
する評価データを収集しておき，そのデータに基づい
て類似した評価を行っているユーザーを近似し、それら
の評価情報に基づいて行列を補完する手法である．従
来の方法では，実際のニーズと推薦結果との間にギャッ
プが存在しているため，本研究では行列の各要素に重
みを設定し，そのギャップを埋める試みを検討する．
2
推薦システム
推薦システムには大きく分けて二つの種類がある．
ユーザーの行動履歴中のコンテンツの内容に基づくフィ
ルタリング方法と，コンテンツの内容は見ずに，前もっ
て多数のユーザーのコンテンツに対する評価を収集し
ておき，そのデータに基づく協調フィルタリングとい
う手法である．本研究では協調フィルタリングを扱う．
協調フィルタリング手法で解く問題を行列補完問題と
いう．収集されているデータは，ユーザー m 人の情報
を行方向に，コンテンツ n 個の情報を列方向に格納し
ている行列である．しかし，コンテンツは無数にある
のに対し，ユーザーがレビューしているコンテンツの
数は行列の全体の要素のうち, 数% に過ぎない. 行列
補完問題とは，そのような行列の値が解らない要素の
値を，解っている要素の値から予測し補完する問題で
ある. 協調フィルタリングの手法はユーザー間の類似
度から予測値を計算する User-based 方式 [3]，コンテ
ンツ間の類似度から予測値を計算する Item-based[4] 方
式等があり，本研究では Low-Rank Approximation[1]
という手法を用いる．
3
Low-Rank Approximation
Low-Rank Approximation(LRA) とは，行列の階数
（ランク）を小さくすることで少ない特徴量で行列を
表現し，予測行列の近似を行う手法である．しかし，
通常はランク最小化は定式化しにくいため，特異値の
合計であるフロベニウスノルム最小化に近似する．フ
ロベニウスノルム最小化は特異値分解（Singular Value
Decomposition）を用いて解を導く．入力行列を Y ，求
める行列を X とし，値が観測されているインデックス
の集合を Ω とおき，過学習を抑えるための正則化項を
加えた SVD の主問題は
min
(Y − X)2Ω + γ∥X∥∗
(1)
と表される．学習は求める行列 X とそれを分解し定義
される K を交互に更新する EM アルゴリズムを利用
する．ここで，K とは，X = U V , K = V V T で表さ
れる行列である．
4
重み付き Low-Rank Approximation
Low-Rank Approximation を重み付き問題に拡張し
たものは，
min
W Ω (Y − X)2Ω + γ∥X∥∗
(2)
と表される．重み W は，行列の各要素それぞれに与
えられるものである．しかし学習に用いるのは元々値
が観測されるインデックスのみである．
4.1
提案法：NDCG@k を重視した重み付け
本研究はコンテンツを推薦する際，ユーザーのニーズ
はランキング上位に集中していることを考慮している
ため，評価基準に Normalized Discounted Cumulative
Gain（NDCG）を用いることにする．NDCG@k は学
習によって得られた行列で作成するランキングを評価
するための手法であり，第 k 位までのドキュメントの
み評価に用いることでランキング上位の整合性を重視
するものである．ランキング学習では上位にランクさ
れるドキュメントが下位にランクされるドキュメント
に比べて重視されるため，それを考慮した重み付けが
望ましい．具体的には推薦システムでは検索結果を 5
平成 23 年度創成シミュレーション工学専攻修士論文梗概集
計算応用科学分野
件または 10 件表示するものが一般的であるため，評価
実験では特に k = 5 および 10 の NDCG 値に注目する
のがよいと考えられる．NDCG 値は 0 から 1 までの値
を取り，以下のように算出される．
{
k
1 ∑ 2yr − 1, r = 1，
NDCG@k =
(3)
2yr −1
r ≠ 1．
Zk
log (r) ,
i=1
validNDCG
1
nowait
flat
adapt
0.95
0.9
0.85
NDCG@k
0.8
0.75
0.7
2
0.65
r は順位であり，Z k の値は，評価用レビューまたは，
学習前のレビューの値をそのまま並び変えた正解の値
である．この評価基準を考慮した重み付けの算出方法が
{
W i,j = 1.0，
r = 1，
yr
(4)
W i,j = 2ymax2 log (r) ， r ≠ 1，
図 2: 検証データでの NDCG@k ，各アプローチの比較
であり，提案法では，以下のアルゴリズムによって行
列補完を行う．
'
$
5.1
2
0.6
0.55
0.5
1
2
3
4
5
6
7
8
9
10
k
考察
X, K を初期化する
学習データでも検証データでも，提案法である重み
X の各行に関して，
付き SVD の NDCG 値が大きいことから，この評価基
minxi W Ω (yi − xi )2Ω + γxTi K −1 xi
準では重みを付けたアプローチが有効であることがわ
を解き，xi を求める．
かった．
(3) 特異値分解を用い，K を更新する．
要素の順位を求め，W を更新する．
(4) 収束するまで (2)(3) を繰り返す．
&
% 6
まとめ
(1)
(2)
5
実験
実験の元となるデータセットは，6040 人のユーザー
の 3952 タイトルの映画に対するレビューが入力されて
いる MovieLens データであり，得られているレビュー
は全体の 10% 以下，レビューの段階は 5 段階である．
その行列を，訓練データと検証データに分割し，実験
を行う．正則化パラメータ γ = 1.0 とし，データセッ
ト 5 パターンで，重み無し LRA，固定重み LRA と提
案法である適応的重み付き LRA の訓練データと検証
データでの NDCG@k を比較する．
trainNDCG
1
nowait
flat
adapt
0.98
0.96
本研究では，ランキングを考慮した推薦システムの
ための行列補完問題に個別重みを導入した．評価基準
に NDCG@k を用いた重み付き Low-Rank Approximation を提案し，従来法よりもよい NDCG 値を出すこと
に成功した．今後の課題として，膨大なデータを扱う
ためアルゴリズムの高速化，汎化性能を向上させるた
めに交差検証などを用いて最適な正則化パラメータ選
択等が挙げられる．
参考文献
[1] Kai Yu，Shenghuo Zhu，John Lafferty and Yihong Gong．Fast Nonparametric Matrix Factorization for Large-scale Collaborative Filtering．SIGIR，2009．
[2] 土方嘉徳，情報推薦・情報フィルタリングのため
のユーザプロファイリング技術，人工知能学会誌，
Vol.19, No.3, pp.365–372, 2004.
NDCG@k
0.94
0.92
0.9
0.88
1
2
3
4
5
6
7
8
9
10
k
図 1: 訓練データでの NDCG@k ，各アプローチの比較
[3] Wang, Jun and de Vries, Arjen P. and Reinders,
Marcel J. T. Unifying user-based and item-based
collaborative filtering approaches by similarity fusion.SIGIR ,2006.
[4] Sarwar, Badrul and Karypis, George and Konstan,
Joseph and Reidl, John. Item-based collaborative
filtering recommendation algorithms.WWW,2001.