はじめに

　Pythonで使えるベイズ最適化ライブラリを用いてOpenFOAMのケーススタディを行っていく流れを説明していきます。

　OpenFOAMを用いることである一つの条件下での流体の挙動を計算することができます。ただ実際に用いる際には"複数の条件のうち最適なものはどれか"を探し出したいことが多いかと思います。総当りで試してみるのもよいですが、計算時間やリソースが限られている中ではできる限り計算パターン数を減らしたいものです。 OpenFOAMの最適条件の探索にはDAKOTAがよく用いられており分かりやすく紹介いただいているサイトも多くあります。

Howto DAKOTA OpenFOAM - OpenFOAMWiki

　ただ近年の機械学習ブームにより多くの最適化ツールなどがたくさん登場してきています。せっかくだからこれらをCAEに活用できないかなーと思い今回はDAKOTAではなく機械学習のハイパーパラメータの最適化などで用いられるベイズ最適化ツール（Python）を用いて最適値探索を行いました。以下の流れで最適化を行っていきます。

OpenFOAMのケースの作成
OpenFOAMで計算し値を返すPython関数を作成
その関数を繰り返し実行しベイズ最適化により探索

　本記事では全体の流れをかいつまんで書いていきます。それぞれの詳細は別記事で書く予定です。本当はこの流れで強化学習もやってみたかったのですが、かなり難しかったのでその辺はまた今度ということで。

　ちなみにこのベイズ最適化で遊んでいたのが冬季オリンピックの頃でしたので、その影響でパシュートの空気抵抗を題材にしています。以下のように二人の選手が縦に並んで滑っている時に、後ろの選手は【どのくらい離れて】【どのくらいタイミングをずらして】【どのくらいの幅で】滑ると最も空気抵抗が少なくなるのかを探索してみました。

f:id:inabower:20190404202706g:plain

はじめに
チームパシュート
- 免責事項
計算環境
1. OpenFOAMのケースの作成
2.OpenFOAMの実行を行うためのPython関数を作成
3. Pythonのライブラリを使ってベイズ最適化により探索
最後に

チームパシュート

　チームパシュートはスピードスケートの競技の一つであり、４人または３人で列を作りそのスピードを競う競技です。

チームパシュート - Wikipedia

　４人で縦に列を作るとスリップストリームが形成されて二人目以降の空気抵抗が軽減され、体力の温存などの効果が見込まれます。

スリップストリーム - Wikipedia

　スピードスケートの場合は自転車や自動車とは異なり左右に動きますのでこのスリップストリームを最大限に活かすことができるタイミングなどが存在するのではないかと考えました。

免責事項

計算結果に対して一切の責任を負いかねます
以下の理由から計算結果と実現象との間に許容できない誤差が存在する可能性が高いです
- 私はスピードスケートに関する知見を一切持っていません
- 限られた条件(正弦振動、速度・振動数固定)のみで検討しています
- 計算マシンの関係上、計算不可の少ないメッシュ及び乱流モデルを使用しています

計算環境

OS : Ubuntu 18.04
ソルバー：pimpleFoam(OpenFOAM v1812)
メッシャー：SALOME v8.4.0
Python 3.6.6
- PyFoam : 0.6.10
- paraview : 5.6.0 conda-forge

1. OpenFOAMのケースの作成

　まずは普通に二人の選手の空気抵抗を計算するケースを作成します。

　今回はPCのスペックの都合上以下のモデルで計算を行います。

ソルバー : pimpleFoam
移動メッシュ : motionSolver
乱流モデル : k-εモデル（RANS）
非圧縮ニュートン流体
二次元メッシュ（ヘキサメッシュ、23718セル）

　こんな大きさで

f:id:inabower:20190114142714p:plain

　こんなメッシュを使います。

f:id:inabower:20190210092820p:plain

　こんな感じで動きを制御すると

0/pointDisplacement

FoamFile
{
    version     2.0;
    format      ascii;
    class       pointVectorField;
    object      pointDisplacement;
}

dimensions      [0 1 0 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    ".*"
    {
        type            uniformFixedValue;
        uniformValue    (0 0 0);
    }
    human
    {
        type            oscillatingDisplacement;
        amplitude       (0.2 0 0);
        omega           6.28;
        value           uniform (0 0 0);
    }
    human2
    {
        type            oscillatingDisplacement;
        amplitude       (0.2 0 0);
        omega           6.28;
        value           uniform (0 0 0);
        
    }
}

constant/dynamicMeshDict

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dynamicFvMesh dynamicMotionSolverFvMesh;
motionSolverLibs ("libfvMotionSolvers.so");
solver  displacementLaplacian;

displacementLaplacianCoeffs
{
    diffusivity     inverseDistance ( human human2 );
}

　こんな感じで計算できました。

f:id:inabower:20190329200406g:plain

　今回は空気抵抗の指標として抗力係数Cdが最も小さい条件を探索します。

2.OpenFOAMの実行を行うためのPython関数を作成

　上のケースを使って二人目の動きをPythonのライブラリで最適化することを目指してみます。

