概要

simulinkによる1次遅れシステムのシミュレーション結果を用いて、システム同定を行った。
具体的な流れを示す。
1：1次遅れ系伝達関数を作成
2：作成した伝達関数の時定数からシステム同定するためのM系列信号を生成
3：M系列信号を入力(input)としてsimulinkにより伝達関数の出力(output)を保存
4：inputとoutputを用いて、System Identificationによりシステム同定
5：システム同定結果と1で作成した伝達関数を比較

システム同定した伝達関数と1で作成した伝達関数を比較したところほぼ一致していることを確認した。
詳細を以下に示す。

使用ソフト

MATLAB バージョン 9.9 (R2020b)
Simulink バージョン 10.2 (R2020b)
System Identification Toolbox バージョン 9.13 (R2020b)

1：1次遅れ系伝達関数を作成

機械的モデルをベースに以下1次遅れ伝達関数を設定した

G(s) = \frac{1}{\frac{J}{D}s+1}

ただし
D = 0.0005
J = 0.001

時定数はtau = J/Dより
2sとなる
実際にシミュレーションしてみると確かに時定数は2sであることが確認できる

今回、本伝達関数が同定対象となる

2：作成した伝達関数の時定数からシステム同定するためのM系列信号を生成

同定対象の伝達関数の時定数は2sであるためM系列信号の条件は以下に設定した
n = 9; %シフトレジスタ数
N = 2^n; %M系列信号1周期当たりのデータ数
Ts = 0.1; %データサンプリング周期[s]
p = 10; %Mクロック周期の倍数：
Tm = Ts * p; %M系列サンプリング周期[s]

条件決めにあたっての方向性は以下2つの制約に基づいている
制約1：
M系列信号の入力によって同定対象の応答は定常値になる必要がある。
つまり以下関係を満足する必要がある

nTm > tau

tau：時定数
※nTmはM系列信号が最大で同じ値をとる時間になる。

制約2：
同定対象のステップ応答の間に5~8サンプル点が入るくらいの間隔をサンプリング周期Tsとする

同定対象の時定数は2sであることからデータサンプリング周期Tsは0.1とした。
※20サンプルも入ってしまうという突っ込みはあるが。。。
シフトレジスタ数n = 9は決め打ち。
※参考文献が9であった。
制約1からTmを1に設定

補足
Tm=1なのでサンプリング周期を10倍(p)させた値である。10以下にするとM系列信号が
すぐに切り替わり定常値になりずらかった。

M系列信号生成コードを以下に示す。

n = 9; %シフトレジスタ数
N = 2^n; %M系列信号1周期当たりのデータ数
Ts = 0.1; %データサンプリング周期[s]
p = 10; %Mクロック周期の倍数
Tm = Ts * p; %M系列サンプリング周期[s]


u = idinput(N,'prbs',[n,1/p],[-1,1]);
t =0:Tm:( length (u) -1) *Tm;
subplot(311)
stairs(t,u), axis([t(1) t(end) -1.2 1.2]);
xlabel('time[sec]')
ylabel('信号')
title('M系列信号');

subplot(312),plot(covf(u, N)); axis([0 N -0.5 1]);
xlabel('ラグ')
ylabel('相関')
title('自己相関関数');
subplot(313), periodogram(u);
title('パワースペクトル');

%ary_m_signal = [t,transpose(u)]; %配列生
m_signal = timeseries(transpose(u),t);
save('mat_m_signal','m_signal','-v7.3'); %matファイル生成

3：M系列信号入力(input)と伝達関数の出力(output)を保存

simulink(下図)よりinputとoutputを保存
Tsは0.1なのでzohはTsにしている

ワークスペース保存ブロックを用いた。
※ワークスペース保存ブロック設定画面

保存形式は2次元配列にすることで以下のようにワークスペースに保存される

上画面はoutputのみだがinputも同様に保存できる。

また、inputとoutputの波形は以下となった

4：inputとoutputを用いて、System Identificationによるシステム同定

inputとoutputデータをSystem Identificationに読み込ませる。

systemIdentification

inputとoutputはsimulinkで出力した変数名をいれる
Startは0
SampleTimeは0.1

伝達関数による推定を選択する
Estimate->Transfer function
同定対象の伝達関数は1次遅れ系なので
ゼロ点の数は0
極の数は1
で設定し同定を開始する

5：システム同定結果と1で作成した伝達関数を比較

同定結果を示す。

今回1で設定した伝達関数は

G(s) = \frac{1}{\frac{J}{D}s+1}

J/D = 2なので

G(s) = \frac{0.5}{s+0.5}

1で設定した伝達関数とほぼ一致していることがわかる

参考

本：Matlabによる制御のためのシステム同定
本：https://booth.pm/ja/items/1180601

今後

いったんシステム同定系はここで終了する