Notebook

古田の振り子のパラメータのキャリブレーションを行う¶

JModelica.org Users Guide - Version 2.10 6.6. Derivative-Free Model Calibration of FMUs の実習を行います。本文書は Users Guid の例題を解説するものであり、モデルのソースコードやオペレーションで使用する Python コードなどの著作権は Modelon AB にあります。

概要¶

古田の振り子 (Furuta pendulum) は、水平面内で回転するアームとアームに対して垂直な面内で自由に回転する振り子から構成されています。アームの回転角度 $\varphi$ とアームと振り子の間の回転角度 $\theta$ の２つの自由度があります。実験で測定した $\varphi$ と $\theta$ の値を用いて、モデルのパラメータであるアームの摩擦トルク係数 armFriction と振り子の摩擦トルク係数 pendulumFriction を求めます。最適化アルゴリズムはNeider-Mead simplex optimization algorithm を使用します。

準備¶

必要なものは、

FurutaData.mat : $\varphi$ と $\theta$ の測定データ
Furuta.fmu : 古田の振り子のモデルのFMU

です。測定データは、MATLABの matファイルになっていますが、Python で時系列データを読み込めれば他のファイル形式でもOKです。

Windows の場合は、以下のようにして作業フォルダにコピーします。

%copy C:\JModelica.org-2.14\install\Python_64\pyjmi\examples\files\FurutaData.mat .
%copy C:\Jmodelica.org-2.14\install\Python_64\pyjmi\examples\files\FMUS\Furuta.fmu .

MacでDockerを使用している場合は、FMUが提供されていないのでソースコードからコンパイルして作成します。

%cp /usr/local/jmodelica/Python/pyjmi/examples/files/FurutaData.mat .
%cp /usr/local/jmodelica/Python/pyjmi/examples/files/Furuta.mo .
from pymodelica import compile_fmu
fmu = compile_fmu('Furuta', 'Furuta.mo')

In [1]:

%cp /usr/local/jmodelica/Python/pyjmi/examples/files/FurutaData.mat .
%cp /usr/local/jmodelica/Python/pyjmi/examples/files/Furuta.mo .
from pymodelica import compile_fmu
fmu = compile_fmu('Furuta', 'Furuta.mo')

必要なモジュールをインポートします。

In [2]:

%matplotlib notebook
from scipy.io.matlab.mio import loadmat
import matplotlib.pyplot as plt
import numpy as N
from pyfmi import load_fmu
from pyjmi.optimization import dfo

測定値とキャリブレーション前のモデルを比較する¶

まず、測定値を読み込みます。

In [3]:

# Load measurement data from file
data = loadmat('FurutaData.mat',appendmat=False)
# Extract data series
t_meas = data['time'][:,0]
phi_meas = data['phi'][:,0]
theta_meas = data['theta'][:,0]

キャリブレーション前のモデルのシミュレーションを実行して結果を取得します。

In [4]:

# Load model
model = load_fmu("Furuta.fmu")
# Simulate model response with nominal parameters
res = model.simulate(start_time=0.,final_time=40)
# Load simulation result
phi_sim = res['armJoint.phi']
theta_sim = res['pendulumJoint.phi']
t_sim = res['time']

Final Run Statistics: --- 

 Number of steps                                 : 1416
 Number of function evaluations                  : 2232
 Number of Jacobian evaluations                  : 90
 Number of function eval. due to Jacobian eval.  : 360
 Number of error test failures                   : 70
 Number of nonlinear iterations                  : 1872
 Number of nonlinear convergence failures        : 0
 Number of state function evaluations            : 2205
 Number of state events                          : 89

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 - 40.0 seconds.
Elapsed simulation time: 0.217043876648 seconds.

実験による測定値とシミュレーション結果をプロットします。

In [5]:

# Plot measurements
plt.figure (1)
plt.subplot(2,1,1)
plt.plot(t_meas,theta_meas,label='Measurements')
plt.plot(t_sim,theta_sim,'--',label='Simulation nominal parameters')
plt.ylabel(r"$\theta$ [rad]")
plt.legend(loc=1)
plt.grid ()
plt.subplot(2,1,2)
plt.plot(t_meas,phi_meas,label='Measurements')
plt.plot(t_sim,phi_sim,'--',label='Simulation nominal parameters')
plt.ylabel(r"$\varphi$ [rad]")
plt.legend(loc=1)
plt.grid ()
plt.show ()

