Notebook

4つのタンクに液体を分配するシステムのパラメータの推定¶

JModelica.org Users Gauide - Ver.2.10 6.5.2.4. Parameter estimation の例題を実習します。

この例題のモデルは、JModelica の配布パッケージの中の QuadTankPack.mopに含まれています。また、以下のように入力することによってデモを実行することができます。

%matplotlib inline
from pyjmi.examples import qt_par_est_casadi
qt_par_est_casadi.run_demo()

ここでは、QuadTankPack.mop からシステムモデルとパラメータ推定用モデルを抜粋しています。

本文書は、Users Guide の例題を解説するものであり、モデルのソースコードやオペレーションで使用する Python スクリプトの著作権は Modelon AB にあります。

4つのタンクに液体を分配するシステム (quadruple tank system) の概要¶

このシステムは、ポンプ２個、３方型バルブ２個、タンク４個から構成されています。

ポンプ１とポンプ２の吐出流量は、電圧によって制御します。
ポンプ１から供給される液体の流路は３方型バルブ１で分岐し、タンク１とタンク４に至ります。
ポンプ２から供給される液体の流路は３方型バルブ２で分岐し、タンク２とタンク３に至ります。
各タンクには底部に穴があり、タンク内の液体が流出します。
ダンク１の上にタンク３があり、タンク3から流出した液体はタンク１に流れ込みます。
タンク２の上にタンク４があり、タンク４から流出した液体はタンク２に流れ込みます。

状態変数¶

$x_1, x_2, x_3, x_4$ : タンク1〜タンク4の液位 [m]

初期条件¶

$x_1(0), x_2(0), x_3(0), x_4(0)$ : タンク１〜タンク4の液位の初期値 [m]

制御変数 (入力変数)¶

$u_1, u_2$ : ポンプ1、ポンプ２の電圧 [V]

方程式¶

タンク1〜タンク4の流量バランスから以下の方程式が得られます。方程式の導出は最後の方に記述します。

$\begin{align} \dot{x_1} &= -\frac{a_1}{A_1} \sqrt{2 g x_1} + \frac{a_3}{A_1} \sqrt{2 g x_3} + \frac{\gamma_1 k_1}{A_1} u_1\\ \dot{x_2} &= -\frac{a_2}{A_2} \sqrt{2 g x_2} + \frac{a_4}{A_2} \sqrt{2 g x_4} + \frac{\gamma_2 k_2}{A_2} u_2\\ \dot{x_3} &= -\frac{a_3}{A_3} \sqrt{2 g x_3} + \frac{(1-\gamma_2) k_2}{A_3} u_2 \\ \dot{x_4} &= -\frac{a_4}{A_4} \sqrt{2 g x_4} + \frac{(1-\gamma_1) k_1}{A_4} u_1 \end{align}$

パラメータ¶

$k_1, k_2$ : ポンプの流量と電圧の関係を表す係数 [m3/s.V]
$\gamma_1, \gamma_2$ : ３方型バルブの流量比
$A_1, A_2, A_3, A_4$ : タンク1〜タンク4の断面積 [m2]
$a_1, a_2, a_3, a_4$ : 圧力損失を考慮して補正したタンク出口の断面積 [m2]
$g$ ：　重力加速度

パラメータ推定問題の概要¶

ポンプの電圧 $u_1, u_2$ を変化させてタンクの液位 $x_1, x_2, x_3, x_4$ の時刻歴を測定したデータを使用して、上記システムのパラメータの一部を推定します。

次のような手順を行います。

測定データを読み込みます。
測定データの初期条件 $x_1(0), x_2(0), x_3(0), x_4(0)$ と制御変数 $u_1(t), u_2(t)$ を使用して、タンクシステムモデルのシミュレーションを行います。
結果を比較し、推定するパラメータを決めます。
パラメータ推定用のモデルを作成します。
パラメータ推定計算を行います。
測定データと推定したパラメータを使用したシミュレーション結果を比較します。

準備¶

以下が必要です。

測定データ
Modelica で作成した上記のタンクシステムのモデル

MATLAB の mat ファイル形式の測定データが JModelica.org に含まれています。

Windows の場合は

