伝達関数のパラメタをベイズ推定する試み

今以下の伝達関数のような2次のシステムで表現できるブロックがある。

$$G(s)=\frac{\omega_n^2}{s^2+2\zeta \omega_n s+\omega_n^2}$$

いま2次のシステムに対するステップ応答を考える(ステップの高さは1とした)

これは以下のような逆ラプラス変換計算によって計算できる

$$v_{\text{out}}(t)=\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right]$$

計算結果は以下のようになり、以降この応答波形データを中心に扱う

二次システムのステップ応答

$0<\zeta<1$のとき

$$v_{\text{out}}(t)=1 - \frac{\zeta}{\sqrt{1 - \zeta^2}} e^{-\zeta \omega_n t} \sin \left(\omega_n \sqrt{1 - \zeta^2} t \right) - e^{-\zeta \omega_n t} \cos \left(\omega_n \sqrt{1 - \zeta^2} t \right)$$

1曲線からパラメータの推定

いま、$(t_i,y_i)_{i=1,2,...,N}$が与えられている。この応答が二次システムのステップ応答によるもので、なんかしらの観測誤差が乗っているとして、二次システムのパラメータ$\omega_n$や$\zeta$をベイズ推定することを考えてみた

使用するデータ

真の応答に対して分散未知の正規ノイズが乗っているデータを作成

データ生成プログラム

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import lti, step
import japanize_matplotlib

# システムのパラメータ
omega_n = 2.0  # 自然角振動数
zeta = 0.5     # 減衰比

# 伝達関数の分子と分母
numerator = [omega_n**2]
denominator = [1, 2*zeta*omega_n, omega_n**2]

# 伝達関数を作成
system = lti(numerator, denominator)

# シミュレーションの時間ベクトル
t = np.linspace(0, 10, 50)

# ステップ入力信号
u = np.ones_like(t)

# シミュレーション実行
t, y = step(system, T=t)

sigma=0.02

y += np.random.randn(*t.shape)*sigma

def calculated_step_response(t, omega_n, zeta):
    omega_d = omega_n * np.sqrt(1 - zeta**2)  # 減衰振動数
    exp_term = np.exp(-zeta * omega_n * t)
    return 1 - (zeta / np.sqrt(1 - zeta**2)) * exp_term * np
.sin(omega_d * t) - exp_term * np.cos(omega_d * t)

# 修正された応答式によるステップ応答
y_calc = calculated_step_response(t, omega_n, zeta)

# グラフのプロット
plt.figure()
plt.plot(t, y, "b.-",label="observed")

plt.plot(t, y_calc, color="red",label='ideal', linestyle='-')

plt.xlabel('時間')
plt.ylabel('応答')
plt.title('システムのステップ応答')
plt.legend()
plt.grid()
plt.show()

omega_n = 2.0  # 自然角振動数
zeta = 0.5     # 減衰比

を真値のパラメータとしてデータ生成。(以降このパラメータを推定する)

モデリング

以下の仮定のもと、推定対象とするパラメータ$\zeta$,$\omega_n$を推定する

観測誤差の仮定

• 観測誤差は正規分布に従うとした。すなわち

$$y_{\text{observed}}\sim \mathcal{N}(\mu,\sigma^2)$$

• 観測データの平均は$v_{\text{out}}(t)$である。すなわち

$$\mu=1 - \frac{\zeta}{\sqrt{1 - \zeta^2}} e^{-\zeta \omega_n t} \sin \left(\omega_d t \right) - e^{-\zeta \omega_n t} \cos \left(\omega_d t \right)$$

ただし、

$$\omega_d = \omega_n {\sqrt{1 - \zeta^2}}$$

事前分布の設定

• $\omega_n$ は0から3の範囲で一様分布に従うものとする。正の数という情報以外に得られる情報はないとみて暫定的に当範囲に設定

$$\omega_n \sim \mathcal{U}(0,3)$$

• 減衰比$\zeta$ は0から1の範囲で一様分布に従う。これは波形が減衰振動であることから当範囲でカバーできるとみて当範囲に設定

$$\zeta \sim \mathcal{U}(0,1)$$

上の組み立てた仮説および事前分布をpymc
によって以下のように記述することができる

python

import pymc as pm
import numpy as np

# 観測データ
y_observed = y