目的関数を作成する¶

最適化アルゴリズムは、目的関数 (objective function)

$f(x) = \sum_{i=1}^M{\left( \varphi^{sim}(t_i, x)-\varphi^{meas}(t_i) \right)^2} + \sum_{i=1}^M{\left( \theta^{sim}(t_t, x)-\theta^{meas}(t_i) \right)^2}$

を最小化するパラメータ $x$ を求めます。 $\varphi$ と $\theta$ に付けられた上付きの添え字 $sim$ と $meas$ は、それぞれ、シミュレーションの結果と測定結果を表しています。 $t_i$ は測定時刻を表しています。

まず、パラメータの初期値、下限値、上限値の配列を作成します。このパラメータは、最適化アルゴリズムが適切に動作する大きさになるようなスケールファイクター $=10^3$ をかけています。

In [6]:

# Choose starting point
x0 = N.array([0.012,0.002])*1e3
# Choose lower and upper bounds (optional)
lb = N.zeros (2)
ub = (x0 + 1e-2)*1e3

目的関数を含む Python のスクリプトファイルを作成します。独立したファイルにする理由は、パラメータスタディをマルチプロセスで並列的に行うためです。

スクリプトファイル名は、目的関数名.py とします。
目的関数の引数 $x$ は、スケーリングしたパラメータの配列です。
スクリプトを読み込むと、必要なモジュールと実験による測定値が読み込まれます。

目的関数内の主な処理内容を示します。

引数 $x$ をスケールファクターで割ってパラメータの値を得ます。
FMU をロードします。
パラメータとシミュレーションオプションを設定してシミュレーションを実行します。
シミュレーション結果を取得します。
目的関数の値を計算します。

def furuta_cost(x): から始まる部分が目的関数です。目的関数の計算は、作業フォルダ内のサブディレクトリ dir_1, dir_2,... で並列的に実行されるので、作業フォルダにある Furuta.fmu は

model = load_fmu('../Furuta.fmu')

のようにして読み込みます。

In [7]:

%%writefile furuta_cost.py
# Import Modules
from pyfmi import load_fmu
from pyjmi.optimization import dfo
from scipy.io.matlab.mio import loadmat
import numpy as N

# Load measurement data from file
data = loadmat('FurutaData.mat',appendmat=False)
t_meas = data['time'][:,0]
phi_meas = data['phi'][:,0]
theta_meas = data['theta'][:,0]

# Define the objective function
def furuta_cost(x):
    # Scale down the inputs x since they are scaled up
    # versions of the parameters (x = 1e3*[param1,param2])
    armFrictionCoefficient = x[0]/1e3
    pendulumFrictionCoefficient = x[1]/1e3
    # Load model
    model = load_fmu('../Furuta.fmu')
    # Set new parameter values into the model 
    model.set('armFriction', armFrictionCoefficient)
    model.set('pendulumFriction', pendulumFrictionCoefficient)
    # Create options object and set verbosity to zero to disable printouts
    opts = model.simulate_options()
    opts['CVode_options']['verbosity'] = 0
    # Simulate model response with new parameter values
    res = model.simulate(start_time=0., final_time=40, options=opts)
    # Load simulation result
    phi_sim = res['armJoint.phi']
    theta_sim = res['pendulumJoint.phi']
    t_sim  = res['time']
    # Evaluate the objective function
    y_meas = N.vstack((phi_meas, theta_meas))
    y_sim = N.vstack((phi_sim, theta_sim))
    obj = dfo.quad_err(t_meas, y_meas, t_sim, y_sim)
    return obj

Overwriting furuta_cost.py

最適化問題を解く¶

最適化アルゴリズムよる計算を実行します。

In [8]:

# Solve the problem using the Nelder-Mead simplex algorithm
x_opt,f_opt,nbr_iters,nbr_fevals,solve_time = dfo.fmin("furuta_cost.py",
xstart=x0,lb=lb,ub=ub,alg=1,nbr_cores=4,x_tol=1e-3,f_tol=1e-2)

 
Number of iterations: 0
Number of function evaluations: 3
Current time: 0.016255 s
 
Current x value: [ 12.   2.]
Current function value: 209.331105292
 
Termination criterion for x: 3.64965751818
Termination criterion for the objective function: 123.157423108
 
 
Number of iterations: 1
Number of function evaluations: 10
Current time: 0.044058 s
 
Current x value: [ 12.   2.]
Current function value: 209.331105292
 
Termination criterion for x: 1.82482875909
Termination criterion for the objective function: 119.22142411
 
 
Number of iterations: 2
Number of function evaluations: 17
Current time: 0.071913 s
 
Current x value: [ 12.   2.]
Current function value: 209.331105292
 
Termination criterion for x: 1.82482875909
Termination criterion for the objective function: 95.8448736809
 
 
Number of iterations: 3
Number of function evaluations: 24
Current time: 0.098484 s
 
Current x value: [ 13.36862157   1.08758562]
Current function value: 38.360389733
 
Termination criterion for x: 1.64487841496
Termination criterion for the objective function: 190.054639792
 
 
Number of iterations: 4
Number of function evaluations: 31
Current time: 0.126799 s
 
Current x value: [ 11.54379281   1.08758562]
Current function value: 20.5482250917
 
Termination criterion for x: 1.82482875909
Termination criterion for the objective function: 188.782880201
 
 
Number of iterations: 5
Number of function evaluations: 38
Current time: 0.153931 s
 
Current x value: [ 11.54379281   1.08758562]
Current function value: 20.5482250917
 
Termination criterion for x: 1.82482875909
Termination criterion for the objective function: 107.524431612
 
 
Number of iterations: 6
Number of function evaluations: 45
Current time: 0.182587 s
 
Current x value: [ 12.57025899   0.85948203]
Current function value: 18.7380020422
 
Termination criterion for x: 1.05150561696
Termination criterion for the objective function: 19.6223876908
 
 
Number of iterations: 7
Number of function evaluations: 52
Current time: 0.213943 s
 
Current x value: [ 10.74543023   0.85948203]
Current function value: 7.53256646059
 
Termination criterion for x: 1.82482875909
Termination criterion for the objective function: 13.0156586311
 
 
Number of iterations: 8
Number of function evaluations: 59
Current time: 0.243384 s
 
Current x value: [ 10.74543023   0.85948203]
Current function value: 7.53256646059
 
Termination criterion for x: 1.82482875909
Termination criterion for the objective function: 11.2054355816
 
 
Number of iterations: 9
Number of function evaluations: 66
Current time: 0.270603 s
 
Current x value: [ 9.77598995  0.97353382]
Current function value: 1.96608170793
 
Termination criterion for x: 1.82482875909
Termination criterion for the objective function: 6.64950588566
 
 
Number of iterations: 10
Number of function evaluations: 73
Current time: 0.297702 s
 
Current x value: [ 9.77598995  0.97353382]
Current function value: 1.96608170793
 
Termination criterion for x: 1.1551264072
Termination criterion for the objective function: 5.56648475266
 
 
Number of iterations: 11
Number of function evaluations: 80
Current time: 0.328624 s
 
Current x value: [ 10.15735065   1.00917501]
Current function value: 1.61555414395
 
Termination criterion for x: 0.776069960793
Termination criterion for the objective function: 1.59256106871
 
 
Number of iterations: 12
Number of function evaluations: 87
Current time: 0.35713 s
 
Current x value: [ 10.15735065   1.00917501]
Current function value: 1.61555414395
 
Termination criterion for x: 0.672748641737
Termination criterion for the objective function: 0.350527563986
 
 
Number of iterations: 13
Number of function evaluations: 94
Current time: 0.387477 s
 
Current x value: [ 9.79848845  0.99269096]
Current function value: 1.30125665815
 
Termination criterion for x: 0.359240589922
Termination criterion for the objective function: 0.379279413464
 
 
Number of iterations: 14
Number of function evaluations: 101
Current time: 0.417535 s
 