　今回は以下のベイズ最適化ライブラリを使用します。

fmfn/BayesianOptimization

/* ここからREADMEの抜粋

　以下の関数がもしもブラックボックスであった時に

def black_box_function(x, y):
    return -x ** 2 - (y - 1) ** 2 + 1

　以下のように実行すると

from bayes_opt import BayesianOptimization

# Bounded region of parameter space
pbounds = {'x': (2, 4), 'y': (-3, 3)}

optimizer = BayesianOptimization(
    f=black_box_function,
    pbounds=pbounds,
    random_state=1,
)

optimizer.maximize(
    init_points=2,
    n_iter=3,
)

|   iter    |  target   |     x     |     y     |
-------------------------------------------------
|  1        | -7.135    |  2.834    |  1.322    |
|  2        | -7.78     |  2.0      | -1.186    |
|  3        | -19.0     |  4.0      |  3.0      |
|  4        | -16.3     |  2.378    | -2.413    |
|  5        | -4.441    |  2.105    | -0.005822 |
=================================================

　以下のように関数の最大値を見つけ出してくれます。

print(optimizer.max)

 {'target': -4.441293113411222, 'params': {'y': -0.005822117636089974, 'x': 2.104665051994087}}

ここまでREADMEの抜粋 */

　このBayesianOptimizationを使って空気抵抗の最小値を目指したいと思います。

　READMEから以下のような関数を作ることができればbayes_optが使えるということになります。

def func(orgCase, width, dt, dist):
    # 2-1. ベースとなるケースをコピー
    # 2-2. 引数の移動幅、時間差、距離を使ってケースの計算条件を変更
    # 2-3. ソルバーを実行
    # 2-4. Cdの平均値を算出
    return result

なお引数は以下の条件を設定しました。

width : 移動する横幅 [m]
dt : 前の選手とのタイミングのズレ [秒]
dist : 前の選手との距離 [m]

　関数の中身についてかいつまんで説明していきます。

2-1. ベースとなるケースをコピー

　まずshutil.copytreeなどを使いつつ名前が重複しないように好きな場所に自動でコピーしてくれるように作ります。

例) 重複しないようなコピー

import os, shutil

baseName = '../test/{}_'.format(orgCase.split('/')[-1])
n = 0
newName = baseName + str(n)

while os.path.exists(newName):
    n += 1 # 該当する名前のディレクトリが存在する限りループ
    newName = baseName + str(n)
shutil.copytree(orgCase, newName) # ループから抜け出したらコピー

2-2. 引数の移動幅、時間差、距離を使ってケースの計算条件を変更

　条件の変更にはPyFoamを用います。PyFoamではOpenFOAMのファイルをPythonのdictとして読み込んだり書き換えたりすることができます。

例) 0/pointDisplacementの一部を変更する

boundaryField
{
  ".*"
  {
    type noSlip;
  }
  human
  {
    omega      6.28;
    amplitude (0.2 0 0);
  }
}

というファイルがあったとして

from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile

with ParsedParameterFile(caseName + '/0/pointDisplacement') as pointDisplacement
    pointDisplacement.content['boundaryField']['human']['amplitude'] = [1.0, 0, 0]
    pointDisplacement.content['boundaryField']['human']['omega'] = 3.14
    pointDisplacement.writeFile()

と実行すると

boundaryField
{
  ".*"
  {
    type noSlip;
  }
  human
  {
    omega      3.14;
    amplitude (1.0 0 0);
  }
}

というように書き換えることができます。

2-3. ソルバーを実行

　ソルバーの実行はPyFoamやsubprocessなどを使います。

例) PyFoamのBasicRunnerでpimpleFoamを実行する

from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile
from PyFoam.Execution.BasicRunner import BasicRunner

with ParsedParameterFile(caseName+'/system/controlDict') as controlDict:
    app = controlDict.content['application']

runner = BasicRunner([app,'-case',caseName], silent=True)
runner.start()

　もしくはPythonの標準ライブラリのsubprocessを使ってAllrunなどのスクリプトをそのまま実行するという手もあります。

例) subprocessでAllrunを実行する

import os, subprocess
curDir = os.getcwd()  # 現在のディレクトリを取得
os.chdir("path/to/case") # ケースの場所に移動
run = subprocess.check_output(['./Allrun']) # Allrunスクリプトを実行
os.chdir(curDir) # 元の場所に戻って来る

2-4. Cdの平均値を算出

　最後にOpenFOAMのfunctionObjectで出力されるファイルをpandasで読み込み計算します。

import pandas as pd

fileName = caseName + '/postProcessing/forceCoeffs/0/coefficient.dat' # functionObjectで得られたファイル
df = pd.read_csv(fileName, header=9, sep='\t', index_col=0, names=['time', 'Cm', 'Cd', 'Cl', 'Cl(f)', 'Cl(r)'])
meanCd = df[df.index>1.0].Cd.mean() # 1秒以降のCdの平均値
return -meanCd # bayes_optでは最大値を求めるため出力の符号を逆転している

　というわけでPythonの関数ができました。

3. Pythonのライブラリを使ってベイズ最適化により探索

