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) の概要

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

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

状態変数

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

初期条件

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

制御変数 (入力変数)

$u_1, u_2$ : ポンプ1、ポンプ2の電圧 [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$: 3方型バルブの流量比
  • $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$ の時刻歴を測定したデータを使用して、上記システムのパラメータの一部を推定します。

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

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

準備

以下が必要です。

  • 測定データ
  • 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

測定データを読み込む

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

  • 0.01秒刻みで120秒分のデータを間引いて、60 秒から 120 秒まで1秒毎のデータを抽出します。
  • タンク液位の測定値は 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 について測定値とモデルの値の差の自乗の積分という目的関数を生成します。2つの入力変数 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()