%copy c:\JModelica.org-2.14\install\Python_64\pyjmi\examples\files\qt_par_est_data.mat .

でコピーします。

Mac で Docker にインストールした JModelica.org を使ってる場合は

%cp /usr/local/jmodelica/Python/pyjmi/examples/files/qt_par_est_data.mat .

でコピーします。

In [1]:

%cp /usr/local/jmodelica/Python/pyjmi/examples/files/qt_par_est_data.mat .

タンクシステムのモデルは QuadTankPack.mop から以下の部分を抜粋しました。

In [2]:

%%writefile QuadTankSystem.mo
package QuadTankSystem

  model QuadTank
    
    // Parameters
    parameter Modelica.SIunits.Area A1 = 4.9e-4, A2 = 4.9e-4, A3 = 4.9e-4, A4 = 4.9e-4;
    parameter Modelica.SIunits.Area a1 = 0.03e-4, a2 = 0.03e-4, a3 = 0.03e-4, a4 = 0.03e-4;
    parameter Modelica.SIunits.Acceleration g = 9.81;
    parameter Real k1_nmp(unit="m3^/s/V") = 0.56e-6, k2_nmp(unit="m^3/s/V") = 0.56e-6;
    parameter Real g1_nmp = 0.30, g2_nmp = 0.30;
    
     // Initial tank levels
    parameter Modelica.SIunits.Length x1_0 = 0.04102638;
    parameter Modelica.SIunits.Length x2_0 = 0.06607553;
    parameter Modelica.SIunits.Length x3_0 = 0.00393984;
    parameter Modelica.SIunits.Length x4_0 = 0.00556818;

    // Tank levels
    Modelica.SIunits.Length x1(start=x1_0);
    Modelica.SIunits.Length x2(start=x2_0);
    Modelica.SIunits.Length x3(start=x3_0);
    Modelica.SIunits.Length x4(start=x4_0);
  
    // Inputs
    connector InputVoltage = input Modelica.SIunits.Voltage;
    InputVoltage u1;
    InputVoltage u2;

  equation
    der(x1) = -a1/A1*sqrt(2*g*x1) + a3/A1*sqrt(2*g*x3) + g1_nmp*k1_nmp/A1*u1;
    der(x2) = -a2/A2*sqrt(2*g*x2) + a4/A2*sqrt(2*g*x4) + g2_nmp*k2_nmp/A2*u2;
    der(x3) = -a3/A3*sqrt(2*g*x3) + (1-g2_nmp)*k2_nmp/A3*u2;
    der(x4) = -a4/A4*sqrt(2*g*x4) + (1-g1_nmp)*k1_nmp/A4*u1;
  end QuadTank;
end QuadTankSystem;

Overwriting QuadTankSystem.mo

測定データを読み込む¶

ファイルから測定データを読み込みます。

０．０１秒刻みで１２０秒分のデータを間引いて、６０　秒から　120　秒まで１秒毎のデータを抽出します。
タンク液位の測定値は cm 単位なので 100 で割って m 単位に変換します。

In [3]:

from scipy.io.matlab.mio import loadmat
# Load measurement data from file
data = loadmat("qt_par_est_data.mat", appendmat=False)
# Extract data series
t_meas = data['t'][6000::100, 0] - 60
y1_meas = data['y1_f'][6000::100, 0] / 100
y2_meas = data['y2_f'][6000::100, 0] / 100
y3_meas = data['y3_d'][6000::100, 0] / 100
y4_meas = data['y4_d'][6000::100, 0] / 100
u1 = data['u1_d'][6000::100, 0]
u2 = data['u2_d'][6000::100, 0]

パラメータ修正前のモデルのシミュレーションを行う¶

モデルをコンパイルして FMU を生成し、ロードします。

In [4]:

from pymodelica import compile_fmu
from pyfmi import load_fmu
# Compile and load FMU, which is used for simulation
model = load_fmu(compile_fmu('QuadTankSystem.QuadTank', "QuadTankSystem.mo"))

測定データの初期条件をモデルに設定します。

In [5]:

import numpy as N
# Set initial states in model, which are stored in the optimization problem
x_0_names = ['x1_0', 'x2_0', 'x3_0', 'x4_0']
x_0_values = [y1_meas[0], y2_meas[0], y3_meas[0], y4_meas[0]]
model.set(x_0_names, x_0_values)

