FMUのパラメータ推定ライブラリmodestpyを使う

概要

表題のとおりmodestpyはFMU(Functional Mockup Unit)のパラメータ推定をするpythonのライブラリで、ライセンスはBSD 2-clauseライセンスである。FMUの計算自体はJModelicaのpyFMIを使っている。ある時系列データに対してシミュレーション結果とのRMSEが小さくなるようなパラメータを探索するようになっている。資料はAPIのドキュメントとカンファレンスのプレプリントのPDFがある。

使用バージョン
- modestpy 0.0.7
- JModelica2.4 →Modelica標準ライブラリ3.2.2

インストール

condaからインストールするのが推奨とサイトには書いてあるけれど、JModelicaのpython上で使えたほうがよさそうなのでコマンドプロンプトで以下のようにしてgithubの最新版のmedestpyをインストールした。

C:\JModelica.org-2.4\Python27\Python_64\python.exe -m pip install https://github.com/sdu-cfei/modest-py/archive/master.zip

そうするとC:\JModelica.org-2.4\Python27\Python_64\Lib\site-packages\modestpyにインストールされる。インストールしたらJModelicaのIPythonを起動してテストファイルを実行してみる。とりあえず何かしら計算をしているようだがよくわからない。

from modestpy.test import run
run.tests()

サンプルの実行

exampleフォルダがgithubにあるのでとりあえず作業フォルダにコピーしてくる。example\simpleの中にスクリプトがいくつかあり、resourcesフォルダの中にいくつか設定用のファイルがある。とりあえず計算に使うモデルSimple2R1C_ic.moをOpenModelicaで開いてみる。フォルダの中にはいくつかのプラットフォーム用のfmuファイルがあるものの、すべての元はこのファイルである。以下のようなモデルになっている。

熱容量1000[J/W]の点に熱抵抗１[K/W]が二つつながっており、抵抗の反対側は片方がTi1[K]、もう片方がTi2[K]になっている。熱容量の点は初期温度Tstart=20で、温度Tはセンサーで取得している。
このモデルのパラメータを推定するらしい。モデルのなかでparameterとして定義されているのは二つの熱抵抗R1とR2、初期温度Tstart、熱容量Cの４つである。温度の時系列データを入力して、Tstart以外の3つのパラメータを推定するのだが、まずは広く探索するソルバーで計算して途中から局所的に探索していくらしい。よくわからないのでとりあえず実行してみる。
まずはexample\simple\simple.pyを作業フォルダにコピーしてくる。 simple.pyを開いてif __name__ == "__main__":をコメントアウトし、そこから先を一括選択してShift+tabでインデントを元に戻す。（操作はエディタにもよるかもしれない）
あとはJModelicaのIPythonで以下のようにすれば実行される。(とくに編集せずにコマンドプロンプトからpython simple.pyでもよいと思う。)

cd '作業フォルダ'
import simple

なにやらエラーのような終わり方をしている。ただし計算はできていそうだった。

計算結果は作業フォルダ/examples/simple/workdirに保存されるようになっている。
summary_1.csv(左)とsummary_2.csv(右)に以下のようにパラメタスタディのログが出力されている。GAが広域探索として遺伝的アルゴリズムを、PSが局所探索としてパターンサーチを使用している部分である。局所探索は他にもScipyにある探索が使えるらしい。

best_per_run.csvにそれぞれのなかでもっともよいパラメータを書いている。

そしてfinal.csvにもっともよい結果が出力される。

それからいくつかの画像が出力される。
GAの結果（計算回数の枚数だけ出力）

PSの結果（計算回数の枚数だけ出力）

RMSEの推移

推定したパラメータでの計算結果と推定に用いた実測値の比較

とりあえずうまくパラメータを推定できたようである。

設定の確認

simple.pyは以下のようになっている。

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import json
import os
import pandas as pd
from modestpy import Estimation
from modestpy.utilities.sysarch import get_sys_arch

#if __name__ == "__main__":
"""
This file is supposed to be run from the root directory.
Otherwise the paths have to be corrected.
"""
# DATA PREPARATION ==============================================
# Resources
platform = get_sys_arch()
assert platform, 'Unsupported platform type!'
fmu_file = 'Simple2R1C_ic_' + platform + '.fmu'
fmu_path = os.path.join('examples', 'simple', 'resources', fmu_file)
inp_path = os.path.join('examples', 'simple', 'resources', 'inputs.csv')
ideal_path = os.path.join('examples', 'simple', 'resources', 'result.csv')
est_path = os.path.join('examples', 'simple', 'resources', 'est.json')
known_path = os.path.join('examples', 'simple', 'resources', 'known.json')

# Working directory
workdir = os.path.join('examples', 'simple', 'workdir')
if not os.path.exists(workdir):
    os.mkdir(workdir)
    assert os.path.exists(workdir), "Work directory does not exist"
# Load inputs
inp = pd.read_csv(inp_path).set_index('time')
# Load measurements (ideal results)
ideal = pd.read_csv(ideal_path).set_index('time')
# Load definition of estimated parameters (name, initial value, bounds)
with open(est_path) as f:
    est = json.load(f)
# Load definition of known parameters (name, value)
with open(known_path) as f:
    known = json.load(f)

# MODEL IDENTIFICATION ==========================================
session = Estimation(workdir, fmu_path, inp, known, est, ideal,
                        lp_n=2, lp_len=50000, lp_frame=(0, 50000),
                        vp=(0, 50000), ic_param={'Tstart': 'T'},
                        methods=('GA', 'PS'),
                        ga_opts={'maxiter': 5, 'tol': 0.001, 'lhs': True},
                        ps_opts={'maxiter': 500, 'tol': 1e-6},
                        scipy_opts={},
                        ftype='RMSE', seed=1,
                        default_log=True, logfile='simple.log')
estimates = session.estimate()
err, res = session.validate()

はじめの方は読み込みファイルのパスを変数に入れたりしており、MODEL IDENTIFICATION 以降が計算の実行部となっている。
session = Estimation(~)で変数を一括で設定し、
estimates = session.estimate()でパラメータ推定、
err, res = session.validate()で推定した値を使用して実験値と計算値の比較をしている。

設定に使用するファイルは以下の通り。

• fmu_file：FMUファイル、JModelicaなどで作成できる。
• inputs.csv：Ti1とTi2の時系列データ
• result.csv：inputs.csvの条件の時のTの実測データ
• est.json：パラメータ探索の設定
• known.json：Tの初期値Tstartの設定

パラメータ探索の設定est.jsonは以下のようになっている。だいたいの目安の値、最小値、最大値の順番に記載している。 "R1": [0.08, 0.001, 0.5],
"R2": [0.08, 0.001, 0.5],
"C": [1000.0, 100.0, 10000.0]

その他は以下のように設定している。

• lp_n=2：探索の回数
• lp_len=50000：１回の探索の長さ
• lp_frame=(0, 50000)：探索回数の最小値と最大値
• vp=(0, 50000)：探索後の検証の計算の長さ
• ic_param={'Tstart': 'T'}：初期値の対応関係
• methods=('GA', 'PS')：使用するソルバーの指定、現状GAは固定
• ga_opts={'maxiter': 5, 'tol': 0.001, 'lhs': True}：GAのソルバーの設定、詳細はmodestpy/estim/ga/ga.pyのコメント参照
• ps_opts={'maxiter': 500, 'tol': 1e-6}：PSのソルバーの設定、詳細はmodestpy/estim/ps/ps.pyのコメント参照
• scipy_opts={}：Scipyのソルバーの設定、詳細はmodestpy/estim/scipy/scipy.pyのコメント参照
• ftype='RMSE'：目的関数現状はNRMSEかRMSEのみ。
• seed=1：乱数使用時の種、指定しなければ現在時間などOS指定の値を使う。
• default_log=True, logfile='simple.log'：ログの設定

終わりに

今回は第6回Modelicaライブラリ勉強会でパラメタ最適化している話を聞いて、そういえば何かあったよなぁと思って調べてみた。まだ発展途上中のシンプルなライブラリだがしっかり計算はできるようなのでなかなかよさそうなライブラリだと思った。