【Julia入門】modelicaとJuliaで水タンクをモデリングする

※今回の内容は以下論文を再現しただけになります

https://www.researchgate.net/publication/330809787_OMJulia_An_OpenModelica_API_for_Julia-Modelica_Interactionwww.researchgate.net

前回に引き続きmodelicaとJuliaを繋ぐOMJulia関連の記事になります。

github.com

関連するOMJulia記事は以下です。

otepipi.hatenablog.com

今回のモデリング対象

今回は下図のような水の入ったタンクをモデリングします。タンク上部から水が流入し、タンク底からは重力に従い水が流出します。

f:id:Otepipi:20210202192231p:plain — 水タンク概略図

記号の意味を以下にまとめます。

記号	意味	Head3
$$\dot{m_i}$$	水の流入速度	kg/s
m	タンク内の水量	kg
$$\dot{m_e}$$	水の流出速度	kg/s
$$V$$	タンク内の水体積	$dm³$
$$\rho$$	水の密度	kg/L
$$h^{\zeta}$$	代表高さ	dm
$$h$$	水位	dm

モデルの数式

　上記モデルは以下のように数式で表すことができます。

f:id:Otepipi:20210202193717p:plain

modelicaによるモデリング

水タンクモデルをmodelicaで書き下した場合以下のようになります。今回は水の流入量md_iを入力に、タンク水位hを出力に指定しています。この段階では入力md_iの値を指定していないのがポイントです。

package WaterTank
// Package for simulating


    model ModWaterTank
        constant Real rho = 1 "Density";
        parameter Real A = 5, K=5, h_s=3;
        parameter Real h_0=0.5, m_0=rho*h_0*A;
        Real m(start=m_0,fixed=true) "Mass in tank, kg";
        Real V, md_e;
        input Real md_i;
        output Real h;
    equation
        der(m) = md_i - md_e;
        m = rho*V;
        V = A*h;
        md_e = K*sqrt(h/h_s);
    end ModWaterTank;

//End package
end WaterTank;

OMJuliaによるシミュレーション

作成したmodelicaモデルをシミュレーションします。

using Plots
using OMJulia
using ControlSystems

tnk = OMJulia.OMCSession();
tnk.ModelicaSystem("WaterTank.mo","WaterTank.ModWaterTank");

tnk.setInputs("md_i=3");
tnk.setSimulationOptions(["stopTime=1e4","stepSize=10"]);
tnk.simulate();
tm,h = tnk.getSolutions(["time","h"]);
plot(tm,h)
plot!(title="Water tank level")
plot!(xlabel="time [s]")
plot!(ylabel="h [dm]")

今回tnk.setInputs("md_i=3");とすることで流入速度を3一定としました。

f:id:Otepipi:20210202195129p:plain

シミュレーション結果によるとモデルはすぐに定常状態h=1.08になります。

次に、パラメーターh_0の値を先程の定常状態の値に置き換えてみます。加えて、md_iの値を一定ではなく時系列で変動するようにしてみます。

tnk.setParameters("h_0=$(h[end])");
tnk.setSimulationOptions(["stopTime=10","stepSize=0.02"]);
tnk.setInputs("md_i = [(0,3),(2,3),(2,4),(6,4),(6,2),(10,2)]");
tnk.simulate();
tm,h = tnk.getSolutions(["time","h"]);
plot(tm,h,title="Water tank level",xlabel="time [s]",ylabel="h [dm]")

結果は以下のようになります。

f:id:Otepipi:20210202201738p:plain

今回はtnk.setInputs("md_i = [(0,3),(2,3),(2,4),(6,4),(6,2),(10,2)]");で流入量を時系列で変動させています。これをグラフにすると以下のようになります。

f:id:Otepipi:20210202201907p:plain

モデルの線形化

モデルの線形化により以下のような状態空間モデルを構築することができます。$x$は状態変数、$y$が出力変数、$u$が入力変数です。$A,B,C,D$は状態空間行列です。

$$\dot{x} = Ax + Bu $$ $$y = Cx + Du $$

OMJuliaではlinearise()でmodelicaの計算モデルを線形化することができます。 ※2021年2月現在 open modelica v1.13では動作しますが最新のv1.16では動作しません。

tnk.setLinearizationOptions("stopTime=1e-6");
A,B,C,D = tnk.linearize();

このとき状態空間行列は以下のようになります。

f:id:Otepipi:20210202202804p:plain

また、状態変数・出力変数・入力変数が何を示しているかはそれぞれgetLinearStates() 、getLinearOutputs()、getLinearInputs()でわかります。

f:id:Otepipi:20210202203253p:plain

この線形化モデルをJuliaの制御系パッケージ ControlSystemsに渡すこともできます。

using ControlSystems
sys = ss(A,B,C,D)
tf(sys)

f:id:Otepipi:20210202203558p:plain

今回はここまでにします。やはりmodelicaとJuliaの組み合わせには可能性を感じます。

システムとモデリング

modelica, Julia, Design Structure Matrix, SysML, 他モデリング全般について。