Current x value: [ 9.7312714   1.00772709]
Current function value: 1.22584476847
 
Termination criterion for x: 0.426081709334
Termination criterion for the objective function: 0.389709375475
 
 
Number of iterations: 15
Number of function evaluations: 108
Current time: 0.443055 s
 
Current x value: [ 9.96111529  1.00469202]
Current function value: 1.16460504035
 
Termination criterion for x: 0.229863925257
Termination criterion for the objective function: 0.136651617807
 
 
Number of iterations: 16
Number of function evaluations: 115
Current time: 0.466012 s
 
Current x value: [ 9.8223409   0.99945026]
Current function value: 1.15532030535
 
Termination criterion for x: 0.138873350249
Termination criterion for the objective function: 0.0705244631182
 
 
Number of iterations: 17
Number of function evaluations: 122
Current time: 0.491536 s
 
Current x value: [ 10.05218478   0.99641519]
Current function value: 1.14559534571
 
Termination criterion for x: 0.229863925257
Termination criterion for the objective function: 0.0190096946348
 
 
Number of iterations: 18
Number of function evaluations: 129
Current time: 0.517574 s
 
Current x value: [ 10.05218478   0.99641519]
Current function value: 1.14559534571
 
Termination criterion for x: 0.229863925257
Termination criterion for the objective function: 0.00972495964301
 
 
Solver: Nelder-Mead
 
Optimization terminated successfully.
Terminated due to sufficiently close function values at the vertices of the simplex.
 
Total number of iterations: 18
Total number of function evaluations: 130
Total execution time: 0.523607 s

最適化の結果¶

In [9]:

# Optimal point (don't forget to scale down)
[armFrictionCoefficient_opt, pendulumFrictionCoefficient_opt] = x_opt/1e3
# Print optimal parameter values and optimal function value
print('Optimal parameter values:')
print('arm friction coeff = ' + str(armFrictionCoefficient_opt))
print('pendulum friction coeff = ' + str(pendulumFrictionCoefficient_opt))
print('Optimal function value: ' + str(f_opt))

Optimal parameter values:
arm friction coeff = 0.0100521847827
pendulum friction coeff = 0.000996415186427
Optimal function value: 1.14559534571

測定値とキャリブレーション後のモデルを比較する¶

最適化したパラメータを設定したモデルのシミュレーションを行います。

In [10]:

# Load model
model = load_fmu("Furuta.fmu")
# Set optimal parameter values into the model
model.set('armFriction',armFrictionCoefficient_opt)
model.set('pendulumFriction',pendulumFrictionCoefficient_opt)
# Simulate model response with optimal parameter values
res = model.simulate(start_time=0.,final_time=40)
# Load simulation result
phi_opt = res['armJoint.phi']
theta_opt = res['pendulumJoint.phi']
t_opt = res['time']

Final Run Statistics: --- 

 Number of steps                                 : 2246
 Number of function evaluations                  : 3538
 Number of Jacobian evaluations                  : 141
 Number of function eval. due to Jacobian eval.  : 564
 Number of error test failures                   : 118
 Number of nonlinear iterations                  : 2978
 Number of nonlinear convergence failures        : 0
 Number of state function evaluations            : 3499
 Number of state events                          : 139

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 - 40.0 seconds.
Elapsed simulation time: 0.320708036423 seconds.

測定値と最適化したパラメータによる計算結果を比較します。

In [11]:

# Plot simulation result
plt.figure (2)
plt.subplot(2,1,1)
plt.plot(t_meas,theta_meas,label='Measurements')
plt.plot(t_opt,theta_opt,'--',linewidth=2, label='Simulation optimal parameters')
plt.legend(loc=1)
plt.ylabel(r"$\theta$ [rad]")
plt.subplot(2,1,2)
plt.plot(t_meas,phi_meas,label='Measurements')
plt.plot(t_opt,phi_opt,'--',linewidth=2, label='Simulation optimal parameters')
plt.legend(loc=1)
plt.ylabel(r"$\varphi$ [rad]")
plt.show ()

In [ ]: