OpenMDAO を使ってみる
2014年6月8日
はじめに
OpenMDAO を使ってみる。
使用バージョン
OpenMDAO 0.9.7
OpenMDAO とは？
OpenMDAO (Opensource Multidisciplinary Design Analysis and Optimization) は、さまざまなコンポーネントから構成されたシステムの設計や、解析、最適化を行うためのプラットフォームである。
問題
関数 f(x) = x**2 - 2 を作って、2 の平方根を求めてみる。
プロジェクトの作成
OpenMDAO を有効化、GUI を起動。
$ . ~/opt/openmdao-0.9.7/bin/activate $ openmdao gui
"New Project" ボタンで、適当な名前 (たとえば "sqrt") でプロジェクトを作成する。
関数の作成
左側の "Files" タブのメニュー[New File] で "square.py" というファイルを作成する。
"square.py" をダブルクリックして編集モードにし、以下の内容にする。
from openmdao.main.api import Component
from openmdao.lib.datatypes.api import Float
class Square(Component):
x = Float(0.0, iotype='in', desc='x')
f = Float(0.0, iotype='out', desc='f')
def execute(self):
x = self.x
self.f = x**2
ここで、class 名 "Square" はコンポーネント名に、"x"、"f" はそれぞれコンポーネントの入力と出力になる。ここでは "f = x**2" としている。必要な関数は f(x) = x**2 - 2 だが、操作方法を確認するため、わざとめんどうな作り方をしている。このあと "- 2" 分のコンポーネントを作り、2 つのコンポーネントをつないで f(x) = x**2 - 2 とする。
"square.py" の編集が済むと、右側の "Library" タブに "Square" という項目が追加される。"
"Square" を "Dataflow" タブの領域にドラック & ドラッグする。名前を聞かれるので、適当に "s" にでもしておく。
"Square" のアイコンを右クリックするとポップアップメニューが出るので、そこから "Properties" を選ぶと、プロパティが表示される (右側の "Properties" タブでもよい)。そこに入力 "x" や 出力 "f" が表示されている。
コンポーネントの動作を確認するため、"x" に 2 を入れ、"Square" のアイコンを右クリックし、ポップアップメニューから "Run" を選ぶ。コンポーネントが実行され、"f" の値が 2 の 2 乗である 4 になることが確認できる。
コンポーネントの動作が確認できたので、本番の関数作りをはじめる。新しいファイル "minus.py" を作り、以下のような内容にする。
from openmdao.main.api import Component
from openmdao.lib.datatypes.api import Float
class Minus(Component):
x = Float(0.0, iotype='in', desc='x')
a = Float(0.0, iotype='in', desc='a')
f = Float(0.0, iotype='out', desc='f')
def execute(self):
x = self.x
a = self.a
self.f = x - a
これは入力 "x"、"a"、出力 "f" をもつコンポーネント "Minus" で、"f = x - a" を計算する。"Square" とこの "Minus" をつなげて f(x) = x**2 - 2 を構成する。
このままだとコンポーネントをつなげられないので、一旦 "Square" のアイコンを右クリックし、ポップアップメニュー "Remove" で削除する。そして右の "Library" から "Assembly" を "Dataflow" にドラッグ & ドロップする。名前は "a" とでもしておく。
続いて、"Square"、"Minus" を "Assembly" の内部にドラッグ & ドロップする。名前はそれぞれ "s"、"m" とでもしておく。
"Square" と "Minus" をつなぐ。"Square" の周りにある緑の点のどれかから "Minus" の周りの点のどれか (大きな丸が表示されるもの)" へとドラッグする。
出力と入力をつなぐためのダイアログが表示されるので、"s" 側の "f" から "m" 側の "x" までドラッグし、線をつなぎ、ダイアログを閉じる。
"Square" のアイコンと "Minus" のアイコンにコネクタが追加される。
動作確認のため、"Square" の "Properties" で "x" に 2 を設定して "Square" を実行、続いて "Minus" の "a" に 2 を設定して "Minus" を実行する。"Minus" の実行時に "x" が "Square" の "f" の値 4 になり、"f" が 2**2 - 2 = 2 となることが確認できる。
この状態では "Square" と "Minus" はバラバラにしか実行されないので、"ドライバ" に入れる。真ん中の "Workflow" タブを表示し、左の "Objects" タブのリストから "Components" を選び、"s" と "m" を "Workflow" のドライバ "a.driver" の右下の箱の中にドラッグ & ドロップする。
この状態で "Dataflow" にて "driver" ("Run_Once")を実行すると、まず "Square" を計算し、その結果を用いて "Minus" を実行する、という一連の流れ (ワークフロー) が実行される。"Square" の "x" を変えて "driver" を実行すると "Minus" の "f" の値が変わることが確認できる。
これで必要な関数が完成した。
実験計画法
ドライバは、ワークフローに対して何かをするものである。ここでは実験計画法 (DOE) ドライバを使ってみる。
右の "Library" から "DOEdriver" を、"driver" ("Run_Once") にドラッグ & ドロップする。ドライバを入れ替えてよいか聞かれるので、"Ok" を押す。
ドライバに入力パラメタ (ここでは "s" の "x") を設定する。"driver" のアイコンをダブルクリックする。"Parameters" タブの "Add Parameter" を押す。"Target" に "s.x" を入力する。"Low" は 0、"Heigh" は 10 にでもしておく。他は空欄のまま "Ok" を押す。
"driver" から "s" にコネクタができる。
再度 "driver" のアイコンをダブルクリックし、"Slots" タブを選ぶ。左の "Library" の一番上の "Search" で "DOE" と打ち込むと、"DOE" 関係の項目が表示される。この中の "FullFactorial" を "Slots" の "DOEgenerator IDOEgenerator" にドラッグ & ドロップする。
"num_leves" を聞かれるので 11 でも入れておく (あとで変更可能)。この数字は、入力を何分割して入れるかを指定する。たとえば入力値が 0〜10 で、分割数が 3 だったら、入力に 0、0.5、10 を入れてそれぞれ計算される。
今度は "Library" の "Search" に "recorder" と打ち込む。その中から "DumpCaseRecorder" と "CSVCaseRecorder" を、"Slots" の "recorders ICaseRecorder" にドラッグ & ドロップする。何か聞かれるが、どちらも "Ok" でよい。
これらは、出力の設定である。"DumpCaseRecorder" は、計算結果をを下の出力ウインドウに出力する。"CSVCaseRecorder" は計算結果を CSV 形式で出力する。CSV の設定は "driver" の "Inputs" タブの "case_outputs" で行う。この値を ["m.f"] としておけば、入力 ("s.f") に対する "m.f" の値のリストが得られる。
"driver" を "Run" する。下の出力ウインドウに入力値 "s.x" と出力値 "m.f" の結果が出力される。左のタブ "Files" には "cases.csv" というファイルが追加されている。
"timestamp","label","/INPUTS","s.x","/OUTPUTS","m.f","/METADATA","retries","max_retries","parent_uuid","msg" 1402238247.420257,"","",0.0,"",-2.0,"",0,1,"","" 1402238247.42517,"","",1.0,"",-1.0,"",0,1,"","" 1402238247.429745,"","",2.0,"",2.0,"",0,1,"","" 1402238247.434337,"","",3.0000000000000004,"",7.0000000000000036,"",0,1,"","" 1402238247.438705,"","",4.0,"",14.0,"",0,1,"","" 1402238247.440307,"","",5.0,"",23.0,"",0,1,"","" 1402238247.441839,"","",6.0000000000000009,"",34.000000000000014,"",0,1,"","" 1402238247.443291,"","",7.0000000000000009,"",47.000000000000014,"",0,1,"","" 1402238247.444886,"","",8.0,"",62.0,"",0,1,"","" 1402238247.4464,"","",9.0,"",79.0,"",0,1,"","" 1402238247.447885,"","",10.0,"",98.0,"",0,1,"",""
f(x) = x**2 - 2 で f = 0 となる x を求めたいので、この結果から x = 1〜2 ということがわかる (答えはわかっているが、わかっていない体で)。
最小化
ワークフローの出力を最小化するドライバもある。そのひとつは "SLSQPdriver" である ("Library" の "Search" で "optimizer" と打ち込むと出てくる)。
"driver" に "SLSQPdriver" をドラッグ & ドラッグして、ドライバを入れ替える。設定項目は、"Parameters"、"Objectives"、必要に応じて "Constraints" である。
"Parameters" の設定は "ODEdriver" と同じである。ドライバを入れ替えたなら、"ODEdriver" の設定がそのまま残っている。
"Objectives" の "Add Objective" で目的のパラメタ (目的関数) を設定する。ここでは "Expression" に "m.f" と設定する。
このドライバは目的関数を最小化する。ここでは目的関数を 0 にしたいが、このままだと目的関数が負に行ってしまうので、制約条件を加える。"Constraints" の "Add Constraint" で、"Expression" に "m.f >= 0" を設定する。
"s" の "Properties" で "x" に初期値として 1 を設定する。そして "driver" を "Run" する。そうすると、"s" の値として 1.41421356...、つまり 2 の平方根が得られる。
Python による実行
GUI を使わずに Python だけで上記の操作を行うこともできる。たとえば、実験計画法の手順と同じことを行うスクリプトは以下のようになる。
assembly.py
from openmdao.main.api import Assembly
from openmdao.lib.drivers.api import DOEdriver
from openmdao.lib.doegenerators.api import FullFactorial
from openmdao.lib.casehandlers.api import ListCaseRecorder
from square import *
from minus import *
class Analysis(Assembly):
def configure(self):
self.add('s', Square())
self.add('m', Minus())
self.add('driver', DOEdriver())
self.driver.DOEgenerator = FullFactorial(11)
self.driver.add_parameter('s.x',low=0, high=10)
self.driver.case_outputs = ['m.f',]
self.driver.recorders = [ListCaseRecorder(),]
self.driver.workflow.add('s')
self.driver.workflow.add('m')
self.connect('s.f', 'm.x')
if __name__ == "__main__":
import time
analysis = Analysis()
analysis.m.a = 2.
tt = time.time()
analysis.run()
print "Elapsed time: ", time.time() - tt, "seconds"
from pylab import *
x = []
y = []
for c in analysis.driver.recorders[0].get_iterator():
print "%f %f" % (c['s.x'], c['m.f'])
x.append(c['s.x'])
y.append(c['m.f'])
plot(x, y, 'r-o')
xlabel('$x$')
ylabel('$f$')
title('$x^2 - %d$' % analysis.m.a)
show()
square.py、minus.py を置いておく必要がある。
OpenMDAO の環境を有効にして、実行。
$ python assembly.py Elapsed time: 0.0429668426514 seconds 0.000000 -2.000000 1.000000 -1.000000 2.000000 2.000000 3.000000 7.000000 4.000000 14.000000 5.000000 23.000000 6.000000 34.000000 7.000000 47.000000 8.000000 62.000000 9.000000 79.000000 10.000000 98.000000
つぎのようなグラフが表示される。