制御変数のオブジェクトを作成し、入力変数として設定してモデルのシミュレーションを実行します。

In [6]:

u = N.transpose(N.vstack([t_meas, u1, u2]))
# Simulate model response with nominal parameter values
res_sim = model.simulate(input=(['u1', 'u2'], u), start_time=0., final_time=60.)

Final Run Statistics: --- 

 Number of steps                                 : 183
 Number of function evaluations                  : 309
 Number of Jacobian evaluations                  : 4
 Number of function eval. due to Jacobian eval.  : 16
 Number of error test failures                   : 36
 Number of nonlinear iterations                  : 305
 Number of nonlinear convergence failures        : 0

Solver options:

 Solver                   : CVode
 Linear multistep method  : BDF
 Nonlinear solver         : Newton
 Linear solver type       : DENSE
 Maximal order            : 5
 Tolerances (absolute)    : 1e-06
 Tolerances (relative)    : 0.0001

Simulation interval    : 0.0 - 60.0 seconds.
Elapsed simulation time: 0.0289959907532 seconds.

シミュレーション結果を取得します。

In [7]:

# Load simulation result
x1_sim = res_sim['x1']
x2_sim = res_sim['x2']
x3_sim = res_sim['x3']
x4_sim = res_sim['x4']
t_sim = res_sim['time']
u1_sim = res_sim['u1']
u2_sim = res_sim['u2']

測定データと修正前のモデルのシミュレーション結果を比較する¶

タンクの液位とポンプの電圧をプロットします。

In [8]:

%matplotlib notebook
from matplotlib import pyplot as plt
# Plot measurements and inputs
plt.figure(1)
plt.subplot(2, 2, 1)
plt.plot(t_meas, y3_meas)
plt.plot(t_sim, x3_sim)
plt.title('x3')
plt.grid()
plt.subplot(2, 2, 2)
plt.plot(t_meas, y4_meas)
plt.plot(t_sim, x4_sim)
plt.title('x4')
plt.grid()
plt.subplot(2, 2, 3)
plt.plot(t_meas, y1_meas)
plt.plot(t_sim, x1_sim)
plt.title('x1')
plt.xlabel('t[s]')
plt.grid()
plt.subplot(2, 2, 4)
plt.plot(t_meas, y2_meas)
plt.plot(t_sim, x2_sim)
plt.title('x2')
plt.xlabel('t[s]')
plt.grid()

plt.figure(2)
plt.subplot(2, 1, 1)
plt.plot(t_meas, u1)
plt.plot(t_sim, u1_sim, 'r',linestyle='dashed')
plt.ylabel('u1')
plt.grid()
plt.subplot(2, 1, 2)
plt.plot(t_meas, u2)
plt.plot(t_sim, u2_sim, 'r',linestyle='dashed')
plt.ylabel('u2')
plt.xlabel('t[s]')
plt.grid()

パラメータ推定用のモデルを作成する。¶

タンクの液位 x1 と x2 は、 x3 や x4 と比べると、測定値とシミュレーションの結果のずれが大きくなっています。そこで x1 と x2 に与える影響が大きいパラメータ a1 と a2 の推定を行うことにします。

まず、測定データの液位の初期値を確認します。

In [9]:

print(x_0_values)

[0.062715498260467531, 0.060437812711968755, 0.023949266862170089, 0.023249081133919845]

タンクシステムのモデルを継承して、次のようなパラメータ推定用のモデルを作成します。このモデルも　QuadTankPack.mop からの抜粋です。

初期値 x1_0 〜 x4_0 は、上で確認した測定データの初期値を設定します。
推定を行うパラメータ a1と a2 は、free=true を設定します。

In [10]:

%%writefile QuadTankOpt.mop
optimization QuadTank_ParEstCasADi(startTime=0, finalTime=60)
    extends QuadTankSystem.QuadTank(
        x1(fixed=true), x1_0=0.062715498260467531,
        x2(fixed=true), x2_0=0.060437812711968755,
        x3(fixed=true), x3_0=0.023949266862170089,
        x4(fixed=true), x4_0=0.023249081133919845,
        a1(free=true, min=0, max=0.1e-4),
        a2(free=true, min=0, max=0.1e-4));