with pm.Model() as model:
    # 角周波数ω_nの事前分布
    omega_n = pm.Uniform('omega_n', lower=0, upper=3)

    # 減衰比ζの事前分布
    zeta = pm.Uniform('zeta', lower=0, upper=1)

    # 伝達関数の計算
    # 丁寧にDeterministicを使用。npを使わずpmで完結させた
    omega_d = pm.Deterministic('omega_d', omega_n * pm.math.sqrt(1 - zeta**2))
    exp_term = pm.Deterministic("exp_term", pm.math.exp(-zeta * omega_n * t))
    mu = pm.Deterministic('mu', 1 - (zeta / pm.math.sqrt(1 - zeta**2)) * exp_term * pm.math.sin(omega_d * t) - exp_term * pm.math.cos(omega_d * t))

    # 誤差のモデル
    sigma = pm.HalfNormal('sigma', sigma=10)

    # 観測データとの関連
    Y_obs = pm.Normal('Y_obs', mu=mu, sigma=sigma, observed=y_observed)

    # サンプリング
    trace = pm.sample(3000, chains=5)

pm.model_to_graphviz(model)

推定結果の確認

pm.plot_posterior(trace,var_names=["omega_n","zeta"])

サンプルされた曲線を図示

omega_n_posterior=trace.posterior["omega_n"].values.reshape(-1)
zeta_posterior=trace.posterior["zeta"].values.reshape(-1)

plt.figure()
plt.plot(t, y_observed, "b.",label="sample")

for om,zt in zip(omega_n_posterior,zeta_posterior):

  y_calc = calculated_step_response(t, om, zt)

  plt.plot(t, y_calc, color="red", linestyle='-',alpha=0.01)

plt.xlabel('時間')
plt.ylabel('応答')
plt.title('システムのステップ応答')
plt.legend()
plt.grid()
plt.show()

予測区間、平均のベイズ信用区間の分布

プログラム

# 各サンプリングされたパラメータセットに対して予測を生成
pred = []
for omega_n_sample, zeta_sample in zip(trace.posterior['omega_n'].values.flatten(), 
                                       trace.posterior['zeta'].values.flatten()):
    omega_d_sample = omega_n_sample * np.sqrt(1 - zeta_sample**2)
    exp_term_sample = np.exp(-zeta_sample * omega_n_sample * t)
    y_pred_sample = 1 - (zeta_sample / np.sqrt(1 - zeta_sample**2)) * exp_term_sample * np.sin(omega_d_sample * t) - exp_term_sample * np.cos(omega_d_sample * t)
    pred.append(y_pred_sample)

# 予測をnumpy配列に変換
pred = np.array(pred)

# 信用区間を計算
y1, med, y2 = np.percentile(pred, [2.5, 50, 97.5], axis=0)

import plotly
import plotly.graph_objects as go

# 信用区間と中央値の予測をプロットするためのデータを準備
trace1 = go.Scatter(x=t, y=med, mode='lines', name='中央値の予測', line=dict(color='blue'))
trace2 = go.Scatter(x=t, y=y1, mode='lines', name='下限信用区間', line=dict(color='blue'), showlegend=False)
trace3 = go.Scatter(x=t, y=y2, mode='lines', name='上限信用区間', line=dict(color='blue'), fill='tonexty', fillcolor='rgba(0,0,255,0.2)', showlegend=False)
trace4 = go.Scatter(x=t, y=y_observed, mode='markers', name='観測データ', marker=dict(color='black'))

# プロットのレイアウトを設定
layout = go.Layout(
    title='PyMC 5による応答のベイズ信用区間95%範囲',
    xaxis_title='時間 (t)',
    yaxis_title='応答 (y)',
    legend_title='凡例',
    showlegend=True
)

# プロットを作成
fig = go.Figure(data=[trace1, trace2, trace3, trace4], layout=layout)

# プロットを表示
fig.show()

制御対象の個体差ばらつきを表現する

先ほどと状況が違い、複数の曲線データを利用する。

仮定

• 制御対象は同じ製造工程で製造されたものの、多少部品ごとのばらつきの影響上$\zeta$や$\omega_n$などのパラメータが多少部品ごとでばらついてしまうことを想定

こうした個体差のばらつきを考慮したモデリングを行いたい。

使用するデータ

今度は5つのパラメータに未知のばらつきを与えて5本の応答波形を作成。

プログラム

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
np.random.seed(42)  # 乱数生成のためのシードを設定

def calculated_step_response(t, omega_n, zeta):
    """
    指定された自然角振動数と減衰比に基づいてステップ応答を計算する関数。
    """
    omega_d = omega_n * np.sqrt(1 - zeta**2)  # 減衰振動数
    exp_term = np.exp(-zeta * omega_n * t)
    return 1 - (zeta / np.sqrt(1 - zeta**2)) * exp_term * np.sin(omega_d * t) - exp_term * np.cos(omega_d * t)

