JuliaパッケージModelingToolkit.jlでタンクモデルを作成する

久しぶりの投稿です。
今回は非因果モデリングのためのJuliaパッケージModelingToolkit.jlを使用してタンクモデルを作成してみます。

ModelingToolkit.jlについて

ModelingToolkitはJuliaネイティブのモデリングツールで、方程式ベースでモデルを記述できます。方程式ベースのため入力と出力を区別しない非因果的なモデリングが可能です。その他の非因果モデリングツールはmodelicaや同じくJuliaベースのmodia.jlが有名ですが、ModelingToolkitとそれらの相違点は以下の記事に詳しく書かれています。

https://www.stochasticlifestyle.com/modelingtoolkit-modelica-and-modia-the-composable-modeling-future-in-julia/

また、ModelingToolkitのドキュメントは以下になります。

https://mtk.sciml.ai/stable/

今回の題材

今回の題材は以下のようなものになります。

f:id:Otepipi:20210724145919p:plain

ポンプから一定流量qinでタンクに水が流れ込みます。タンクの下部には開度Rのバルブが取り付けられており、水位hに応じて水が流出します。
これをModelingToolkitでモデリングしてみます。

モデリングの考え

以下のようにモデリングしてみます。

f:id:Otepipi:20210724150016p:plain

なお、タンクのこれらの関係式は以下のようになります。

f:id:Otepipi:20210724150041p:plain

$C \frac{dh}{dt} = q_{in} - q_{out} \\ q_{ou$

実装

では、この簡単なモデルを実装してみましょう。

使用パッケージの宣言

以下のように使用するパッケージを宣言します。

using ModelingToolkit
using DifferentialEquations, Plots

モデリングに必要なModelingToolkitと、微分方程式を解くのに必要なDifferentialEquations、そして結果の描画に必要なPlotsですね。

タンクのモデリング

タンクは以下のようにモデリングできます。

function tank(;name)
  @parameters C R
  @variables t h(t) qin(t) qout(t)
  D = Differential(t)
  eqs = [ D(h) ~ (qin - qout)/C,
        qout ~ R * h]
  ODESystem(eqs,t,[h,qin,qout],[C,R];name=name)
end

@parametersでパラメーターを宣言しています。ここでパラメーターとはシミュレーション中に値が変化しないもの（前後では変化する）としておきます（明確な定義はありません。しかしmodelicaはそのようなイメージでパラメーターを扱っています)。また、@variablesは変数で、シミュレーション中に値が変化する可能性があるものとします。

D = Differential(t)は時間tの微分を定義しています。続くeqsで方程式を定義します。qout ~ R * hのようにModelingToolkitでは等号は~で表現します。

ここでいくつか注意事項です。

ModelingToolkitではパラメーターと変数の違いが不明確で、パラメーターもシミュレーション中に値が変化しても実は問題ありません。ただし、パラメーターの経時変化はシミュレーション結果に記録されません。

D(h) ~ (qin - qout)/Cのように、微分項は必ず数式の左辺に単独で配置する必要があります。たとえばこれをC*D(h) ~ (qin - qout)とする場合、微分項に係数があるためエラーになります。

ポンプのモデリング

続いてポンプのモデリングは以下のようになります。

function pump(;name)
   @variables t q(t)
   eqs = q ~ 1
   ODESystem(eqs,t,[q],[];name=name)
end

注意点はタンクのモデリングの際と特に違いありません。

システムのモデリングと簡略化

これまではタンクとポンプという各コンポーネントのモデリングでしたが、次にこれらを繋ぎ合わせてシステムとしてモデリングしていきます。

モデルの呼び出し

タンクとポンプはまだ関数を作成しただけですので、それらを呼び出してみます。

@named tankA = tank()
@named pumpA = pump()

それぞれtankA、pumpAと名付けました。

モデルの接続

タンクとポンプを接続します。モデル図を見てわかりますように、ポンプの変数qとタンクの変数qinは等しいです。

f:id:Otepipi:20210724150016p:plain

eqs = tankA.qin ~ pumpA.q

@named connected = compose(ODESystem(eqs),tankA,pumpA)

これはモデルconnectedがtankAとpumpAをサブシステムとして内包していると解釈できます。

モデルの簡略化

ここでモデルconnectedの数式を見てみましょう。

equations(connected)

で確認できます。

4-element Vector{Equation}:
 tankA₊qin(t) ~ pumpA₊q(t)
 Differential(t)(tankA₊h(t)) ~ (tankA₊C^-1)*(tankA₊qin(t) - tankA₊qout(t))
 tankA₊qout(t) ~ tankA₊R*tankA₊h(t)
 pumpA₊q(t) ~ 1

この式にはtankA₊qin(t)とpumpA₊q(t)が両方存在していますが、両者には等号の関係があるのでどちらかは削除できるはずです。他にもpumpA.q(t)やtankA.qoutも削除できそうです。このように変数を消去して簡略化する命令がstructural_simplifyになります。

simplified_sys = structural_simplify(connected)

ちゃんと簡略化されたか確認してみましょう。

equations(simplified_sys)

1-element Vector{Equation}:
 Differential(t)(tankA₊h(t)) ~ (tankA₊C^-1)*(1.0 - (tankA₊R*tankA₊h(t)))

4つあった数式が1つにまで簡略化されました。

シミュレーション

最後に、このモデルを微分方程式ソルバーで解いていきます。

x0 = [tankA.h => 1] ##初期値
p = [tankA.C => 1,
     tankA.R => 1/5] ##パラメーターの値

prob = ODEProblem(simplified_sys, x0,(0.0, 20.0),p)
sol = solve(prob,Tsit5())

plot(sol) ## 結果の描画

f:id:Otepipi:20210724150202p:plain

無事、結果の描画までたどり着くことができました。

さいごに、今回のコードを置いておきます。

using ModelingToolkit
using DifferentialEquations, Plots

function tank(;name)
  @parameters C R
  @variables t h(t) qin(t) qout(t)
  D = Differential(t)
  eqs = [ D(h) ~ (qin - qout)/C,
        qout ~ R * h]
  ODESystem(eqs,t,[h,qin,qout],[C,R];name=name)
end

function pump(;name)
    @variables t q(t)
    eqs = q ~ 1
    ODESystem(eqs,t,[q],[];name=name)
end

@named tankA = tank()
@named pumpA = pump()

eqs = tankA.qin ~ pumpA.q

@named connected = compose(ODESystem(eqs),tankA,pumpA)


simplified_sys = structural_simplify(connected)
x0 = [tankA.h => 1] ##初期値
p = [tankA.C => 1,
     tankA.R => 1/5] ##パラメーターの値

prob = ODEProblem(simplified_sys, x0,(0.0, 20.0),p)
sol = solve(prob,Tsit5())

plot(sol) ## 結果の描画