end QuadTank_ParEstCasADi;

Overwriting QuadTankOpt.mop

パラメータ推定用モデルをロードします。

In [11]:

from pyjmi import transfer_optimization_problem
op = transfer_optimization_problem("QuadTank_ParEstCasADi",[ "QuadTankOpt.mop","QuadTankSystem.mo"])

パラメータ推定用モデルは、他の最適化問題のように目的関数 (objective や objectiveIntegrand など) を設定しません。その代わり ExternalData オブジェクトを使用します。このオブジェクトは、x1 と x2 について測定値とモデルの値の差の自乗の積分という目的関数を生成します。２つの入力変数 u1 と u2 についてはペナルティ法を導入します。入力変数は ExternalData の eliminated パラメータで消去することも可能です。

data_x1, data_x2, data_u1, data_u2 は、測定データです。
quad_pen には、ペナルティ法を適用する変数と測定データの組み合わせを設定します。
Q は、他の要素に対する変数の重みを表す行列で、対角行列で制御変数 u1, u2 の重みが大きくなるように設定しています。

In [12]:

from collections import OrderedDict
from pyjmi.optimization.casadi_collocation import ExternalData
Q = N.diag([1., 1., 10., 10.])
data_x1 = N.vstack([t_meas, y1_meas])
data_x2 = N.vstack([t_meas, y2_meas])
data_u1 = N.vstack([t_meas, u1])
data_u2 = N.vstack([t_meas, u2])
quad_pen = OrderedDict()
quad_pen['x1'] = data_x1
quad_pen['x2'] = data_x2
quad_pen['u1'] = data_u1
quad_pen['u2'] = data_u2
external_data = ExternalData(Q=Q, quad_pen=quad_pen)

パラメータ推定計算を行う¶

オプションで　ExternalData オブジェクトを設定してパラメータ推定用問題を解きます。 initial trajectory と nominal trajectory に先程行ったパラメータ修正前のシミュレーションの結果ファイルを設定しています。

In [13]:

# Set optimization options and optimize
opts = op.optimize_options()
opts['n_e'] = 60 # Number of collocation elements
opts['external_data'] = external_data
opts['init_traj'] = res_sim
opts['nominal_traj'] = res_sim
res = op.optimize(options=opts) # Solve estimation problem

Warning: Could not find initial trajectory for variable temp_1. Using initialGuess attribute value instead.
Warning: Could not find initial trajectory for variable temp_2. Using initialGuess attribute value instead.
Warning: Could not find nominal trajectory for variable temp_1. Using nominal attribute value instead.
Warning: Could not find nominal trajectory for variable temp_2. Using nominal attribute value instead.
Warning: Could not find initial trajectory for variable temp_1. Using initialGuess attribute value instead.
Warning: Could not find initial trajectory for variable temp_2. Using initialGuess attribute value instead.

Initialization time: 0.45 seconds

Total time: 0.52 seconds
Pre-processing time: 0.00 seconds
Solution time: 0.05 seconds
Post-processing time: 0.02 seconds

変数 temp_1 と temp_2 の　initial trajectory や nominal trajectory が見つからないという Warning が出ていますが、このような名前の変数は使用していません。詳細は不明です。この Warning は、オリジナルのデモ qt_par_est_casadi でも発生します。

測定結果とパラメータ修正を行ったシミュレーション結果を比較する¶

推定計算結果のパラメータを表示します。

In [14]:

# Extract estimated values of parameters
a1_opt = res.initial("a1")
a2_opt = res.initial("a2")
# Print estimated parameter values
print('a1: ' + str(a1_opt*1e4) + ' cm^2')
print('a2: ' + str(a2_opt*1e4) + ' cm^2')

a1: 0.0265770557815 cm^2
a2: 0.0271381846631 cm^2

推定計算の結果から、液位とポンプ電圧を抽出します。

In [15]:

# Load state profiles
x1_opt = res["x1"]
x2_opt = res["x2"]
x3_opt = res["x3"]
x4_opt = res["x4"]
u1_opt = res["u1"]
u2_opt = res["u2"]
t_opt = res["time"]