def generate_step_response_with_noise(t, omega_n_mean, omega_n_std, zeta_mean, zeta_std, noise_std):
    """
    ノイズを含むステップ応答を生成する関数。
    """
    omega_n = np.random.normal(omega_n_mean, omega_n_std)
    zeta = np.random.normal(zeta_mean, zeta_std)
    y = calculated_step_response(t, omega_n, zeta)
    y += np.random.normal(0, noise_std, len(t))
    return y,omega_n,zeta

# シミュレーションの時間ベクトル
t = np.linspace(0, 10, 50)


# N_sampleの曲線データを生成
N_sample = 5

data_sets = []
omega_n_fluctuated=[]
zeta_fluctuated=[]
for _ in range(N_sample):
    y,om,zt = generate_step_response_with_noise(t, omega_n_mean=2.0, omega_n_std=0.1, zeta_mean=0.5, zeta_std=0.1, noise_std=0.02)
    data_sets.append(y)
    omega_n_fluctuated.append(om)
    zeta_fluctuated.append(zt)

# プロット
plt.figure(figsize=(10, 6))
for i, y in enumerate(data_sets):
    plt.plot(t, y, label=f'データセット {i+1} $\omega_n$:{omega_n_fluctuated[i]:.3f}_$\zeta$:{zeta_fluctuated[i]:.3f}')
plt.xlabel('時間')
plt.ylabel('応答')
plt.title('複数のシステムのステップ応答（ノイズ付き）')
plt.legend()
plt.grid()
plt.show()

モデリング

5つの制御パラメータがそれぞれ共通の事前分布によってサンプルされているよう扱うことで個体ごとのばらつきが表現できると考えた。(結果的に階層ベイズモデルになる)

import pymc as pm

cl = [i for i in range(N_sample)]

%%time

with pm.Model() as model:
    cl_data = pm.ConstantData('cl_data', np.arange(N_sample))

    # 共通の事前分布の平均と標準偏差
    omega_n_mu = pm.Normal('omega_n_mu', mu=1.5, sigma=0.5)
    omega_n_sd = pm.HalfNormal('omega_n_sd', sigma=0.5)
    zeta_mu = pm.Normal('zeta_mu', mu=0.5, sigma=0.2)
    zeta_sd = pm.HalfNormal('zeta_sd', sigma=0.2)

    # 各データセットに対するパラメータ
    omega_ns = pm.Normal('omega_ns', mu=omega_n_mu, sigma=omega_n_sd, shape=N_sample)
    zetas = pm.Normal('zetas', mu=zeta_mu, sigma=zeta_sd, shape=N_sample)

    # 伝達関数の計算
    omega_ds = pm.Deterministic('omega_ds', omega_ns[cl_data] * pm.math.sqrt(1 - zetas[cl_data]**2))
    exp_terms = pm.Deterministic("exp_terms", pm.math.exp(-zetas[cl_data] * omega_ns[cl_data] * t[:, None]))
    mu = pm.Deterministic('mu', 1 - (zetas[cl_data] / pm.math.sqrt(1 - zetas[cl_data]**2)) * exp_terms * pm.math.sin(omega_ds * t[:, None]) - exp_terms * pm.math.cos(omega_ds * t[:, None]))

    # 誤差のモデル
    sigma = pm.HalfNormal('sigma', sigma=10)

    # 観測データとの関連
    for i in range(N_sample):
        Y_obs = pm.Normal(f'Y_obs_{i}', mu=mu[:, i], sigma=sigma, observed=data_sets[i])

    # サンプリング
    trace = pm.sample(random_seed=42,target_accept=0.998)

pm.model_to_graphviz(model)

import matplotlib.pyplot as plt
import numpy as np
# import pymc3 as pm

# トレースのプロット
axes = pm.plot_trace(trace, var_names=["omega_ns", "zetas"])

# Matplotlibのデフォルトの色サイクルを取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

# 真値の縦線を追加
for i in np.arange(N_sample):
    color = colors[i % len(colors)]  # 色をループで回す
    axes[0, 0].axvline(omega_n_fluctuated[i], color=color, linestyle='--')
    axes[0, 1].axvline(omega_n_fluctuated[i], color=color, linestyle='--')
    axes[1, 0].axvline(zeta_fluctuated[i], color=color, linestyle='--')
    axes[1, 1].axvline(zeta_fluctuated[i], color=color, linestyle='--')

plt.tight_layout()

真値とは多少ずれているように見えているが概ね推定できている模様