　それではようやく最適化を行います。公式ドキュメントを頑張って読みながら以下のように実行します。

width = 移動する横幅 [m] （前の選手は0.2m幅で振幅している）
dt = 前の選手とのタイミングのズレ [秒]
dist = 前の選手との距離 [m] - 2.0

bayes = BayesianOptimization(func, {
        'width': (0.19,0.21),
        'dt': (-0.25, 0.25),
        'dist': (0.0, 2.0)
        })
logger = JSONLogger(path="./logs.json")
bayes.subscribe(Events.OPTMIZATION_STEP, logger)

start_time = timer(None)
with warnings.catch_warnings():
    warnings.filterwarnings('ignore')
    bayes.maximize(init_points=2, n_iter=100, acq='ei', xi=0.0)
timer(start_time)

　すると以下のようなjsonファイルが吐き出されます。

{"target": -0.12488855661832597, "params": {"dist": 1.0123251860356828, "dt": 0.11468345268645169, "width": 0.20300313735322253}, "datetime": {"datetime": "2019-02-09 11:53:45", "elapsed": 0.0, "delta": 0.0}}
{"target": -0.11350154262227637, "params": {"dist": 0.23562374526490482, "dt": 0.011811599547604579, "width": 0.20843863840949256}, "datetime": {"datetime": "2019-02-09 12:15:00", "elapsed": 1275.108646, "delta": 1275.108646}}
{"target": -0.12666777807360347, "params": {"dist": 1.0903919796249917, "dt": 0.21091296197416276, "width": 0.2021055109659244}, "datetime": {"datetime": "2019-02-09 12:35:56", "elapsed": 2530.687802, "delta": 1255.579156}}
//中略
{"target": -0.12650695477299723, "params": {"dist": 1.0099823018034515, "dt": 0.24783059218889988, "width": 0.19}, "datetime": {"datetime": "2019-02-11 12:01:44", "elapsed": 173278.716431, "delta": 1290.549481}}
{"target": -0.12570835542030734, "params": {"dist": 1.2199388210282855, "dt": 0.15748704600468058, "width": 0.19}, "datetime": {"datetime": "2019-02-11 12:22:31", "elapsed": 174526.325504, "delta": 1247.609073}}
{"target": -0.12685274399778948, "params": {"dist": 1.1768319169307369, "dt": 0.22925668553124134, "width": 0.21}, "datetime": {"datetime": "2019-02-11 12:43:41", "elapsed": 175796.589145, "delta": 1270.263641}}
{"target": -0.12692054790359833, "params": {"dist": 0.8905474973894918, "dt": 0.2421151423565294, "width": 0.21}, "datetime": {"datetime": "2019-02-11 13:04:33", "elapsed": 177047.836131, "delta": 1251.246986}}
{"target": -0.12636646810634, "params": {"dist": 1.2624964305717332, "dt": 0.23739564249127137, "width": 0.19}, "datetime": {"datetime": "2019-02-11 13:25:59", "elapsed": 178333.773986, "delta": 1285.937855}}
{"target": -0.10762280026261561, "params": {"dist": 0.41815714564885614, "dt": -0.1763589800626097, "width": 0.21}, "datetime": {"datetime": "2019-02-11 13:46:26", "elapsed": 179560.734808, "delta": 1226.960822}}
{"target": -0.1270092530489997, "params": {"dist": 1.6665504479534357, "dt": 0.1686015459885953, "width": 0.21}, "datetime": {"datetime": "2019-02-11 14:07:45", "elapsed": 180839.915713, "delta": 1279.180905}}

"target"が関数から得られたスコアになります。この値が最大になるときが空気抵抗最小ということになります。

bayes.max

{'target': -0.09517227800080204,
 'params': {'dist': 0.0,
  'dt': -0.12149279231617968,
  'width': 0.2072435416122845}}

というわけで最適値を得ることができました。検証した中では以下の値が空気抵抗の最も小さい条件ということになります。

距離は2.0m（検討した中で最も近い値）
タイミングは0.12秒遅れ
幅は前の人より7cm外側

以下のような動きになります。

f:id:inabower:20190404202706g:plain

　ちなみに100回計算の各パラメータを座標として、スコアを球の大きさと色で、ParaViewで描画してみました。

ベイズ最適化（3変数）の各スコアをParaViewで可視化してみた pic.twitter.com/rXvqzscxLc
— inabower (@public_inabower) 2019年2月26日

　こうやって見ると距離は近ければ近いほどよいことが分かりますしタイミングや幅も最適な場所がありそうに見えます。

最後に

　相当駆け足でしたがOpenFOAMでベイズ最適化を行う流れを説明してみました。同じ感じでPythonのツールであるOptunaとかKerasとかも使用できると思います。もしよかったら試してみてください。

ただし空気抵抗は無視しないものとする

流体解析の勉強の題材として日頃気になったことをまとめたり確かめたりしていく予定です。OSSしか使わないつもりです。

【OpenFOAM×ベイズ最適化】パシュート（スケート）の二人目の空気抵抗が最も小さい条件を探索する