プロットして測定データと比較します。

In [16]:

# Plot measurements and inputs
plt.figure(3)
plt.subplot(2, 2, 1)
plt.plot(t_meas, y3_meas)
plt.plot(t_opt, x3_opt)
plt.title('x3')
plt.grid()
plt.subplot(2, 2, 2)
plt.plot(t_meas, y4_meas)
plt.plot(t_opt, x4_opt)
plt.title('x4')
plt.grid()
plt.subplot(2, 2, 3)
plt.plot(t_meas, y1_meas)
plt.plot(t_opt, x1_opt)
plt.title('x1')
plt.xlabel('t[s]')
plt.grid()
plt.subplot(2, 2, 4)
plt.plot(t_meas, y2_meas)
plt.plot(t_opt, x2_opt)
plt.title('x2')
plt.xlabel('t[s]')
plt.grid()

plt.figure(4)
plt.subplot(2, 1, 1)
plt.plot(t_meas, u1)
plt.plot(t_opt, u1_opt, 'r',linestyle='dashed')
plt.ylabel('u1')
plt.grid()
plt.subplot(2, 1, 2)
plt.plot(t_meas, u2)
plt.plot(t_opt, u2_opt, 'r',linestyle='dashed')
plt.ylabel('u2')
plt.xlabel('t[s]')
plt.grid()

検証のため、次のようにタンクシステムモデルに推定したパラメータを設定してシミュレーションを行ってみましたが全く同じ結果が得られました。

model = load_fmu('QuadTankSystem_QuadTank.fmu')
model.set('a1', a1_opt)
model.set('a2', a2_opt)
model.set(x_0_names, x_0_values)
res_optsim = model.simulate(input=(['u1', 'u2'], u), start_time=0., final_time=60.)

付録. タンクシステムモデルの方程式の導出¶

ポンプから各タンクに供給される流量¶

タンク1〜タンク4にポンプから供給される流量は次にようになります。

$\begin{align} Q_{p1} &= \gamma_1 k_1 u_1\\ Q_{p2} &= \gamma_2 k_2 u_2 \\ Q_{p3} &= (1-\gamma_2) k_2 u_2 \\ Q_{p4} &= (1-\gamma_1) k_1 u_1\\ \end{align}$

$u_1, u_2$ ポンプの電圧 [V] (制御変数)
$k_1, k_2$ ポンプの流量と電圧の関係を表す係数 [m3/s.V]
$\gamma_1, \gamma_2$ ３方型バルブの流量比

タンク出口の流量¶

タンクの出口において、内部と外部の圧力差を $\Delta p$ [Pa]、圧力損失係数を $\zeta$ 、密度を $\rho$ [kg/m3]、流速を $u$ [m/s]とすると、圧力と流速の関係は次のように表せます。

$\Delta p = \zeta \frac{\rho u^2}{2}$ 圧力差として

$\Delta p = \rho g x$ [Pa] を用い、流速

$u$ [m/s] の代わりに体積流量

$Q = a_a u$ [m3/s] を用いて整理すると、

$Q = \frac{a_a}{\sqrt{\zeta}} \sqrt{2 g x}$ となります。

$a_a$ [m2] は実際の断面積です。圧力損失を考慮して補正した断面積

$a \equiv a_a / \sqrt{\zeta}$ [m2] を導入すると、

$Q = a \sqrt{2 g x}$ となります。出口が理想的な急拡大の場合

$\zeta \sim 1$ 、

$a \sim a_a$ [m2] となります。

各タンクの出口の体積流量は、

$Q_i = a_i \sqrt{2 g x_i}, i=1,2,3,4$ となります。

タンクの流量バランス¶

液位の変化と、タンクに出入りする流量の関係は次のようになります。

$\begin{align} A_1 \dot{x_1} &= -Q_1 + Q_3 + Q_{p1} \\ A_2 \dot{x_2} &= -Q_2 + Q_4 + Q_{p2} \\ A_3 \dot{x_3} &= -Q_3 + Q_{p3} \\ A_4 \dot{x_4} &= -Q_4 + Q_{p4} \end{align}$

$Q_{pi}$ や $Q_i$ に上に述べた式を代入して整理すると、

が得られます。