Modelicaで最適化を試す

Modelicaで最適化(動的計画法)を行う場合いくつか選択肢があるが
OpenModelica, JModelicaにはOptimizationクラスが用意されており最適化を行うことが出来る

最適化の方法の概略は以下を参考
Modelicaで最適化をする方法をまとめてみた(OSS編)
https://qiita.com/UedaShigenori/items/aae444326b28ccc8b20f

JModelicaの場合

実行方法

このUser's guideを参考にする
https://jmodelica.org/downloads/UsersGuide.pdf

Jmodelica起動
"C:\JModelica.org-2.4\IPython64.bat"

JModelicaはつまりはPythonのようだ。
残念ながら私はPythonを使ったことがないのでなんとなくでは使いこなせない
ディレクトリの移動からググる

# ディレクトリの移動
import os
os.chdir("C:\hoge")

# カレントディレクトリの確認
os.getcwd()

HelloWorld.moを読み込み計算を実行する
以下の投稿を参考にさせていただいた
https://qiita.com/TatsuyaKatayama/items/a660535fc31d784d70a4

# モデルをFMUsにコンパイル
In [7]: from pymodelica import compile_fmu
In [8]: hello_fmu = compile_fmu("HelloWorld","he.mo")

# FMUsを読み込む
In [9]: from pyfmi import load_fmu
In [10]: hello_model=load_fmu(hello_fmu)

# 読み込んだFMUsを計算する
In [11]: res = hello_model.simulate(final_time=1)

スクリプトの実行方法

毎回FMUを作って読み込んで実行する、というのは骨が折れるためスクリプトとして実行するのが良い
スクリプトの実行方法は以下

exec(open(“hogehoge.py").read())

スクリプトの例
cla, moを定義することで実行できる
claはクラス名
moはファイル名

fmu = compile_fmu(cla,mo)
loadModel=load_fmu(fmu)
res = loadModel.simulate()

最適化実行

以下のモデルを実行する

optimization VDP_Opt (objectiveIntegrand = x1^2 + x2^2 + u^2,
startTime = 0,
finalTime = 20)
// The states
Real x1(start=0,fixed=true);
Real x2(start=1,fixed=true);
// The control signal
input Real u;
equation
der(x1) = (1 - x2^2) * x1 - x2 + u;
der(x2) = x1;
constraint
u<=0.75;
end VDP_Opt;

vdp = transfer_optimization_problem("VDP_Opt_Min_Time", "VDP_Opt_Min_Time.mop")
res = vdp.optimize()

計算結果の確認

from pyjmi.common.plotting import plot_gui
plot_gui.startGUI()

使いやすいが、OMEditで出力結果matファイルを読み込んでも良い
上の画像は特に何も表示していないが、ちゃんと結果を見ても最適化されているか元の式が複雑なためよく分からない。。。

分かりやすい問題作成

最適な値が自明な問題を作成する
どうやらJModelicaは目的関数が最小となるようにデフォルトで設定されているようだ
以下の答えはx=1のとき、yが最小(81)となる。
OpenModelicaでは何故かここで躓いた。
OpenModelicaでテストしたときはy=-(x-10)^2としていた。
これはOpenModelicaの最適化機能が目的関数を最大化するからだった気がする。

ss.mop


optimization ff(objectiveIntegrand = y)
  Real y;
  input Real x;
equation
  y=(x-10)^2;
constraint
  x>=-1;
  x<=1;
end ff;

以下は結果で、うまくx=1, y=81となっている。

インスタンスを使う場合

本筋ではないが、モデルをインスタンスにして使う場合が多いと思うため
方法を確認した
ちゃんと計算できた

package nn

model mm

  Real y;
  input Real x;
equation
  y=(x-10)^2;
constraint
  x>=-1;
  x<=1;

end mm;

optimization ff(objectiveIntegrand = ii.y)
mm ii;
input Real u=ii.x;
end ff;
end nn;

パラメータ変数を計算式とする

最適化はパラメータ同定に使うのでパラメータを計算式で扱えないと不便である
まずは上記が出来るか確認する

package nn

model mm

  Real y;
  input Real x;
parameter Real o=x;
equation
  y=(x-10)^2;
constraint
  x>=-1;
  x<=1;

end mm;

optimization ff(objectiveIntegrand = ii.y)
mm ii;
input Real u=ii.x;
end ff;
end nn;

計算実行メッセージ

res=llf.optimize()
・・・

RuntimeError: An exception occured while evaluating the expression:  on line 449 of file "r:/jenkins-internal/workspace/OCTTeam/P538-JM-2.4_pipeline/JModelica/ThirdParty/CasADi/CasADi/casadi/core/function/mx_function_internal.cpp"
Cannot evaluate "function("unnamed_mx_function")" since variables [MX(ii.x)] are free.

Are the parameter values calculated?

パラメータ値を計算しようとしてない？って怒られた
どうやら仕様として無理らしい
不便だが仕方ない

どうしてもパラメータ変数を計算式にして最適化したい場合、パラメータ宣言を外してモデルを作り直すなど対策を取らないといけない

Modelicaで最適化 - JModelica編