あらすじ
倒立振子ロボットをPIDで立たそうとしたが、中々立たない。
やはり、現在制御でやってみるのが良さそうだ。
そのために、まず倒立振子の物理モデルを構築する必要がある。
どうせ物理モデルを作るなら、パラメータはえいやで決めるのではなく、
やはり実際の値に近いものを算出したい。
そこで、システム同定で物理モデルのパラメータを算出することにした。
まず気になるのが、制御対象への入力値。
同定方法を書籍やネットで調査すると、単なるステップ信号の他に、
「M系列信号」なるものを使用しているケースが目立つ。
M系列信号とは何なのか？どのような性質を持つのか？どうやって生成するのか？
それらの答えや調査の過程で学んだ考え等を、以下記載していく。
M系列信号とは？
M系列(maximum-length sequence)とは、疑似的な不規則信号のことである。システム同定にはPE性（持続的励振）の観点から白色雑音が望ましいが、理想的な白色雑音は物理的に実現不可なため、本信号を用いる。特に、線形システム同定を行うには、取扱の簡易さから二値信号を用いる場合が多い。生成には、直列に繋げた複数のシフトレジスタと排他的論理和（xor）を用いる。
基本的に、信号の周期は2^[シフトレジスタ数]-1で決定されるが、そのためには
xorの入力をとる場所（タップ位置）に気を付けなければならない。
例えば入力数2のxorを1つだけ用いる場合、2つの入力のうち1つは必ず最下段の
シフトレジスタからとる必要があるが、もう1つは所定のシフトレジスタの出力を入力として
とらないと、2^[シフトレジスタ数]-1の周期を実現できない。
※ 詳細は後述のソース「maximum_length_sequence.m」を参照
シフトレジスタとは？
Dフリップフロップを複数繋げることで、入力信号を左か右に１ビットずつシフトして出力する回路のこと。これを使用すれば、シリアル信号をパラレルの回路に変換することができる。
因みに「4ビットのシフトレジスタ」とは、4つのDフリップフロップを保持している。
下記の記事など、シフトレジスタを有効活用して回路に用いる配線の簡易化に成功している。
https://qiita.com/elyoshio92/items/b207d8e534e7169854f7
Dフリップフロップとは？
入力値を保持し、クロック信号のエッジ立ち上がりのタイミングで出力する回路。
Dとはdelayを意味し、文字通り入力値を「遅れて」出力するイメージとなる。
上の図（simulink）で言うと、①が入力、②がクロック信号、③がDフリップフロップを有効/無効化するためのスイッチ、④が出力、⑤が出力を逆パルスにしたもの。
Dフリップフロップの中身は上記のように、NANDゲート４つを組み合わせたものになっている。
D-Latch（クロックが1の時は入力を素通しし、0の時は値を保持）と混同しないように注意が必要。
因みに、他には以下のようなフリップフロップが存在する。
| 名前 | 機能 |
|---|---|
| RS-FF | Reset-Setの意味。Setに入力で出力が1に、Resetに入力で出力が0になり、切り替えスイッチのように使用できる。ただしSetとReset両方が1の場合、出力が不定になるので注意。 |
| JK-FF | Jack-Knifeの意味。RS-FFがRとS両方1の時に出力が不定になる欠点を補ったもので、両方1の時は入力を反転させる。入力クロックの立下り時に状態が変化する。 |
| T-FF | Toggleの意味。JK-FFの入力JとKをひとまとめにしたもので、入力が1の時に出力が反転する。1を入力し続けると周期毎に出力が切り替わるので注意。 |
因みにT-FFはsimulinkに用意されていないので、XORを使用して以下のように作成する。
M系列信号を実際に生成してみた
では本題に戻り、matlab/Simulinkを使用して実際にM系列信号を生成してみる。
今回は、シフトレジスタ数:7、周期:127(2^[シフトレジスタ数]-1) の
M系列信号を生成する。
3つのグラフは、上から順に以下を意味する。
１．M系信号
２．M系列信号の自己相関関数
３．M系列信号のパワースペクトル
１．M系列信号のみ見ると、様々な周波数の信号を網羅していることがなんとなく分かる。。程度だが
２．自己相関関数を見ると、周期の127（とその倍数）の時のみ1を、それ以外はほぼ0を示しており、
周期127の様々な周波数を含んだ信号であることが分かる。
因みに分かりづらいが、赤い点線は127とその倍数の値の時にプロットされており、自己相関関数が
1になるタイミングとピッタリ一致しているのが分かる。
３．パワースペクトルを見てみても、様々な周波数が含まれており比較的平坦な周波数特性を示していることが分かる。
因みに自己相関関数とは、ある波形とそれが$τ$の時間分だけずれたもの同士がどれぐらい整合するか（＝近いか）、を示してる。理想的な白色雑音の自己相関関数は、$τ$ が0以外の時は0に近い値となるが、前述のように理想的な白色雑音は実現不可であり、$τ$が$0$及び周期以外の時は$0$近くになっていることが確認できたため、OKとした。
ソースファイル
srnum = 7; % シフトレジスタの数。これにより出力信号の周期が設定される
helz = 2^srnum - 1; % 周期
len = 1000; % 長さ
x = 0:(len-1);
% M系列信号を生成＆表示
msig = maximum_length_sequence(len, srnum);
subplot(3, 1, 1);
plot(x, msig);
title('１．M系列信号');
% M系列信号の自己相関関数を生成＆表示
ac = auto_correlation(0:len, msig);
subplot(3, 1, 2);
plot(x(1:length(ac)), ac);
title('２．自己相関関数');
% M系列の周期毎に縦線を表示（自己相関関数が妥当か、の確認のため）
hold on
for i = 1 : len
if helz * i <= (len/2)
plot([helz*i helz*i] - 1, [min(ac) max(ac)], 'r--'); % x軸の-1は、ズレる分の補正(黒線は'k')
else
break;
end
end
hold off
% パワースペクトルを表示
subplot(3, 1, 3);
[f, power] = disp_power_spectrum(msig);
plot(f, power);
title('３．パワースペクトル');
function output = maximum_length_sequence(time, srnum)
% time:時間 srnum:シフトレジスタ数
% シフトレジスタの数毎の、xorのタップ位置（最後のシフトレジスタと、どれかのシフトレジスタのxor）
eachtap = [1 0 0 0 0 0 0 0 0;
1 1 0 0 0 0 0 0 0;
1 0 1 0 0 0 0 0 0;
1 0 0 1 0 0 0 0 0;
0 1 0 0 1 0 0 0 0;
1 0 0 0 0 1 0 0 0;
1 0 0 0 0 0 1 0 0;
0 0 0 0 0 0 0 0 0; % シフトレジスタ数が8の時は、xor1つのみでは（n^2-1周期の信号を生成できない）
0 0 0 1 0 0 0 0 1];
z = eachtap(srnum, 1:srnum); % 選択したシフトレジスタの数に合わせたxorのタップ数を抽出
postap = find(z==1, 1); % 最後のシフトレジスタ以外のタップ位置を抽出
output = 0;
for i = 1 : time
output(i) = xor(z(postap), z(srnum));
% それぞれのシフトレジスタの出力を更新
for j = srnum : -1 : 2
z(j) = z(j - 1);
end
z(1) = output(i); % xor出力と直接繋がるシフトレジスタのみ、xorの出力で初期化
end
output(output == 0) = -1; % -1と1の二値信号にするため、0を-1に置換
% plot(0:length(output) - 1, output); % テスト用
end
function yout = auto_correlation(xin, yin)
T = length(xin)/2;
% τを少しずつずらしていく
for a = 1 : T % a = τ のつもり
s = 0;
for t = 2 : T % τを一定値ずらした状態で、2つの波形の積を求める（t = time）
s = s + yin(t) * yin(t + a); % 右の式を±2/Tの範囲で定積分：x(t)x(t+τ)dt
end
yout(a) = s; % x(t)をτの分だけずらした時の積算値を格納
end
yout = yout / T;
% plot(xin(1:length(yout)), yout);
end
function [f, power] = disp_power_spectrum(sig)
fs = 1; % 帯域
y = fft(sig);
n = length(sig); % サンプル数
f = (0:n-1) * (fs/n); % 表示周波数
power = abs(y) .^ 2/n; % 離散フーリエ変換
% plot(f,power)
grid on;
end
参考にさせていただいた情報
大変お世話になりました。誠にありがとうございます。
[ネット]
[書籍]
システム同定の基礎（足立 修一 著）
最後に
ここはこうしたほうがいいぞ！ここ間違っているぞ！等あれば、お手数ですがご指摘いただけますと大変喜びます。