【Julia入門】modelicaとJuliaで感度分析 - システムとモデリング

前回に続いてopen modelicaのJulia APIOMJuliaを使用して計算モデルを解析していきます。前回の記事は↓こちら

otepipi.hatenablog.com

パラメーターの感度分析

パラメーターの感度分析とは、あるパラメーターの変化がシステムのアウトプット変数にどれだけ影響を与えるかを分析する手法です。仮にアウトプット変数を $y(t)$ 、パラメーターを $a$ とすると、Open Modelicaは感度として $\frac{\partial {y(t)}}{\partial a}$ を計算します。この機能を実際に使用してみます。

計算モデル自体は前回と同様のものを使用します。

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

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での感度解析コマンドはsensitivity()ですので、その命令を使用するJuliaコードを書きます。

using Plots
using OMJulia

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

tnk.setInputs("md_i=3");

a = ["A"];
y = ["h"];

r1,r2 = tnk.sensitivity(a,y);

plot(tm,r2,title="sensitivity analysis",xlabel="time [s]",ylabel="sensitivity");

今回、パラメーターをA(タンクの断面積）、出力変数をh（水位）として、シミュレーション時間(10秒)時々の $\frac{\partial {h}}{\partial A}$ の値を計算してプロットします。