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 にあります。
このシステムは、ポンプ2個、3方型バルブ2個、タンク4個から構成されています。
$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} $$
ポンプの電圧 $u_1, u_2$ を変化させてタンクの液位 $x_1, x_2, x_3, x_4$ の時刻歴を測定したデータを使用して、上記システムのパラメータの一部を推定します。
次のような手順を行います。
以下が必要です。
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 .
でコピーします。
%cp /usr/local/jmodelica/Python/pyjmi/examples/files/qt_par_est_data.mat .
タンクシステムのモデルは QuadTankPack.mop から以下の部分を抜粋しました。
%%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;
ファイルから測定データを読み込みます。
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 を生成し、ロードします。
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"))
測定データの初期条件をモデルに設定します。
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)
制御変数のオブジェクトを作成し、入力変数として設定してモデルのシミュレーションを実行します。
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.)
シミュレーション結果を取得します。
# 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']
タンクの液位とポンプの電圧をプロットします。
%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 の推定を行うことにします。
まず、測定データの液位の初期値を確認します。
print(x_0_values)
タンクシステムのモデルを継承して、次のようなパラメータ推定用のモデルを作成します。このモデルも QuadTankPack.mop からの抜粋です。
%%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;
パラメータ推定用モデルをロードします。
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 パラメータで消去することも可能です。
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 に先程行ったパラメータ修正前のシミュレーションの結果ファイルを設定しています。
# 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
変数 temp_1 と temp_2 の initial trajectory や nominal trajectory が見つからないという Warning が出ていますが、このような名前の変数は使用していません。詳細は不明です。この Warning は、オリジナルのデモ qt_par_est_casadi でも発生します。
推定計算結果のパラメータを表示します。
# 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')
推定計算の結果から、液位とポンプ電圧を抽出します。
# 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"]
プロットして測定データと比較します。
# 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()