【Julia入門】JuliaでOpen Modelicaのパラメーターを最適化する

今回は久々にmodeica関連の記事になります。

JuliaでOpen Modelicaを操作するパッケージ「OMJulia」を用いてopen modelicaのモデルのパラメーターを最適化します。関連する記事は以下になります。

otepipi.hatenablog.com

今回のモデル

f:id:Otepipi:20210129192217p:plain

今回モデリングする対象は上図です。物理の授業でよくある、ボールをある角度で投擲して地面(y=0)に落下するまでのモデルになります。今回、ModelicaとJuliaを使って、地面にぶつかるまでｘ方向に最も長く飛ぶ角度θを探索します。

Open Modelicaでモデリング

まず、適当にθ=0.5としてModelicaでモデリングしてみます。初期速度v0=0.5としています。

model RocketBall

 parameter Real v0 = 5;//初期速度
 parameter Real theta = 0.5;//打ち上げ角度
 constant Real g = 9.81;//重力加速度
          Real x;
          Real y;
          Real vy(start=v0*sin(theta));//y方向速度

equation

der(y) = vy;
der(vy) = -g;
der(x) = v0*cos(theta);

algorithm
  when y < 0 then
    terminate("simulation end");
    end when;

end RocketBall;

これをOpenModelicaでシミュレーションします。

シミュレーション結果のうち、高さｙの時間変化をプロットしたものが以下のグラフです。与えられた初速によりある程度の高さに到達した後、重力に従って地面(y=0)に落下するまでの過程が描かれています。

f:id:Otepipi:20210129192927p:plain — OpenModelica y方向

また今回最大化する、横方向ｘの時間変化を見てみます。

f:id:Otepipi:20210129193259p:plain — Open Modelica x方向

角度θ=0.5の時はx方向に2.1~2.2進むことができたようです。

Juliaを使用してパラメーターの最適化

　では、x方向に最も進む角度θを探していきます。人手でOpenModelicaのパラメーターを変えながら逐一シミュレーションしていくのは非常に手間がかかりますので、この過程をJuliaで自動化します。パッケージとしてOMJuliaを使用します。

github.com

原始的な総当り方法を使用します。 Juliaのfor~end構文でパラメーターθの値を変えながらmodelicaのシミュレーションを繰り返します。 コードは以下になります。

using OMJulia
using Plots


RocketBall = OMJulia.OMCSession() 
RocketBall.ModelicaSystem("RocketBall.mo","RocketBall")

np = 100;
par = range(0,pi/2,length=np);
position = zeros(np);


for i = 1:np
    RocketBall.setSimulationOptions(["stopTime=10", "stepSize=0.01"])
    RocketBall.setParameters("theta=$(par[i])");
    RocketBall.simulate()
    position[i] = RocketBall.getSolutions("x")[end][end];
end


plot(par,position)
xlabel!("θ")
ylabel!("position")