[Python]伝達関数モデルの同定

記事の内容

Pythonでシステム同定を行うサンプルスクリプト。
モデルは伝達関数。
同定用入力信号にステップ入力およびChirp信号を用いてみたが、単純な伝達関数ではどちらの信号を用いても問題なくシステム同定できた。

スクリプトの例１

モデル：伝達関数
同定用入力信号：ステップ入力
検証用入力信号：ステップ入力

a03_SystemIdentification.ipynb

from control.matlab import*
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# -------入力信号の生成------------
duration = 100  # 信号の継続時間 (秒)
sampling_rate = 100  # サンプリングレート (Hz)

# 時間軸の作成
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)

#入力信号の生成（入力信号:ステップ入力)
u = np.zeros_like(t)
u[int(1 * sampling_rate):] = 1

#-------伝達関数モデル(実)の定義----------
s = tf('s')
params = [1,2,1]
print("実モデルのパラメータ",params)
a, b, c = params
P = 1/(a*s**2+b*s+c)

#-----------lsimで出力（応答）計算---------
def calc_y(P, u, t):
    x0 = [0,0] #初期値
    y, t, _ = lsim(P, U=u, T=t, X0=x0)
    return y

y_act = calc_y(P,u,t)

#---------パラメータ推定の推定を行う-------------
#目的関数
def objective_function(params, t, u, y_act):
    s = tf('s')
    a, b, c = params
    P_pred = 1/(a*s**2+b*s+c)
    y_pred = calc_y(P_pred, u, t)
    sse = sum((y_act - y_pred)**2)
    return sse

#-------- 最適化を持ちいてパラメータを推定する------------------
initial_params = [1, 1, 1]  # 初期推定値

#最適化実施
result = minimize(objective_function, initial_params, args=(t, u, y_act))
print('推定さえたモデルのパラメータ',result.x)
a_pred, b_pred, c_pred = result.x

#応答の予測値
P_pred = 1/(a_pred*s**2+b_pred*s+c_pred)

#-----------lsimで出力（応答）計算---------
y_pred = calc_y(P_pred, u, t)

#---------入出力信号のグラフ表示---------
plt.figure(figsize=(10, 10))

#全体表示
plt.subplot(2, 1, 1)
plt.axhline(0, color="gray", linestyle="--")
plt.plot(t, u, label="Input")
plt.plot(t, y_act, label="Output")
plt.plot(t, y_pred, linestyle='--', label="Output_Prediction")
plt.xlabel("t (s)")
plt.ylabel("Signal [-]")
plt.title("Response")
plt.legend()
plt.grid(True)

#拡大表示
plt.subplot(2, 1, 2)
plt.axhline(0, color="gray", linestyle="--")
plt.plot(t, u, label="Input")
plt.plot(t, y_act, label="Output")
plt.plot(t, y_pred, linestyle='--', label="Output_Prediction")
plt.xlim(0,20)
plt.xlabel("Time (s)")
plt.ylabel("Signal [-]")
plt.title("Response")
plt.legend()
plt.grid(True)

plt.show()

実行結果

スクリプトの例２

モデル：伝達関数
同定用入力信号：Chirp信号
検証用入力信号：ステップ入力

a03_SystemIdentification.ipynb

from control.matlab import*
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# -------入力信号の生成------------
duration = 100  # 信号の継続時間 (秒)
sampling_rate = 100  # サンプリングレート (Hz)
frequency_start = 0.1  # チャープ信号の開始周波数 (Hz)
frequency_end = 1 # チャープ信号の終了周波数 (Hz)
# 時間軸の作成
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)


#訓練データの入力信号(chirp信号)
u_train = np.sin(2 * np.pi * np.cumsum(np.linspace(frequency_start, frequency_end, len(t)) / sampling_rate)) 

#検証データの信号(step入力)
u_vali = np.zeros_like(t)
u_vali[int(1 * sampling_rate):] = 1


#-------伝達関数モデル(実)の定義----------
s = tf('s')
a = 1
b = 2
c = 1
P = 1/(a*s**2+b*s+c)

#-----------lsimで出力（応答）計算---------
def calc_y(P, u, t):
    x0 = [0,0] #初期値
    y, t, _ = lsim(P, U=u, T=t, X0=x0)
    return y

#訓練データ
y_train_act = calc_y(P,u_train,t)

#---------パラメータ推定の推定を行う-------------
#目的関数
def objective_function(params, t, u, y_train_act):
    s = tf('s')
    a, b, c = params
    P_pred = 1/(a*s**2+b*s+c)
    y_train_pred = calc_y(P_pred, u, t)
    sse = sum((y_train_act - y_train_pred)**2)
    return sse

# 最適化を持ちいてパラメータを推定
initial_params = [1, 1, 1]  # 初期推定値

#最適化実施
result = minimize(objective_function, initial_params, args=(t, u_train, y_train_act))
print('予測されたパラメータ',result.x)
a_pred, b_pred, c_pred = result.x

#予測されたモデル
P_pred = 1/(a_pred*s**2+b_pred*s+c_pred) 


#------------訓練データに対する予測値の確認------------------
#応答計算
y_train_act = calc_y(P, u_train, t)
y_train_pred = calc_y(P_pred, u_train, t)

#グラフ表示
print("訓練データ");
plt.figure(figsize=(10, 10))

plt.subplot(2, 1, 1)
plt.axhline(0, color="b", linestyle="--")
plt.plot(t, u_train, label="Input")
plt.plot(t, y_train_act, label="Output")
plt.plot(t, y_train_pred, linestyle='--', label="Output_Prediction")
plt.xlabel("t (s)")
plt.ylabel("Signal [-]")
plt.title("Training Data")
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.axhline(0, color="b", linestyle="--")
plt.plot(t, u_train, label="Input")
plt.plot(t, y_train_act, label="Output")
plt.plot(t, y_train_pred, linestyle='--', label="Output_Prediction")
plt.xlim(0,20)
plt.xlabel("Time (s)")
plt.ylabel("Signal [-]")
plt.title("Training Data")
plt.legend()
plt.grid(True)

plt.show()

#------------検証データに対する予測値の確認------------------
#応答計算
y_vali_act = calc_y(P, u_vali, t)
y_vali_pred = calc_y(P_pred, u_vali, t)

#グラフ表示
print("検証データ");
plt.figure(figsize=(10, 10))

plt.subplot(2, 1, 1)
plt.axhline(0, color="b", linestyle="--")
plt.plot(t, u_vali, label="Input")
plt.plot(t, y_vali_act, label="Output")
plt.plot(t, y_vali_pred, linestyle='--', label="Output_Prediction")
plt.xlabel("t (s)")
plt.ylabel("Signal [-]")
plt.title("Validation Data")
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.axhline(0, color="b", linestyle="--")
plt.plot(t, u_vali, label="Input")
plt.plot(t, y_vali_act, label="Output")
plt.plot(t, y_vali_pred, linestyle='--', label="Output_Prediction")
plt.xlim(0,20)
plt.xlabel("Time (s)")
plt.ylabel("Signal [-]")
plt.title("Validation Data")
plt.legend()
plt.grid(True)

plt.show()