2024-12-06

OpenModelicaのPythonからの変数の変更とシミュレーション実行

Modelica Python

Modelicaでモデルを開発していると、パラメトリックスタディなど、モデルを外部から操作・データ取得したいと思うことがあります。

OpenModelicaでもパラメータを変更してシミュレーションを行う仕組みはある*1のですが、可能ならPythonからパラメータを指定して、出力データも取得・描画したいです。

以前、FMPyというPythonライブラリを使って、OpenModelicaで作成したモデルをFMUに変換し、それをPythonから実行するということを行いました。

ohtayo.hatenablog.com

しかし、モデルによってはFMUへの変換ができないことがあるため、OpenModelica自体を使って実行できるとより便利です。

OpenModelicaには、コマンドで外部からシミュレーション実行などが可能なOpenModelica Shellと、それをPythonから操作できるOMPythonというライブラリが付属しています。

本記事ではOMPythonを使って、PythonからOpenModelicaで作成したモデルのパラメータを変更して、シミュレーション結果を取得する方法を記載します。

1. 開発環境構築

1.1 バージョン

Windows 11 Pro 23H2
OpenModelica v1.24.2 (64bit)
Python 3.9.13 64bit
Modelica Standard Library 4.0.0

1.2インストール

OMPythonのインストール方法は、READMEに記載があります。

github.com

Windowsでは、OpenModelicaをインストールする際にOMPythonもついてくるので、コマンドプロンプトを使って、以下のコマンドでインストールが可能です。

cd %OPENMODELICAHOME%\share\omc\scripts\PythonInterface
python -m pip install -U .

しかし、私の環境ではREADME.mdが無いという旨のエラーが出て、インストールできませんでした。上記のREADME.mdをダウンロードしてきて、PythonInterfaceフォルダに置くとインストールすることができました。

しかし、これをやるくらいなら、素直にレポジトリをCloneしてきてインストールしたほうが早そうです。

git clone https://github.com/OpenModelica/OMPython.git
cd OMPython
python -m pip install -U .

2. OMPythonの使い方

OMPythonの使い方の概要は以下にあります。これに則って、BouncingBallモデルを実行してみます。

openmodelica.org

2.1 動作確認

まずは、OMPythonがうまく動いているか確認するために、以下のコマンドをPythonで実行してみます。

from OMPython import OMCSessionZMQ
omc = OMCSessionZMQ()
print(omc.sendExpression("getVersion()"))

OpenModelica v1.24.2 (64-bit)という表記が出ていればOKです。

2.2 モデルの実行

最初に、モデルを実行すると実行ファイル等が多く出力されるので、モデル実行用のフォルダを作ります。

mkdir BouncingBall

次に、以下のコマンドで、PythonからBouncingBallのモデルを実行します。

omc.sendExpression("cd(\"./BouncingBall\")")
omc.sendExpression("loadFile(getInstallationDirectoryPath() + \"/share/doc/omc/testmodels/BouncingBall.mo\")")
omc.sendExpression("simulate(BouncingBall, stopTime=3.0, outputFormat=\"csv\")")

コマンドの内容としては、作成したBouncingBallフォルダにカレントディレクトリを移動し、モデルファイルを読み込んで、シミュレーションする、というものです。正しく実行できていれば、messages: LOG_SUCCESSを含む結果表示が出て、BouncingBall_res.csvがフォルダに出力されます。

ちなみに、loadFile()の実行後は、checkModel(BouncingBall)でモデルのチェックなど、OpenModelica上で可能な操作をスクリプトからも行うことができます。このあたりのスクリプトは以下に説明があります。

build.openmodelica.org

csvが出力されたら、それを読み込んでPythonで描画します。

import pandas as pd
from matplotlib import pyplot as plt
df = pd.read_csv("./BouncingBall/BouncingBall_res.csv", index_col=0)
plt.plot(df.index, df["h"])

実行結果が出力されました。

2.3 パラメータ変更後実行

モデルのパラメータを変更して実行する方法は複数あります。実用上は、2.3.1のXMLを書き換える方法か、2.3.3のモデルの値を書き換える方法を用途に応じて使い分ける形になるかと思います。

2.3.1 初期値XMLを書き換えて実行する

モデルを一度シミュレーションすると、BouncingBall.exeという実行ファイルが作成され、これが初期値をもとにシミュレーション計算を行います。初期値はBouncingBall_init.xmlに記載されているので、これを書き換えてexeを実行すれば、パラメータ変更時のシミュレーションができます。

Modelica Standard Libraryでは、Modelica.Blocks.Sources.Constantの値は初期値扱いになるので、変更したいパラメータをこのクラスで定義しておくと、xmlによる変更がしやすくなります*2。

setInitXmlStartValue()がxmlを書き換える関数なので、以下のようなスクリプトで初期値を書き換えて実行できるはずです*3。

omc.sendExpression("setInitXmlStartValue(\"BouncingBall_init.xml\", \"h\", \"3\", \"BouncingBall_init.xml\")")
omc.sendExpression("system(\"BouncingBall.exe\")")

しかし、私の環境ではどうしてもsetInitXmlStartValue()の結果がFalseになってしまい、xmlの書き換えができませんでした。そこで、Pythonでxmlを書き換えるスクリプトを用意しました。

def set_init_xml_start_value(xml_file, variable_name, new_start_value):
    # OpenModelicaのsimulate()で生成されるinit.xmlのstart値を書き換える
    import xml.etree.ElementTree as ET
    tree = ET.parse(xml_file)
    root = tree.getroot()
    # 要素を全て見つける
    scalar_variables = root.findall(".//ScalarVariable")
    # 指定された名前の要素を全て見つける
    for scalar_variable in scalar_variables:
        if scalar_variable.get("name") == variable_name:
            # 要素内のReal要素を見つける
            real_element = scalar_variable.find("Real")
            if real_element is not None:
                # start属性の値を更新
                real_element.set("start", new_start_value)
                break
    tree.write(xml_file)

xml_file = "./BouncingBall/BouncingBall_init.xml"
set_init_xml_start_value(xml_file, "h", "3")
omc.sendExpression("system(\"BouncingBall.exe\")")

先ほどと同様にCSVをグラフ化すると、結果が変更されていました。

2.3.2 exeの引数としてパラメータを与える方法

BouncingBall.exeという実行ファイルには引数を与えることができ、そこで初期値をxmlで指定した値から上書きすることもできます。

omc.sendExpression("system(\"BouncingBall.exe -override h=4\")")

しかし、こちらは1つしか引数が指定できないため、複数パラメータを変更したい場合などには不向きです。

2.3.3 モデルの値を書き換えて再実行

setComponentModifierValue()という関数を使うと、モデルの値を直接書き換えることができます。その後、最初と同様にsimulate()すれば、異なる値でシミュレーションができます。

omc.sendExpression("setComponentModifierValue(BouncingBall, h.start, $Code(=5))")
omc.sendExpression("simulate(BouncingBall, stopTime=3.0, outputFormat=\"csv\")")

しかし、これだとコンパイル、実行ファイルの作成もやり直すため、時間がかかります。今回のBouncingBallのモデルは軽量なのでそこまで問題ではないですが、モデル規模が大きくなったり、かなりの回数パラメータを変更してシミュレーションしなおすとなると、問題になってきますので、XML書き換えで対応できるならそのほうがいいと思います。

2.3.4 OMPythonの他のAPI(ModelicaSystem)を使った変更

これまで紹介したOMPythonの使い方は、OpenModelica ShellのコマンドをPythonから実行するという、直接的ですがややわかりにくいものでした。

OMPythonはよりモデルを柔軟・簡易に扱うAPIとして、ModelicaSystemを用意しています。これを使うと、以下のようにPython側から見てよりわかりやすい記述で、パラメータ変更や実行ができます。

from OMPython import OMCSessionZMQ
omc = OMCSessionZMQ()
model_path=omc.sendExpression("getInstallationDirectoryPath()") + "/share/doc/omc/testmodels/"
from OMPython import ModelicaSystem
mod=ModelicaSystem(model_path + "BouncingBall.mo","BouncingBall")
mod.buildModel()
mod.setParameters("radius=14")
mod.simulate(resultfile="tmpbouncingBall.mat")

しかし、私の環境だとこのAPIを使うとsimulateで失敗してしまい、うまく実行できませんでした。確実性を考えると、2.3.1のようなやり方がよさそうに思います。

2.4 結果表示

先ほどはスルーしましたが、simulate()実行時に結果出力フォーマットを指定できます。標準はMATLABで使われるmat形式ですが、上記ではCSVを指定して、pandas.read_csvで読み込みしています。読み込んだデータは以下のようになり、pythonで扱いやすい形式となっているので、これをデータ処理するなりグラフ化するなり自由にできます。

time	h	v	der(h)	der(v)	foo	flying
0	3	0	0	-9.81	2	1
0.006	2.9998234	-0.05886	-0.05886	-9.81	2	1
0.012	2.999294	-0.11772	-0.11772	-9.81	2	1
0.018	2.9984107	-0.17658	-0.17658	-9.81	2	1
0.024	2.997175	-0.23544	-0.23544	-9.81	2	1

3. まとめ

OMPythonを使って、Pythonからパラメータの変更、OpenModelicaによるシミュレーション実行、結果取得を行う方法を紹介しました。

今後、これを使って、Pythonの豊富なライブラリによるパラメトリックスタディやModelicaモデルと連携した最適化などをやっていきたいと思います。

*1:Parameter Sensitivities with OpenModelica — OpenModelica User's Guide v1.25.0-dev-220-g494b2520bc documentation

*2:もちろん、他のクラスでも初期値があるものはそれを変更できます

*3:Newbie: initializing state variables from OMShell :: OpenModelica

2024-02-03

遺伝的アルゴリズムを使ったポートフォリオの多目的最適化

Python 多目的最適化遺伝的アルゴリズムポートフォリオ最適化

先日、新NISA関連の調べものをしていた時、以下のような記事を見つけました。

chandraroom.com

リスク-リターン平面上に、数々の資産の組み合わせ＝ポートフォリオの運用結果をプロットし、あるリスクに対して最大のリターンが得られるポートフォリオの集合＝効率的フロンティアという概念があるそうです。

これは、多目的最適化でいうパレートフロントと同じ考え方です。多目的最適化では、複数のトレードオフの関係にある目的関数に対して、これ以上どちらの目的も同時に良くすることはできない解の集合＝パレートフロントを見つけることが目標です。

パレートフロントを厳密に求めるには時間がかかるため、遺伝的アルゴリズムなどを用いて近似的なパレートフロントを探索する方法が取られています。効率的フロンティア上のポートフォリオについても同様の方法を適用し、リスクとリターンを計算する目的関数を用意し探索すれば、良好なポートフォリオ候補を自動で作成できそうだと思ったので、トライしてみました。

ポートフォリオの成績は以下の条件で計算することにします。

資産の種類は現金、日本株、米国株、先進国株、新興国株、米国債券、先進国債券、日本REIT、米国REIT、先進国REIT、コモディティ、金の12個
ドル建資産は日毎為替レートで円換算して計算
ポートフォリオは毎月指定比率となるようにリバランス
配当金・分配金は毎月再投資
売買にかかる手数料や、為替手数料、税金等のコストは考慮しない

計算にはPythonを用います。Pythonの多目的最適化ライブラリpymooを練習として使ってみます。実は、pymooにはケーススタディとしてポートフォリオ最適化の事例があるのですが、リバランスなどは考慮されていないようでしたので、参考程度にとどめ、新規にスクリプトを作成します。

作成したPython スクリプトはGitHubに公開しています。

github.com

本記事の目次です。

1. データ取得
2. ポートフォリオの成績の計算
3. 多目的最適化
4. 結果
5. まとめ

1. データ取得

データはYahoo!ファイナンスから取得します。各資産クラスのデータは、ETF(上場投資信託)から以下の銘柄を選定して利用することにしました。

資産	対象ETF	ティッカー	対象指数	設定
現金	-	-	-	-
日本株	NEXT FUNDS TOPIX連動型上場投資信託	1306.T	TOPIX	2008
米国株	SPDR S&P 500 ETF Trust	SPY	S＆P500指数	1993
先進国株	バンガード・FTSE先進国市場（除く米国）ETF	VEA	FTSE先進国オールキャップ（除く米国）インデックス	2007
新興国株	バンガード・FTSE・エマージング・マーケッツETF	VWO	FTSEエマージング・マーケッツ・オールキャップ（含む中国A株）	2005
米国債券	バンガードトータル債券市場ETF	BND	ブルームバーグ・バークレイズ米国総合浮動調整インデックス	2007
先進国債券	SPDR ブルームバーグ･バークレイズ世界国債（除く米国）ETF	BWX	ブルームバーグ・バークレイズ・グローバル国債（除く米国）キャップド指数	2007
日本REIT	NEXT FUNDS 東証REIT指数連動型上場投信	1343.T	東証REIT指数	2008
米国REIT	バンガードREIT ETF	VNQ	MSCI US Investable Market Real Estate 25/50 Index	2004
先進国REIT	iシェアーズ先進国（除く米国） REIT ETF	IFGL	FTSE EPRA/NAREIT 先進国不動産除く米国	2007
コモディティ	iシェアーズ S&P GSCI コモディティ・インデックス・トラスト	GSG	S&P GSCI 商品指数	2006
金	SPDR ゴールド・シェア	GLD	ロンドン金価格	2004

選定基準としては各資産に該当するETFの中でなるべく長い期間のデータがあることです。日本国債券や新興国 REITといった資産も調べたのですが、適当なETFが見当たらなかったので、今回の対象からは除外しました。*1

Yahoo!ファイナンスからはYahoo!ファイナンス APIを使ってその銘柄に関するデータをダウンロードすることができます。今回選択したETFには分配金が出るものがありますので、株価と分配金のデータをダウンロードすることにします。以下のような形でURLを作成しアクセスすると、株価の時系列データと分配金実績をJSON形式で取得できます。

import requests
from io import StringIO
import json

# urlを作成しjsonを取得
ticker = "SPY" # ここでティッカーシンボルを指定
period = "11y" # 10年データが欲しいので長めにとって切り出す
interval = "1d"
prefix = "https://query1.finance.yahoo.com/v8/finance/chart/"
url = prefix + ticker + "?range=" +period + "&interval=" + interval + "&events=div&includeAdjustedClose=true"
response = requests.get(url)
str = StringIO(response.text)
j = json.load(str)

株価のデータはindicatorsのquote、分配金のデータはeventsのdividendsというフィールドに配置されていますので、それぞれのデータをpandas.DataFrameに変換します。

import pandas as pd
# 株価の取得
quote = pd.DataFrame()
timestamp = j['chart']['result'][0]['timestamp']
quote.index = pd.to_datetime(timestamp, unit="s")
quote.index.name = "timestamp"
quote['close'] = j['chart']['result'][0]['indicators']['quote'][0]['close']

# 分配金情報の取得
div_json = j["chart"]["result"][0]["events"]["dividends"]
div = pd.DataFrame(list(div_json.values()), columns=['amount', 'date'])
div.index = pd.to_datetime(div["date"], unit="s")
div.index.name = "timestamp"
div = div.drop("date", axis=1)
div.columns = ["dividend"]

取得したデータは15～20年ほどのデータ期間がありますが、今回は10年間(2013年末から2023年末)のポートフォリオの成績を見て最適化をすることにします。

それから、ドル建て資産は日本円に換算します。ドル円の為替情報は株価と同様にティッカーをUSDJPY=XとすればYahoo!ファイナンスから取得できます。為替データと株価データをpandas.merge_asofで結合→乗算して、円換算します。

# ここで株価データ取得
price = quote.copy()
# ここで為替データ取得
usdjpy = quote.copy()
usdjpy = usdjpy.dropna()

# 14時間足して日本時間に合わせる
price.index = price.index + timedelta(hours=14)

# 価格を円換算する
merged = pd.merge_asof(price, usdjpy, on="timestamp", direction="nearest") # 最も近いtimestampで結合
merged[price.columns[0]] = merged[merged.columns[1]] * merged[merged.columns[2]] # 円換算
merged = merged.set_index("timestamp") # timestampをindexにする
price = pd.DataFrame(merged[price.columns[0]]) # 円換算したものを新たに株価データとする

2. ポートフォリオの成績の計算

次に、毎月リバランスをしながら運用したときの成績を計算します。ロジックとしては以下です。

ポートフォリオの各資産クラスの割合を与える
初月に定められた割合で資産を購入する
毎月の月末に分配金を再投資し、保有株数×価格で各資産の割合を計算。割合が多い場合は売却、割合が少ない場合は購入する。保有株数とトータル評価額を記録しておく
最終日まで行ったら、トータル評価額の推移から年率平均リターン（複利計算）とリスクを求める

# 初期設定
amount = quotes.copy() # 資産毎の保有額
amount[:] = 0
num_of_stocks = amount.copy() # 資産毎の持株数
portfolio = quotes.iloc[0].copy() # ポートフォリオ比率
portfolio[:] = ratio
portfolio /= portfolio.sum() # ポートフォリオの正規化

# リバランス計算の実行
start_date = quotes.index[0]
for i in range(len(quotes.index)):
    d = quotes.index[i] # 日付

    # (1) 初月は定められた割合で購入。合計10,000円とする
    if d==start_date:
        purchase = portfolio * 10000 # 1万円分買う
        num_of_stocks.loc[d] = purchase / quotes.loc[d] # 購入株数
        amount.loc[d] = num_of_stocks.loc[d] * quotes.loc[d] # 評価額
    else:
        # (2) 毎月得た分配金は再投資。保有株数×価格で割合を計算し、目標と乖離がある場合はリバランスする。
        # 分配金再投資
        reinvested = dividends.loc[d] * num_of_stocks.iloc[i-1] # 該当月の分配金で再投資する金額
        num_of_stocks.loc[d] = num_of_stocks.iloc[i-1] + reinvested / quotes.loc[d] # 金額分の株数を追加購入
        # 各資産クラスの割合を計算
        amount.loc[d] = num_of_stocks.loc[d] * quotes.loc[d]
        ratio = amount.loc[d] / amount.loc[d].sum()
        # ポートフォリオとの差
        diff_portfolio = portfolio - ratio
        # 購入
        purchase = diff_portfolio * amount.loc[d].sum()
        num_of_stocks.loc[d] += purchase / quotes.loc[d]
        amount.loc[d] = num_of_stocks.loc[d] * quotes.loc[d]

amountが資産毎の保有額になりますので、これを使ってリターン・リスクを計算します。

年率平均リターンは、以下の式で求めます。*2

年率平均リターン = (最終日の株価÷起算日の株価)^1/年数-1

df = amount.sum(axis=1) # すべての資産の合計
days=(df.index[-1] - df.index[0]).days # データの日数
n=int(days/365) # 経過年
annualized_return = ((df.iloc[-1] / df.iloc[0])**(1/n) -1) # トータルリターンを年率換算

リスクは、月次株価の標準偏差を以下で求めます。*3

df = amount.sum(axis=1) # すべての資産の合計
monthly_return = df.pct_change(freq='1M') # 月次リターン
monthly_risk = monthly_return.std() # 月次リターンの標準偏差
annualized_risk = monthly_risk * np.sqrt(12) # 年率化

このポートフォリオの成績計算ロジックを各資産の保有割合ratioを入力として計算するよう関数化します。これにより任意のポートフォリオの成績を評価できるようにします。

3. 多目的最適化

最適化には、Pythonの多目的最適化ライブラリであるpymooを使います。まずpymooのクラスとして先ほど定義したポートフォリオ計算関数を評価関数とする最適化問題を定義します。

各資産の割合を表す変数の範囲は、0から1とする
- 変数の数は資産の種類=12とする
- トータルで100%になるように正規化してから評価する
- すべての変数が0にならないように制約を設ける
計算した年率平均のリスクの最小化とリターンの最大化を目的関数とする
マルチプロセスで並列評価する

from pymoo.core.problem import ElementwiseProblem
import multiprocessing
from pymoo.core.problem import StarmapParallelization

class PortfolioOptimizationProblem(ElementwiseProblem):

    def __init__(self, quotes, dividends, **kwargs):
        super().__init__(n_var=12,
                         n_obj=2,
                         n_ieq_constr=1,
                         xl=np.zeros(12),
                         xu=np.ones(12),
                         **kwargs)
        self.quotes = quotes
        self.dividends = dividends

    def _evaluate(self, x, out, *args, **kwargs):
        score = calc_rebalanced_result(self.quotes, self.dividends, x)
        f1 = score.iloc[0,0] # リスクの最小化
        f2 = -1 * score.iloc[0,1] # 平均リターンの最大化
        
        out["F"] = [f1, f2]
        out["G"] = -1.0 * np.sum(x) + sys.float_info.epsilon # 変数の合計が0超過の制約

n_proccess = 8 # プロセス数
pool = multiprocessing.Pool(n_proccess)
runner = StarmapParallelization(pool.starmap)
problem = PortfolioOptimizationProblem(quotes=quotes, 
                                           dividends=dividends, 
                                           elementwise_runner=runner)

次にアルゴリズムとパラメータを選択します。多目的最適化アルゴリズムには、遺伝的アルゴリズムの一つであるNSGA-IIを用い、個体数を50、世代数を100として、残りはデフォルトのパラメータとします。

また、各資産割合が100%となる変数を初期解に混ぜてみます。初期解を指定するには、アルゴリズム定義時にsampling引数を指定すれば良いようです*4。

from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.termination import get_termination
from pymoo.optimize import minimize
# NSGA-IIの設定
algorithm = NSGA2(
    pop_size=50,
    sampling=np.diag(np.ones(12)), # 初期解
    eliminate_duplicates=True
)
# 探索の実行
termination = get_termination("n_gen", 100)
res = minimize(problem,
           algorithm,
           termination,
           seed=1,
           save_history=True,
           verbose=True)

これで最適化を行いました。

4. 結果

上記の条件で最適化を実行したところ、手元のPCで大体3分で計算が終了しました。

最適化で評価したすべての解=ポートフォリオを青点・で、効率的フロンティアにあたるポートフォリオを赤丸●で示したのが以下の図です。数多くあるポートフォリオの中から最もリスクが低く、最大のリターンが出せるポートフォリオが最適化の結果として探索されていることがわかります。

効率的フロンティアの形状を見ると、リスク10%未満の範囲では、リスクに対してリターンがほぼ比例して高くなっています。ところがリスク10%を境にその傾きが低くなっているように見えます。この変化を理解しやすくするため、緑色の補助線を引き、リスクに対するシャープレシオを並べて記載した図を以下に示します。

傾きが緩やかであると言う事は、とっているリスクに対して見返りとして期待できるリターンが少ないということになります。実際シャープレシオもリスク10%を超えると低くなっています。そのためリスク10%程度のポートフォリオがシャープレシオが高いポートフォリオの中でも高いリターンの見込めるバランスの良いポートフォリオと言えると思います。

得られた効率的フロンティアのうち、(1)最もリターンの高いポートフォリオ、(2)リスク10%のバランスの良いポートフォリオ、(3)リスクが低め(5%)のポートフォリオの各資産の割合を見てみると、以下のようになりました。

資産	(1)最高リターン	(2)バランス	(3)リスク低め
現金	0.0	3.5	51.3
日本株	0.0	0.0	0.0
米国株	100.0	46.9	19.6
先進国株	0.0	0.0	0.0
新興国株	0.0	0.0	0.0
米国債券	0.0	0.0	0.6
先進国債券	0.0	0.0	0.1
日本REIT	0.0	0.0	0.0
米国REIT	0.0	0.0	3.0
先進国REIT	0.0	0.0	0.0
コモディティ	0.0	0.0	0.0
金	0.0	49.1	26.0

最もリターンの高いポートフォリオは米国株100%ですが、リスク10%のバランスの良いポートフォリオは米国株45%金55%となっており、リスクが低いポートフォリオはそれに現金を混ぜたものになっています。この結果は冒頭の記事の続きの記事で着目されている金50% 米国株50%のポートフォリオと近いものとなっており、最適化の結果としてはそれなりに妥当なものが得られていそうです。

最後に、探索の推移を、解分布の良さを表す指標Hypervolumeを使って見てみます。pymooではminimize()の引数にsave_history=Trueを入れておくと、各世代の変数・評価結果が得られるので、そこからHypervolumeを計算できます。参照点としてリスク最大・リターン最小の点を基準としました。

from pymoo.indicators.hv import HV
ref_point = np.array([np.max(F_all[:,0]), np.min(F_all[:,1])]) # 参照点はすべての解のリスク最大リターン最小の点
ind = HV(ref_point=ref_point)
hv = [ind(np.array([p.F for p in h.pop])) for h in res.history]
plt.plot(hv)

計算したHypervolumeを世代ごとに並べると以下の図のようになりました。今回50個体100世代=5000個のポートフォリオの評価をしていますが、20世代=1000個くらいの評価でおおむね効率的フロンティアの形状は見えていて、そこから少しずつ割合を調整して数値を上げているような推移になっています。もう少し世代数を増やしても改善できそうですが、実用上はこのくらいの評価でも十分いいポートフォリオが得られそうです。

5. まとめ

遺伝的アルゴリズムの1つであるNSGA-IIをつかって、12の資産クラスからなるポートフォリオのうち、低いリスクで高いリターンが得られるポートフォリオの集合（効率的フロンティア）の探索を試してみました。10年間の実績値を使った結果、米国株・金を中心としたポートフォリオが良好な成績を示すことがわかりました。

ただ、こういったポートフォリオをそのまま使えるかというと微妙で、米国株・金に偏っており分散されていないので、米国株・金が同時に下がるようなシチュエーションには（この2つは相関が低いようなので可能性は低そうですが、もしあれば）弱いでしょう。最低限他の資産にも分散させる、といった制約も考えて最適化したほうがよさそうです。

それから、今回考慮しませんでしたが、リバランスには売買手数料・為替手数料や税金がかかってきますし、今後の市場見通しやその不確実性などの情報も考慮しておく必要があると思います。時々刻々と変わってくる株価・市場に対して最適なポートフォリオも変わってくるので、このような最適化方法を使うにしても、最適化自体やポートフォリオを適宜更新するような運用も同時に考える必要がありそうです。

ここまで考えられているのかはわかりませんが、いくつかの証券会社では、ポートフォリオを提案してくれて、自動で運用するものがあるようです*5 *6 *7。費用が安ければ、こういったサービスを使うのが個人では容易なように思いました。

*1:銘柄については、正直なところ妥当かどうか自信がありません。もしより良い銘柄があれば教えていただきたいです。

*2:リスク・リターン（2）異なる期間の収益率を⽐較する「年率換算」

*3:標準偏差・分散｜証券用語解説集｜野村證券

*4:Biased Initialization

*5:最適ポートフォリオって何ですか？ | ロボアドバイザーならWealthNavi（ウェルスナビ）

*6:マネックスアドバイザーサービスの特長 | マネックスアドバイザー | マネックス証券

*7:投信工房 - ロボアドバイザー | 投資信託 | 松井証券

2022-05-09

Modelica Buildings LibraryモデルのOpenModelicaでの実行・FMU出力可否の状況

以前，以下のようなModelicaに関する記事を書きました．
ModelicaのインストールとModelica Buildings Libraryのサンプルモデルの実行
その後，JModelicaとBuildingsライブラリを使ったモデル開発をしており，その状況をModelicaライブラリ勉強会で紹介させていただいたりしました．

しかし，JModelicaでは以下の点がモデル開発のネックとなっていました．

JModelicaはオープンソースとしての開発が終了しており，無償での最新の開発環境が入手できない
最後のオープンソース版のJModelicaでは，今後Buildingsライブラリが更新された場合に動作するかわからない
JModelica自体にはグラフィカルなモデル作成インターフェースがなく，モデル修正→実行→デバッグが手間

そんな折，finbackさんより，OpenModelicaでもBuildingsライブラリのうち実行できるモデルが増えているようだ，という情報をいただきました．具体的には，以下のOpenModelicaでのtest結果でSimulateできている項目が多くなっている，とのことです． Buildings_latest test using OpenModelica
そこで，BuildingsライブラリのExamplesのうち代表的なものが，OpenModelicaでも実行できるか調査してみました．また，私の場合は最終的には開発したモデルと他のアルゴリズム等の連携による適切な運転設定の探索などをしたいと考えているため，FMUの書き出し・実行ができるかも調査しました．

実行環境

Windows 10 21H2
OpenModelica v1.20.0-dev (Nightly build, 64-bit)
Modelica Standard Library 4.0.0
Modelica Buildings Library (master, 最終コミットFeb 16, 2022)
Python 3.8.7 (FMPy実行用)
FMPy 0.2.27

調査内容

OpenModelica上でBuildingsライブラリモデルの動作を確認する
OpenModelicaでFMU書き出し→FMPyでの実行を確認する
FMU書き出し条件は以下
- バージョン：2.0
- Platforms: Native
- Model Description Filters: internal
- Include Modelica based resources via load Resourceにチェック

結果

調査結果は以下の表のようになりました．MEはModel Exchange, CSはCo-Simulationです．

モデル名	OpenModelicaでの実行	FMU書き出し	FMUの実行(CS)	FMUの実行(ME)
Buildings.Examples.Tutorial.Boiler.System1	○	○	○	○
Buildings.Examples.Tutorial.Boiler.System2	○	○	×	○
Buildings.Examples.Tutorial.Boiler.System3	○	○	×	○
Buildings.Examples.Tutorial.Boiler.System4	○	○	×	○
Buildings.Examples.Tutorial.Boiler.System5	○	○	×	○
Buildings.Examples.Tutorial.Boiler.System6	○	○	×	○
Buildings.Examples.Tutorial.Boiler.System7	○	○	×	○
Buildings.Examples.Tutorial. SpaceCooling.System1	○	○	○	○
Buildings.Examples.Tutorial. SpaceCooling.System2	○	○	×	×
Buildings.Examples.Tutorial. SpaceCooling.System3	○	○	×	×
Buildings.Examples.ChillerPlant. DataCenterContinuousTimeControl	○	×	-	-
Buildings.Examples.ChillerPlant. DataCenterDiscreteTimeControl	○	×	-	-
Buildings.Examples.ChillerPlant. DataCenterRenewables	○	×	-	-
Buildings.Examples.VAVCO2. VAVSystemCTControl	○	○	×	×
Buildings.Fluid.Chillers.Examples. ElectricEIR_AirCooled	○	○	×	○

まず分かったのは，OpenModelica上でのシミュレーションは複雑なモデルでもかなりできるようになってきている，ということです．これは朗報で，Buildingsライブラリを使ったモデルの開発・検証自体はOpenModelicaでも十分やっていけそうだと感じます．

一方，モデルによってはFMU書き出しができない・FMU書き出しができてもFMPyでは実行できない条件があることもわかりました．

CSではFMUとして書き出しできても実行できず，MEで書き出しすると実行できるものがありました．これは，OpenModelicaがFMUに埋め込んだソルバの実行条件がうまく設定されておらず，FMPy側のソルバならうまくいくのが要因のようです．ちなみに，FMPyのソルバはOpenModelicaが埋め込むものと同じCvodeが選択できます．CSとMEどちらも動いたモデルで見る限りは，実行結果に大きな差異は出ないようでした．
Model Exchangeで書き出しをしても，モデルが実行できないものがありました．この場合，Output Intervalを整数としたり，ステップサイズをVariable-step→Fixed-stepとして長めにとったりすることで実行できるものがありました．実際，BoilerモデルではOutput Intervalを300に，SpaceCoolingでは22にすることで実行できました．
Buildingsライブラリでは気象データの取り込みにweaDatが多用されるのですが，この時気象データの指定パスを "modelica://"から始まるパスとすると正しくFMU書き出しできないため，絶対パスで"C:/hoge/huga.mos"などと指定する必要がありました．
理由が不明ですが，WetCoilやCoolingTowerの熱交換器モデルはFMU化すると実行できないようです．そのため，現状はこれらのモデルを使わないか，代替モデル(DryCoil系など)を使う必要があります．

いろいろ制約があるものの，ある程度のモデルはFMUとして書き出し・実行できていることから，今後はJModelicaに代わってOpenModelica上でモデル開発を進めていこうと思います．

2021-08-11

社会人博士課程を振り返って

社会人博士博士課程

私は大学院を修士で修了したあと企業に就職し，その後大学院の社会人博士課程で博士号を取得しました．

学位取得直後は会社業務がピークに来ていたり，博士号以外の資格試験準備で時間がなく，気づいたら1年近くも経過してしまったのですが，社会人博士課程で博士号を取得するまでを振り返ってみたいと思います．

社会人博士課程に入学した理由

前述のとおり，私は修士卒で就職しています．修士のときの指導教官からは博士課程へも誘われていたのですが，以下のような理由から就職を選択しました．

やりたいことが研究というよりは開発寄りであった
私自身は，学部～修士の研究をそれなりに楽しく取り組み，ある程度成果を出せたと思います．しかし，それは研究ができたからというよりは，研究をするために発生するものづくりや開発のプロセスを得意としていたからこそであり，就職して開発をするほうが自身のキャリアとしてはよいと考えていました．
なぜ学部で就職しなかったかというと，当時は理工系の学生は学部よりも修士卒の学生のほうが就職で有利であり，配属先も学部生よりも優遇してもらえることが多い，ということを聞いていたためです．
博士課程進学後の学費と就職の不安があった
博士まで進学することに対して学費の不安がありました．また，学部・修士の研究室は博士課程の学生がおらず，博士号取得後のキャリアパスがあまり見えなかった点も大きかったと思います．企業は博士にあまり門戸を開いていないという印象でしたし，大学教員になろうにも，国立大学の独立行政法人化などでテニュアのポジションが減っている情勢で，就職の不安がありました．合同でゼミをやっていた隣の研究室には博士課程の先輩もいましたが，大変そうにしていたように見えました．
博士課程の研究環境に不安があった
学部・修士の先生は，あまり熱心な指導をされる方ではなく，研究分野の知識や研究に必要なスキルがその研究室にいても十分には身につかない不安がありました．他大学の博士課程へ行くこともできましたが，それよりは就職のほうが有力な選択肢でした．

このように一旦就職したものの，再度博士号取得を考え始めたのは以下の理由でした．

キャリアの不安が解消された
就職した企業の配属先が研究所であり，博士号取得に理解がある(むしろ推奨されている)環境でした．また，研究所の上司・先輩には博士号取得者がおり，キャリアとしても博士号を取る選択肢が現実的になりました．
経済的な不安が解消された
博士課程では，国立大でも年間50万円ほどの学費＋入学料が必要です．就職したことで，これらの学費は賄えるだけの収入ができました．さらに業務に関わるテーマであれば会社から学費補助と博士号取得時の奨励金が出ることもあり，経済的な問題が解消されました．
研究したいテーマができた
配属先が研究所ということもあり，会社業務で応用寄りの研究を行っていました．修士までとは異なる分野でしたが学会誌論文を1本書くことができ，この分野で研究を進められる自信がつきました．また，学会誌論文の内容の発展で博士研究として取り組めそうな研究テーマができました．

社会人博士課程のテーマと研究室の選定方法

社会人になってから博士号を取ろうとすると，通常，テーマと研究室に悩むのではないかと思います．私の場合は，上述の通り就職後に研究を始めた分野でテーマを決めました．研究室は，就職後に業務として参加していた学会で知り合いになった先生のところでお世話になりました．

職場には，私と同様に業務に関連した分野で新たに研究室を選ぶ方もいましたし，業務とは離れるけれども修士まで在籍していた研究室に戻って研究をされている方もいました．

私はあまり修士までの分野にはこだわらないでテーマを決めましたが，博士課程でも修士研究と同じ分野の研究をしたいという方もいると思います．研究を進める上でなかなか成果が出ないこともありますので，自身の興味や研究に対するモチベーションが保てるテーマにするのが良いかと思います．

入試

入学までは以下を行いました．

入学半年前頃まで：上司・先生への入学相談
入学3ヶ月前：出願
入学2ヶ月前：入学試験

私の所属していた大学では社会人入試という枠があり，書類提出と口述試験のみで入学試験が行われました．口述試験では，専攻の先生方3名に対して志望経緯や研究背景，これまでの研究成果，博士課程で取り組む予定の研究内容と大まかなスケジュールについて説明し，質疑応答をしました．

一般入試では英語のTOEICスコアによる評価もありましたが，社会人入試ではありませんでした．
一方，社会人入試では，研究業績一覧，研究計画書と，所属企業の就学承諾書の提出が必要でした．上述のように，会社からの学費補助の申請の必要もあったため，事前に上司にも相談して就学の承諾を得ていました．

学位取得までの大まかな流れ

私は2017/10～2020/9の3年間をかけて学位を取得しました．大まかに以下の流れで学位取得に至りました．

1年目
- 2017/10～2017/12 週1回の集中講義受講
- 2017/12 国内学会での発表(1回目)
- 2018/1 国際会議への論文投稿(1回目)
- 2018/7 国際会議での発表(1回目)
2年目
- 2018/12 国内学会での発表(2回目)
- 2019/1 国際会議への論文投稿(2回目)
- 2019/4 国内論文誌への論文投稿(5月に結果通知：条件付き採録，6月に回答書・訂正論文投稿し採録)
- 2019/4～7 学会のコンペティションへの参加
- 2019/6 国際会議での発表(2回目)
3年目
- 2019/12 国内学会での発表(3回目)
- 2020/1 国際会議への論文投稿(3回目)
- 2020/2 国内研究会で発表
- 2020/3 博士論文の執筆開始，共同研究者への論文利用承諾書の作成依頼
- 2020/5 予備審査
- 2020/7 国際会議での発表(3回目)
- 2020/8 本審査
- 2020/9 学位取得．課程修了．
修了後
- 2020/10 博士論文の製本，審査委員の先生方へ送付
- 2020/12 国内学会での発表(5回目，博士論文記載の内容の発展版)

在学期間中はおおよそ月1回の頻度で進捗打ち合わせと，学会発表前にはゼミで発表練習などを行いました．

また，社会人博士課程といえど博士論文の研究以外に単位が必要だったため，講義を受ける必要がありました．平日は会社業務があるため，土曜午後に開講していた集中講義を1年目に受講しました．

研究を進める環境と1日の過ごし方

私は，研究活動をほとんど自宅で行っていました．これは，以下の理由がありました．

研究テーマが自宅のPCでもできるものだった
研究室にはVPNで接続でき，必要なときに計算資源や論文誌へアクセスできる環境があった
先生が社会人博士に対して理解がある方で，基本的には通学しなくともメール・SlackやSkypeでのやり取りで研究が進められるように配慮いただけた

特に新型コロナウィルスの影響が及ぶ前からテレワークで研究ができており，コロナによる影響はあまり受けずに済みました．この点はとても恵まれていたと思います．

しかしながら，平日昼間は仕事で時間が使えず，時期によってはそれなりに残業もあったことから，仕事以外の時間でなんとか研究の時間を確保するようにしていました．平日は専ら以下のような感じで過ごしていました．

6:30 起床
7:00～9:00 通勤．座れないときは論文を読む．座れればプログラム開発や論文執筆
9:00～18:00 会社業務
18:00～20:00 通勤．座れる電車を選んでプログラム開発や論文執筆
20:00 帰宅．食事の後に研究．先生との打合せがある場合はWeb会議など

休日も家事を除いて多くの時間を研究に費やしていました．リフレッシュの時間も取るようにはしていましたが，たとえば土日のうちどちらかは完全に休日にするなど，もう少し多めに休みを入れたほうが良かったかと思っています．

苦労したこと

博士課程では苦労の連続でした．

ライフイベント
入学前に結婚し，在学中に子供ができました．そのため，研究に費やせる時間が独身のときと比べてかなり少なくなっていました．家族に家事などを協力してもらったり，子供が寝ている時間などを使って，研究の時間を捻出していました．
異動
入学直前に研究所から事業部門への異動があり，業務としては研究から離れる事になりました．異動先では働きながら大学院に行くのはあまり例が無いようだったので，就学を許可してもらえるよう，異動前・異動先の上司に相談をしました．また，勤務先が遠方になったため，通勤時間が増加しさらに時間が減ってしまいました．研究に必要な時間を確保するため異動先の方との飲み会などの交流も減らさざるを得ず，心苦しい思いをすることもありました．
一方で，異動によって良い面もありました．私の研究は応用研究でしたので，事業部門に移ることで応用先のドメイン知識を得ることができました．この知識は研究を進める上で有利な方向に働きました．
実験しても良い結果が得られない
私の研究テーマは，進化計算アルゴリズムの実問題への応用に関するものでしたが，既存手法の改良による性能向上がなかなかできず苦戦していました．改良手法の評価が早く終わるのであれば，トライアンドエラーを繰り返し様々な改良手法のアイデアを試すこともできるのですが，私の扱っていた実問題のシミュレータは計算コストが高く，1つのアルゴリズムの評価結果を得るのに1日以上，下手をすると2週間近くもかかるため，試行数を増やせない状況でした．長期間かけて実験しても良好な結果が得られないことも多かったです．

社会人博士課程で得たもの・良かったこと

苦労したことは多かったものの，それ以上に得たものや良かったことが多くありました．博士課程に進学してよかったと心底思っています。

専門性が高められた
博士では，修士までとは異なる研究分野であったため，大学院では基本的な内容から勉強をし直しました．また，当該研究分野のトップ研究者である先生に師事することができ，課題に対して議論をしたり助言を得ることができました．おかげで研究成果にも恵まれ，トップカンファレンスに採択されて発表・議論することもできましたし，学会の論文賞をいただくこともできました．こういった研究活動を通して，高い専門性を得られたと思います．
論理的思考力が鍛えられた
研究を行う上で，有意義な課題を設定する，課題に対する解決策を立案し検証する，成果を道筋立てて説明する，といったことを行います．これらを論理的に考えて実行する力が養われました．
文章力・英語力・プレゼンテーションスキルが鍛えられた
博士課程では，査読付き論文1編，査読付き国際会議発表3回，国内の学会発表も4回とそれなりの数の論文を書き，発表をしました．また，博士論文と予備審査・本審査時の発表資料では，指導教官の先生からの徹底的なレビューと審査委員の先生方からの指摘を受け，修正や説明の方法を考えました．こういった経験から，文章力や英語力，プレゼンスキルがかなり鍛えられたと思います．
ソフトウェア開発力が向上した
私の研究対象であった進化計算の分野では，深層学習分野でいうTensorFlowやPyTorchのようなフレームワークを用いて，既存のアルゴリズムとの比較や新規アルゴリズムの実装を効率的に行うことがトレンドとなっていました．
これまで私はすべて自作のソフトウェアでアルゴリズムの実装を行っていましたが，フレームワークを使った実験ができるように自作ソフトウェアを改良開発しました．よく抽象化され設計されたフレームワークをカスタムして実装・実験したり，他の研究者の書いたコードを参考にすることで，ソフトウェア設計に関して理解を深める事ができました．また，ソフトウェア開発に必要なツールについての知識も得ることができました．
私は学部・修士と電気電子系の学科であったため情報系の知識が不足しており，実際に開発しながら学ぶことができたのは良い経験でした．
チームで分析・開発する経験が得られた
1年目の集中講義では，ある通販会社の購買データをチームで分析し，分析結果をもとに販売戦略を立てるという作業をしました．また，2年目には学会の開催するコンペティションに，研究室のメンバとチームを組んで参加しました．チームで同じ目標に向かって議論しながら分析や開発をすることで，相手にわかりやすくなるようなコミュニケーションをとったり，チームの意見や成果をまとめたりする経験が得られました．
タイムマネジメント力・根性がついた
職場での業務と博士課程の研究を並行していたため，少ない時間をやりくりしてなんとか研究をこなすための時間の使い方や，投稿締め切り前などここぞというときに集中して研究を進める根性が身についたと思います．

おわりに

現在は，研究所に戻って研究開発業務を行っています．博士号自体が直接業務に役立つことはいまのところありませんが，博士号を取得する過程で得られた能力や経験によって，地力がついたと感じることが多々あります．また，幸いにも指導教員の先生と研究を続けさせていただくこともできています．

今後も研究を続け，会社の製品・サービスはもとより，社会に役立つようなアウトプットが出せるように努力していきたいと思います．

以下に私の博士論文と公聴会で使用したスライドを公開しています．どなたかの参考になれば幸いです．

2021-07-22

第二種電気工事士技能試験を受験しました

第二種電気工事士

タイトルの通り，令和三年度上期の第二種電気工事士の技能試験を受験しました．技能試験で必要となる複線図の書き方や部品などの取り扱いは一から勉強しましたので，受験にあたって取り組んだことや気になった箇所などをまとめておきます．

受験の経緯
購入した工具・部材・参考資料
- 工具
- 部材・参考資料
練習方法
候補問題を解いていて気になった点
受験日当日
おわりに

受験の経緯

電工試験を受けようと思った理由は大きく3つあります．

将来的に自宅内の工事をやりたい
私は趣味で電子工作をしているのですが，電工の資格があると取り組める幅が広がるなと思うことが多々ありました．これまで電子・情報系の学部・学科だったため強電系はあまり扱って来なかったのですが，将来的には資格をとってみたいと思っていました．
業務でも取得しておいたほうが有利に働く場面がありそうだと思った
私自身は研究開発職であり，電気工事自体を行うことは基本的にはありません．しかし，電気工事を含む業務の改善やそれに必要な技術開発をしたり，電気工事を伴う工事の発注などを行う場面があるため，知識と資格を持っておくこと自体は必要だと思っていました．
電験三種に合格していたため，第二種電工試験の筆記が免除となった
大学生の時に技術士（電気電子部門）一次試験に受験して合格していました．この試験と出題範囲が重なる電験三種（第三種電気主任技術者試験）の資格取得を優先しており，昨年無事試験合格し免状取得しました．第三種電気主任技術者免状の取得者は第二種電工試験の筆記試験が免除され，受験するのは技能試験だけで済むため，ここで受けない手はないと考え受験しました．

購入した工具・部材・参考資料

工具

学生時代に電工試験を受験していた友人に工具セットを借りることができたので，ほとんど新規で工具を購入する必要はありませんでした．唯一，借用した工具セットにストリッパが含まれておらず，ストリッパは40分しかない技能試験での利用は必須だと思っていたので買い増ししました．
www.hozan.co.jp

通常は，以下のような工具セットを購入するのが良いのかと思います．実際，試験会場でもこのセットを使っている方がかなり多くいました．
www.hozan.co.jp

部材・参考資料

練習用部材と参考資料として，以下の部材セットのうち，3回練習用セットを購入しました．
www.hozan.co.jp

本セットには技能試験対策DVD，およびハンドブックが付属しており，他に参考書の購入や講習会の受講などをせずとも済みそうだという期待がありました．
これは期待どおりで，本セットに付属のDVDとハンドブックに加え，過去問・HOZANの電工試験対策HPを見るだけで十分な知識をつけ，練習をすることができました．
また，部品がなく十分に練習ができない状態で本番に臨むのは避けたかったので，3回セットを購入しました．しかし，実際には練習にかけられる時間がとれず，各候補問題2回ずつ練習した状態で本番の受験となったので，1回分は不要でした．
もっというと，この部材セットは結構余裕をみて材料が入っているため，2回ずつ各候補問題をやったものの材料消費量は1.5回分程度でした．1回分セットを買って，1回プラスαの練習をし，練習が足りない場合に部材セットを追加する形でも良かったかもしれません．

練習方法

基本的には，DVDを見て手順を確認し，候補問題をNo. 1から13まで順に作成する，という作業を2周行いました．1周目はそもそも扱う部品自体を初めて目にする場合があり施工方法を学んだり，ケーブル被覆や芯線被覆を剥く長さを覚えたりしていたため時間がかかっていました．
2周目になるとだいぶ習熟してきて，作品の完成度向上やなるべく短時間で完成するような工夫，失敗したときのリカバリの練習を主にしていました．1周だけだと不安があったので，本番前に2周できたのは良かったと思います．
また，一周目の途中から，以下のように完成した回路の写真とミスした/しそうだった点・作成時間などを記録するようにしたのですが，これが割と回路作成のペース配分の検討や注意点の見直しに効果的でした．

第二種電気工事士技能試験の練習、候補問題No.12をやった。見直し含めて34分で完成。途中PF管を25と19どちらに接続するか迷ってしまった。サイズ的には19なのだけどボックスコネクタにType25のような表記があったので。確認しておこう。 pic.twitter.com/Q1RbAT5U1e
— Yoshihiro Ohta (@ohtayoshihiro) July 17, 2021

候補問題を解いていて気になった点

上記練習用部材セット付属のDVD・ハンドブックで練習していたのですが，不明瞭な点があって調べたので，いくつか記載します．素人目にはこのあたりがわかりにくいように感じましたが，他の参考書と併用していればとくに気にならなかったのかもしれません．

配線材料が不足していた場合の対処
例えば，上記の参考書では，候補問題No.3ではVVF 1.6mm-2Cのケーブルが1550mm配布される想定でしたが，ケーブル外装を取り除く長さを考慮すると1570mmほど必要で，材料が不足することになります．実際の試験でこのような問題に遭遇した場合どうすべきかわかりませんでした．
→結果から言うと，配線長さの欠陥条件は規定値の半分以下ということなので，どこか多少不足している箇所があっても配線長さが半分以下にならなければ欠陥ではないとのことでした．
ちなみに，本番ではちょうど候補問題No.3が出題されましたが，VVF 1.6mm-2Cケーブルは1650mm配布されたので，不足するということはありませんでした．
圧着をミスして圧着し直した場合，接続に必要な寸法を下回るが欠陥になるのか
→これは，1. に記載したものと同様，所定長さの半分以上を確保できていれば欠陥になりません．したがって，圧着をミスしてその箇所を切断し，少し短くなった配線材料で再度圧着したため，結果として寸法が短くなっても，欠陥にはなりません．
欠陥がいくつかあると不合格なのか，欠陥が1つでもあると不合格なのか
→欠陥が1つでもあると不合格となります．
3極スイッチ，4極スイッチの端子が二つある側に極性はあるのか
→極性はありません．回路が正しければどちらに接続されても良いとされています．3極スイッチの0番端子には非接地（黒線）が来ることが決まっているので，残りの1番と3番は白色・赤色どちらが接続されても良いことになります．私は「3番端子に白色を接続するようにする」というマイルールを定めて，複線図作成と回路の作成をしていました．
候補問題No. 12でPF管のボックスコネクタが大きい穴・小さい穴どちらにも取付可能だが，取付箇所をどう判断するか．
→アウトレットボックスの穴の径は以下3つとなっています．
- G16(C19)＝21.5mm
- G22(C25)=27.1mm
- G28(C31)=34mm
PF管ボックスコネクタはこのうちG16, G22(22mm, 27mm)兼用となっている様なので，どちらでも接続可能です．試験でアウトレットボックスに空いている穴か，施工条件から判断します．

受験日当日

試験地は神奈川だったのですが，大学の教室で試験が行われました．机が普通の講義用の机で作業スペースが狭めだったのでちょっと驚きましたが，自宅の練習スペースもそこまで広いわけではなかったので，特に問題ありませんでした．
試験前には，配布部材に不足がないか確認する時間が設けられたので，そこで候補問題のどれかを特定し，頭の中で複線図や回路作成作業手順をイメージすることができました．出題された候補問題No. 3は記憶に残っていたのと，苦手な問題ではなかったため，試験開始後は苦もなく作成でき，28分ほどで完成，入念な見直しをして35分で概ね完成しました．

おわりに

候補問題の練習にかなり時間を使ったので合格していてほしいですが，欠陥を見落としている可能性もなくはないので，合格発表の8/20を待ちたいと思います．不合格だとショックなのですが，まあその時は下期試験を受けなおすことにします．

2021-04-25

ModelicaのインストールとModelica Buildings Libraryのサンプルモデルの実行

Modelica ModelicaBuildingsLibrary Python

私は研究でビルシミュレーションを行うため，EnergyPlusというソフトウエアを利用していました．
EnergyPlusはビル全体のエネルギーや室内環境をシミュレーション可能で，テキストファイルで入出力をするシンプルな構成だったのですが，建物モデル作成が煩雑であったり，特殊な制御ロジックのシミュレーションが困難でした．

ところで，マルチドメインの物理シミュレーションが可能なモデリング言語として，Modelicaというものがあります． Modelicaのライブラリの1つである，Modelica Buildings Library(MBL)を使うと，建物の水配管や室などのモデルをMATLABのSimulinkライクにグラフィカルに接続してシステムモデルが作成できるほか，制御ロジックのカスタムとシミュレーションも可能なようです．
今回はModelicaを使ってMBLを実行する環境を構築し，サンプルモデルを実際に実行してみます．

1. ModelicaとMBL

1.1 Modelica

Modelicaでは，微分代数方程式のような形で物理モデルを記述してモデリングし，そのモデルの実行・シミュレートができます．また，作成したモデルをモジュールとして扱い，モデル同士を接続してグラフィカルにシステムをモデリングすることも可能です．Modelica自体はオープンソースソフトウエアで，開発環境としては，有償・無償のものがそれぞれ存在します．無償のものとしてはOpenModelicaがよく用いられるようです．

Modelicaの基本的な操作は以下によくまとまっています．
OpenModelicaのチュートリアル資料

その他のOpenModelicaに関する日本語の資料として，以下の勉強会のアーカイブが参考になります．
Modelicaライブラリ勉強会 - connpass

1.2 Modelica Buildings Library(MBL)

MBLは米国ローレンスバークレー国立研究所(LBNL)の開発するオープンソースのModelicaライブラリです．MBLについては，以下の資料が詳しいです．
Modelica Buildings Library の紹介

この資料にも記載があるのですが，MBLは，有償のDymolaという開発環境と無償のJModelicaという2つでしかテストされていないようで，残念ながらOpenModelicaでは使えない部分があります．
以下のページに，テストによるエラー・警告発生とシミュレーション結果ログが記載されています．
Buildings_latest test using OpenModelica

Tutorialは一応OpenModelicaでシミュレーション可能となっていますが，私の環境では実行できないものがありました．そこで，JModelicaを使います．

1.3 JModelica

JModelicaについて，以下の記事が詳しいです．
JModelica.org 使い方メモ | Finback

JModelicaはModelon社の開発したソフトです．元々オープンソースだったようですが，オープンソースとしての開発が終了しています．最後のオープンソース版のJModelicaソースコードおよびバイナリがHPからダウンロード可能です．

しかし，JModelicaにはOpenModelicaにあるようなグラフィカルなモデル作成インターフェースがないようです．そこで，モデル作成はOpenModelicaで行い，JModelicaでFMU(Functional Mock-up Unit)化し，他のシミュレータからでも呼び出し可能な形にして実行します．
実行には，フリーのFMUシミュレーション用PythonライブラリであるFMPyを使うことにします．

2. 開発環境の構築

2.1 バージョン

Windows 10 20H2
OpenModelica v1.17.0 (64-bit)
JModelica 2.14
Modelica Standard Library 3.2.2 (JModelica付属)
Modelica Buildings Library 6.0.0
Python 3.8.7 (FMPy実行用)
Python 2.7.13 (JModelica付属)

2.2 インストール

OpenModelicaを以下からダウンロード・インストールします．
https://openmodelica.org
JModelicaを以下よりダウンロードしてインストールします．
https://jmodelica.org
FMPyの公式HPの通りインストールします．
python -m pip install fmpy
ただし，私の環境ではそのままでは動作せず，別途PyQt5, pyqtgraphもインストールが必要でした．
python -m pip install PyQt5 pyqtgraph

2.3 Modelicaライブラリのフォルダパス

MBLは，Modelicaの標準ライブラリである，Modelica Standard Library(MSL)に依存しています．JModelicaのインストールディレクトリ\install\ThirdParty\MSLに，MSLが格納されています．今回はこちらを用います．
MBLは，OpenModelicaのインストールディレクトリ\lib\omlibraryに格納されています．
OpenModelica v1.17には6.0.0と7.0.0の2つが付属していましたが，今回は6.0.0を使います．(JModelica付属のMSLは古いようで，Buildingsライブラリも最新だと動かないようです．)
同じディレクトリに2つのBuildingsライブラリが混在しているとJModelica側のパス設定が大変なので，Buildings 6.0.0フォルダをJModelicaのインストールディレクトリ\install\ThirdPartyに移動しておきます．

3. サンプルモデルの実行

3.1 モデルの確認

OpenModelica Connection Editor(OMEdit)を開きます．自動でMSLの最新版(3.2.3)と付属ライブラリ(Modelica ,ModelicaServices, Complex)がロードされますが，アンロードしておきます．
ファイル→ライブラリのロードから，JModelicaのMSLとMBL 6.0.0をロードします．
MBLに付属するサンプルモデルを開きます．今回は，Buildings.Examples.Tutorial.Boiler.System6を用います．

f:id:ohtayo:20210425215752p:plain — OpenModelica Connection Editorで表示したBoiler.System6のモデル

このモデルは，ボイラ，熱交換器，空調対象室が1つずつあり，以下の制御が含まれています．
- 室温および外気温に応じたボイラ・熱交換器のポンプ発停制御
- 室温を一定にするための熱交換器入口温度制御
- ボイラ安定運転のためのボイラ入口温度制御
また，外気温データを外部ファイルから読み込み，それを用いて室温変動と制御シミュレーションを行います．

このモデルを変更する際は，System6モデルを右クリック→複製して使用します．Buildingsライブラリ以外の場所においた場合，weaDat, TOut, TSetSupのBuildingsというプレフィックスが省略されているので，コードに付け足す必要があります．

3.2 JModelicaによるFMUの作成

モデルを確認，あるいは変更したら，これをJModelicaでFMU化します．以下の様なPython スクリプトを作成して，JModelicaインストールフォルダにあるIPython64.batで実行します．

# coding: utf-8
# Modelicaライブラリのパスを追加
import os
os.environ['MODELICAPATH'] = r'C:\JModelica.org-2.14\install\ThirdParty\MSL;C:\JModelica.org-2.14\install\ThirdParty'
# FMU化
from pymodelica import compile_fmu
fmu = compile_fmu('Buildings.Examples.Tutorial.Boiler.System6',
                  version='2.0',
                  target='cs',
                  compile_to="Buildings_Examples_Tutorial_Boiler6_v20_cs.fmu",
                  compiler_log_level='w,i:log.txt')

すると，指定したファイル名の.fmuファイルが作成されます．

3.3 FMUの実行

Python3環境でFMPyを実行します．
python -m fmpy.gui
現れたGUIに，さきほど作成したFMUファイルをドラッグドロップします．

f:id:ohtayo:20210425215924p:plain — FMPyのモデルロード結果

右三角のStart Simulationボタンを押すとシミュレーションを実行します．
例えば，熱交換器入口温度設定値TSetSup.y, 熱交換器入口温度計測値temSup.T, 室温vol.T, 熱交換器温度制御用三方弁開度valRad.yにそれぞれチェックをつけると，以下のグラフが見れます．

f:id:ohtayo:20210425215929p:plain — 熱交換器入口温度の制御結果グラフ

外気温に依存して設定値が変動し，それに応じてバルブ開度が変わって熱交換器入口温度計測値と室温が変化していることがわかります．

4. まとめ

OpenModelica, JModelicaとModelica Buildings Libraryを使って，ビル熱源系のサンプルモデルのシミュレーションを実行しました．
ModelicaはMATLAB/Simulinkに使い慣れている私からするとすこしクセがありますが，MBLだけでなくMSLにも豊富な各ドメインのライブラリがあり，非常に有用だと感じました．Examplesもそれなりにあるので，これらを参照しながらいくつかモデル作成をこなしていけば使い方には慣れてくるのではないかと思います．
また，FMPyではFMUモデルの入力を変更しながら1ステップずつ実行などができるようなので，このあたりも試してみたいと思います．

2021-02-17

Juliaで粒子群最適化

Julia 粒子群最適化単目的最適化

最近，Juliaを少しずつ学んでいます．勉強がてら，粒子群最適化(Particle Swarm Optimization, PSO)のアルゴリズムを実装し，Rosenbrock関数の最適化を行ってみました．

github.com

0. 開発環境

Windows 10 20H2
Julia 1.5.3
Visual Studio Code + Julia extension

1. 目標とコンセプト

今回の目標は，単目的の粒子群最適化アルゴリズムを実装し，Rosenbrock関数を最適化することです．ただし，将来的に多目的問題・制約付き問題に拡張可能とするための考慮をします．
以前もブログに記載しましたが，私は普段Javaの多目的最適化ライブラリjMetalを用いており，多目的に拡張する箇所などはjMetalを参考にします．
ただし，JuliaはJavaのようなクラス構造を持つ言語ではなく，型の定義と多重ディスパッチによって振る舞いを変更する特徴を持つ言語なので，そこは言語に合わせた実装とするよう努力します．

2. 問題定義

最適化問題として，今回はRosenbrock関数を対象とします．

function rosenbrock(x)
    [sum([ 100(x[i + 1] - x[i]^2)^2 + (x[i] - 1)^2 for i in 1:length(x) - 1])]
end

将来的に多目的最適化でも使うため，出力はFloat64ではなくVector{Float64}としておきます．
それから，以下の工夫をするコードを追加します．

rosenbrock以外の問題にも対応できるよう，evaluateという別の関数を経由して評価関数を実行する
rosenbrockは連続値最適化問題だが，それ以外の種類の問題ではevaluate関数の動作を変更できるよう，問題を複合型とする
複合型に，変数や目的関数の上下限値，変数の数などのパラメータをもたせる
変数・目的関数値を上下限値を用いて正規化して最適化する

# 抽象型
abstract type AbstractProblem{T <: Number} end
# 具体型
mutable struct DoubleProblem{T} <: AbstractProblem{T}
    objfunc
    numobjs
    numvars
    confunc
    numcons
    limit::Matrix{Float64}
    bound::Matrix{T}
end
# コンストラクタ
function DoubleProblem(objfunc; numvars=1, numobjs=1, confunc=NaN, numcons=0, limit=fill(NaN, 1, 1), bound=repeat([0.0 1.0], numvars))
    DoubleProblem(objfunc, numobjs, numvars, confunc, numcons, limit, bound)
end

# 評価関数
norm(value, min, max) = (value - min) ./ (max - min) # 正規化
denorm(value, min, max) = value .* (max - min) + min # 非正規化
function evaluate!(problem::DoubleProblem, solutions::Vector{<:AbstractSolution})
    for s in solutions
        vars = denorm(s.variables, problem.bound[:,1], problem.bound[:,2])
        objs = problem.objfunc(vars)
        if (isnan(problem.limit[1,1]))
            s.objectives = objs
        else
            s.objectives = norm(objs, problem.limit[:,1], problem.limit[:,2])
        end
    end
end

これで最適化問題が定義できるようになりました．先程のrosenbrock()を使って，以下のように最適化問題を定義します．

numvars = 3
bound = [[-5.0 5.0];[-5.0 5.0];[-5.0 5.0]]
problem = DoubleProblem(rosenbrock, numvars=numvars, bound=bound)

3. 粒子群最適化アルゴリズム

PSOは以下の流れで最適化を行います．

初期化：複数の粒子を作り，粒子の位置=変数を初期化する
評価：粒子の位置=変数から目的関数値を計算する
ベスト粒子の更新：global best(粒子群のうち最も良い評価の位置, gbest)と，personal best(その粒子がこれまで移動した中で最も良い評価の位置, pbest)を目的関数値を用いて更新する
飛翔：粒子の速度と位置を，gbestとpbestを用いて更新する
突然変異：粒子位置に突然変異を加える
所定の繰り返し回数に到達するまで，2に戻る．

これを書き下したのが以下です．

function pso(problem::DoubleProblem; iter, nump=problem.numvars*10, c1=NaN, c2=NaN, w=NaN, vecr=false)
    # 粒子群の初期化
    swarm = [Particle(problem.numvars) for i = 1:nump]
    gbest = Particle(problem.numvars)
    evaluate!(problem, swarm)
    updategbest!(swarm, gbest)
    updatepbest!(swarm)

    # 最適化の実行
    while (iter > 0)
        updateswarm!(swarm, gbest, c1, c2, w, vecr)
        perturbation!(swarm)
        evaluate!(problem, swarm)
        updategbest!(swarm, gbest)
        updatepbest!(swarm)
        iter -= 1
    end
    return gbest
end
# 粒子の初期化
function Particle(numvars, numobjs=1, numcons=0)
    objectives = ones(numobjs) .* Inf
    variables  = rand(numvars)
    velocities = zeros(numvars)
    pbestvars  = copy(variables)
    pbestobjs = ones(numobjs) .* Inf
    constraints = zeros(numcons)
    return Particle(objectives, variables, velocities, pbestvars, pbestobjs, constraints)
end
# 粒子速度・位置の更新
function updateswarm!(swarm::Vector{Particle{Float64}}, gbest::Particle, c1::Float64=NaN, c2::Float64=NaN, w::Float64=NaN, vecr::Bool=false)
    for p in swarm
        if(vecr)
            r1 = rand(length(p.variables))
            r2 = rand(length(p.variables))
        else
            r1 = rand()
            r2 = rand()
        end
        if(isnan(c1)) c1 = rand() * 0.5 + 1.5 end
        if(isnan(c2)) c2 = rand() * 0.5 + 1.5 end
        if(isnan(w))  w = rand() * 0.4 + 0.1 end
        
        # 位置と速度の更新
        p.velocities = w * p.velocities + c1 * r1 * (p.pbestvars - p.variables) + c2 * r2 * (gbest.variables - p.variables)
        p.variables += p.velocities

        # 境界チェック
        for i = 1:length(p.variables)
            if (p.variables[i] > 1.0)
                p.variables[i] = 1.0
                p.velocities[i] *= -1.0
            elseif (p.variables[i] < 0.0)
                p.variables[i] = 0.0
                p.velocities[i] *= -1.0
            end
        end
    end
end

# 突然変異(未実装)
function perturbation!(swarm::Vector{Particle{Float64}})

end

# pbestの更新
function updatepbest!(swarm::Vector{Particle{Float64}})
    for p in swarm
        if (p.objectives < p.pbestobjs)
            p.pbestvars = copy(p.variables)
            p.pbestobjs = copy(p.objectives)
        end
    end
end

# gbestの更新
function updategbest!(swarm::Vector{Particle{Float64}}, gbest::Particle{Float64})
    for p in swarm
        if (p.objectives < gbest.objectives)
            gbest.variables = copy(p.variables)
            gbest.objectives = copy(p.objectives)
        end
    end
end

ここで，粒子は，他の問題と共通の解の抽象型(AbstractSolution)の具体型としています．粒子は，上述のように位置・速度・評価値と，pbestの位置・評価値を持つ必要があるので，これらを型のメンバとしています．

# 抽象型
abstract type AbstractSolution{T <: Number} end
# 具体型
mutable struct Particle{T} <: AbstractSolution{T}
    objectives::Vector{Float64} # 評価値
    variables::Vector{T} # 位置
    velocities::Vector{T} # 速度
    pbestvars::Vector{T} # pbestの位置
    pbestobjs::Vector{Float64} # pbestの評価値
    constraints::Vector{T} # 制約違反量，未使用．
end

PSOの実行の際には，以下のパラメータを必要とします．
* 粒子の数: nump * 繰り返し回数: iter

また，粒子の位置・速度の更新式は次式の様になっているのですが，ここでもパラメータを指定する必要があります．

$v_{t+1} = wv_t +c_1 r_1(x_{pbest}-x_t)+c_2 r_2 (x_{gbest}-x_t)$
$x_{t+1} = x_t + v_{t+1}$

重みw
重み(加速係数)c1, c2

さらに，乱数r1, r2を，1つの解で1つの値とするか，変数ごとに異なる値とする(ベクトル化する)かでも探索の傾向が異なってくるので，これもパラメータとします．

r1, r2をベクトル化するフラグ: vecr

これらのパラメータは，指定する必要がなければデフォルトのパラメータで動くようにします．粒子の数numpは，変数の数×10としています．また，重みw, c1, c2は所定の範囲で乱数としています．今回は，OMOPSOで用いられる乱数範囲と同一としました．

4. 最適化の結果

以下のpso()関数を実行して，最適化を行います．

@time res = pso(problem, iter=1000, nump=100)

何度か実行したところ，以下のような結果となりました．何回かに1回，近似解で止まってしまいますが，ほとんどの試行で最適解(評価値0)に収束しています．
また，Juliaは初回実行時はプリコンパイルがあるので実行時間が長いという話を聞きますが，初回は0.3秒，その後は0.1秒ほどと，その傾向がそのまま出ている気がします．

objectives: [0.0]
variables: [1.0, 1.0, 1.0]
  0.260429 seconds (1.75 M allocations: 165.164 MiB, 14.49% gc time)

objectives: [0.0]
variables: [1.0, 1.0, 1.0]
  0.164344 seconds (1.65 M allocations: 159.511 MiB, 35.80% gc time)

objectives: [2.1172833042747135]
variables: [0.04515555334956112, 0.0031784243458048422, -0.04600862315194654]
  0.088791 seconds (1.83 M allocations: 177.379 MiB, 22.27% gc time)

5. 課題

PSOの実装とRosenbrock関数の最適化ができることを確認しました．今回，アルゴリズムとしては以下ができていないので，今後，少しずつ取り組みたいと思っています．

多目的および制約付き最適化への拡張
冒頭に述べたように，拡張できるような配慮をしましたので，多目的・制約付き最適化でも適用できるようにしたいです．
突然変異の実装
オリジナルのPSOには突然変異はありませんが，何らかの手法を使えるようにしたいです．
w, c1, c2のパラメータ
w, c1, c2の与え方には様々な検討がなされています．別途指定できるようにしたいです．

また，Juliaの学習という意味では以下も行いたいと思います．

Juliaのデザインパターンの学習
今回のような設計が良かったかどうかの感覚があまりつかめていません．他の人のコードを見たりデザインパターンを調べて，実装もしながら慣れていきたいと思います．
探索経過を確認できるような可視化・グラフ描画
グラフ描画の練習も兼ねて，探索経過を可視化する機能を作りたいです．
モジュール分割
プログラム自体がやや長くなってしまったので，問題とアルゴリズムのような形でモジュール化し分割したいと思います．