そのノッチフィルタ、オーバーシュートしますよ

この記事では、制御工学における伝達関数とオーバーシュートとの関係について触れます。間違い、ご指摘等ありましたらコメントいただけると助かります。

はじめに

先日、某SNSで↓の画像をあげたところ、そこそこ反響がありました。私自身も初めてこの挙動を見たときはけっこう驚きました。(以下、本記事では連続時間系を扱いますが、離散時間系でも同様の話です）

これを確認するためのMATLABコードは下記です。
※ Control System Toolboxが必要です(MATLAB Onlineなら無料で使えます)。

Notch_Overshoot_example.m

omega_n = 100*2*pi; %100[Hz]
%% 
% オーバーシュートしないノッチフィルタ

d1 = 0.1; % 深さパラメータ
zeta1 = 1; % 幅パラメータ

Notch1 = tf([1 d1*2*zeta1*omega_n, omega_n^2],[1 2*zeta1*omega_n, omega_n^2]);
text1 = append('(', poly2str(Notch1.Numerator{1},'s'), ')/(',  poly2str(Notch1.Denominator{1},'s'),')');
bodemag(Notch1,'b'); legend('オーバーシュートしないノッチフィルタ','Location','southeast');  grid on;
step(Notch1,'b',0.03); legend('オーバーシュートしないノッチフィルタ','Location','southeast');
pzmap(Notch1,'b'); legend;

%% 
% オーバーシュートするノッチフィルタ
% →幅パラメータζが小さい。

d2 = 0.1; % 深さパラメータ
zeta2 = 0.6; % 幅パラメータ

Notch2 = tf([1 d2*2*zeta2*omega_n, omega_n^2],[1 2*zeta2*omega_n, omega_n^2]);
bodemag(Notch2,'r'); legend('オーバーシュートするノッチフィルタ','Location','southeast'); grid on;
step(Notch2,'r',0.03); legend('オーバーシュートするノッチフィルタ','Location','southeast');

pzmap(Notch2,'r'); legend;

%% 
% 2つのノッチフィルタ比較

bode(Notch1,'b',Notch2,'r'); legend({'①オーバーシュートしない','②オーバーシュートする'},'Location','southeast'); grid on;
pzmap(Notch1,'b',Notch2,'r'); legend({'オーバーシュートしないノッチフィルタ','オーバーシュートするノッチフィルタ'},'Location','southwest');

%%
% ノッチの深さを変えたパターン。これもオーバーシュートしない

Notch1_2 = tf([1 0.01*2*zeta1*omega_n, omega_n^2],[1 2*zeta1*omega_n, omega_n^2]);
bode(Notch1_2,'g'); legend;
step(Notch1_2,'g',0.03); legend;
pzmap(Notch1_2,'g'); legend;

bodemag(Notch1,Notch2,Notch1_2,'g'); legend({'①オーバーシュートしない','②オーバーシュートする','①のノッチ深い版。オーバーシュートしない'},'Location','southeast'); grid on;

特に重要な、青(オーバーシュートしない)と赤(オーバーシュートする)のボード線図とステップ応答について再度載せておきます。

さらに、上記赤線のノッチフィルタは、(SNSに載せる都合上)ある程度はっきりオーバーシュートがわかるものを選びました。パラメータをギリギリオーバーシュートする…程度にまで変える(zeta→0.9)と下図のピンクのようになります。

青とピンクの周波数特性にもはや違いは見受けられません(SNS上では位相特性に関係が…という意見がありましたが、これらを見るとそうとも言えなさそうな気がしています)。以下では、ボード線図上でほぼわからないこれらの特性が、なぜステップ応答のオーバーシュートという違いに表れるのかを確認していきます。

ただし、一番気になる「任意の伝達関数(やその極零配置)を見たとき、ステップ応答がオーバーシュートするか判定する方法は？」の一般論を言うのはなかなか難しそうです。これらの議論は概要だけ述べて後日(できたら)に回すこととします。

オーバーシュートについて

のちのちのために、本記事でいうオーバーシュートについて明確にしておきます。

本記事では、↓の図のようにオーバーシュートを定義することとします。

つまり率直に、

• 制御系の出力が一瞬でも目標値を超えたらオーバーシュートしている
• そうでなかったらオーバーシュートはしていない

という意味とします。

図中の上2つをオーバーシュート、と定義している文献は多数ありますが、下2つのような波形をオーバーシュートではない、と明確にしている文献は少ないようで、ひょっとしたらこの解釈は特殊なものかもしれないのでご注意ください。

なお、以降も基本的に伝達関数のステップ応答を考えていくので、上図でいう目標値はステップ応答の定常値のこと、ととらえてもらえればと思います。

ノッチフィルタについて

下記の記事で大変わかりやすくまとまっています。

本記事に関わる内容を再掲すると、ノッチフィルタの連続時間系での伝達関数は下記となります。ノッチ周波数$\omega_n$、ノッチの幅パラメータ$\zeta$、深さパラメータ$d$で構成されています。
$$G(s) = \frac{s^2+d\cdot2\zeta\omega_ns+\omega_n^2}{s^2+2\zeta\omega_ns+\omega_n^2} $$

ステップ応答との関係についても触れておきます。ノッチフィルタは「ローパスフィルタとハイパスフィルタの組み合わせ」ととらえることができ、その周波数特性（ゲイン特性）とステップ応答の振る舞いの関係は、ざっくり下図のように解釈できます。ローパスフィルタのみだと時定数分の時間遅れだけですが、ハイパスフィルタと組み合わさることで初期値が1、という挙動があらわれていると考えることができます。

振動抑制のためのノッチフィルタ

本記事でなぜそれほどオーバーシュートに着目するのか、について述べたいと思います。

ノッチフィルタのよくある使いどころとして、機械系の振動抑制があります。特に、2慣性系といわれる制御対象を考えるとその効果・課題がわかりやすいので紹介していきます。検証に使ったコードはこの話題が終わったあたりで置いています。

2慣性系の特性

文献¹を参考に、下図のようにモータと回転負荷、間にシャフトを通してトルクが伝達される系の回転速度を制御することを考えます。モータトルク$T_e$ が印加されたときモータが$\omega_m$で回転し、シャフトで繋がった負荷は$\omega_L$で回転する、という様子をモデル化したものです。
(もちろん位置制御を考えても以下の話は同様ですが、話を単純化するために速度制御を考えます)

本記事の問題設定では、モータ回転速度$\omega_m$は計測できる一方、負荷側の速度$\omega_L$が計測できないものとします。しかし本当に精度よく動かしたいのは負荷側…というような状況だとします。

導出過程は省略して結果だけ載せると、モータトルク→モータ角速度の特性$G_{MM}(s)$、モータトルク→負荷の角速度の特性$G_{ML}(s)$はそれぞれ下記で得られます。

\begin{align} 
G_{MM}(s) = 
\frac{\omega_m(s)}{T_e(s)} = \frac{1}{Js}\frac{\left(\frac{1}{\omega_z}\right)^2
 s^2 +2 \left(\frac{\zeta_z }{ \omega_z}\right)s + 1}{ \left(\frac{1}{\omega_a}\right)^2
 s^2 +2 \left(\frac{\zeta_a }{\omega_a}\right)s + 1}, \\
G_{ML}(s) = 
\frac{\omega_L(s)}{T_e(s)} = \frac{1}{Js}\frac{
 2 \left(\frac{\zeta_z }{ \omega_z}\right)s + 1}{ \left(\frac{1}{\omega_a}\right)^2
 s^2 +2 \left(\frac{\zeta_a }{\omega_a}\right)s + 1}
\end{align}

ここで、上式の各定数は↓で定義しています。

\begin{align} 
J  = J_M + J_L \ , \\
\omega_a = \sqrt{K_F \left(\frac{1}{J_M} + \frac{1}{J_L}
 \right)} \ , \\
 \omega_z = \frac{\omega_a}{\sqrt{1+ J_L/J_M}}  \ , \\
 \zeta_a = \frac{C_F (1+ J_L/J_M)}{ 2 \omega_a J_L} \ , \\
\zeta_z = \frac{\zeta_a}{ \sqrt{1 + J_L/J_M}}
\end{align}

適当な定数を設定すると、ボード線図は下図となります。共振のピーク状の特性が$\omega_a$で、反共振といわれるノッチ状の特性が$\omega_z$であらわれていることがわかります。

このような2慣性系の特性は機械系ではよくあらわれます。特に「モータ＋ボールねじ＋直線移動する負荷」という系も2慣性系とみなすことができ、詳しくは文献²などを参照してください。

ここで$G_{MM}(s)$と$G_{ML}(s)$との関係性に着目すると、両者の間には下図のような特性が隠れていることがわかります。

$$ W_{L}(s) = \frac{G_{ML}(s)}{G_{MM}(s)}= \frac{
2 \left(\frac{\zeta_z }{ \omega_z}\right)s + 1}{\left(\frac{1}{\omega_z}\right)^2
s^2 +2 \left(\frac{\zeta_z }{ \omega_z}\right)s + 1} $$

これらの関係性を図でまとめると下図です。各伝達関数のゲイン線図のイメージ図を載せています。

この図に対して、適当なPI制御器$C_v(s)$を配置して速度制御を行います。すると下図となります。

この速度制御系の閉ループのステップ応答を見ると、下図となります。それぞれ、モータの角速度と負荷の角速度です。

適当に作った制御系なので、やはりモータ側、負荷側どちらの回転速度も振動を伴っています。実は、これら振動周波数はどちらもおよそ100Hzで、ほぼ「反」共振周波数$\omega_z$となっています。

なぜ２慣性系の共振周波数$\omega_a$ではなく「反」共振周波数$\omega_z$で振動するのか、その理由をもう少し詳しく見ましょう。先ほど速度制御器$C_v(s)$を挿入して得た閉ループ系(回転速度指令→モータ角速度)の特性を$G_{vm}(s)$、さらにその先の負荷までの系(回転速度指令→負荷角速度)の特性を$G_{vl}(s)=W_L(s)G_{vm}(s)$とおいて、それぞれボード線図で図示すると下記となります。

各々の関係性をブロック図で示すと、下図です。

大事なのは、閉ループを組んだことでもともとの$G_{MM}(s)$の共振周波数$\omega_a$のピークは抑制されている一方、反共振周波数でのノッチが浅くなっています。これによって、打ち消せていた$W_L(s)$のピーク特性が打ち消せなくなり、全体的には、もともと「反」共振周波数だった帯域でピークを持つようになります。

これによって、回転速度指令→負荷角速度の特性$G_{vl}(s)=W_L(s)G_{vm}(s)$ではステップ応答がオーバーシュートし、「反」共振周波数で振動していたのです。

もう一点大事なポイントとして、回転速度指令→モータ角速度の特性$G_{vm}(s)$の方はノッチがあるだけですが、こちらも「反」共振周波数で振動しつつオーバーシュートしています。こうしたケースの解釈は話が長くなりそうなので、やはり後日に回そうと思います。

ノッチフィルタによる振動抑制

上述の特徴を持つ２慣性系に対して振動抑制を考えた場合、

• フィードバック制御器を作り込む
• フィードフォワード制御器で対処する

の2通りやり方がありますが、後者の代表格が「ノッチフィルタを配置する」です。

(余談ですが…)

文献によっては、フィードフォワードだけでなく、フィードバック制御器にノッチフィルタを配置して振動抑制を試みる内容もあります。また、ステップ状ではなく指令カーブを工夫して振動抑制する、ということもよく行われます。指令カーブの工夫…は突き詰めると、やはりそれもフィードフォワード制御器を設計しているのと等価です。このあたりのメリットデメリットは文献³が詳しいです。文献³では指令カーブの整形手段としてプリシェイピング法と呼ばれる手法を紹介していて、ノッチフィルタとの違いなどがまとまっています。

下図のように、フィードフォワード的にノッチフィルタを指令の個所へ挿入します。フィードフォワード的に挿入したことで、指令波形を整形してフィードバック制御系に与える、という形になります。このように「おおよそフィードバック制御器を調整した後、最後に振動抑制のためにフィルタを調整する」というステップは現場で広く行われているのでは、と想像します。フィードバックループの内部を変更しないため発振・不安定化のリスクが少なく、調整が比較的容易です。

ただしこのフィードフォワードのノッチフィルタの選定について、主題の「ノッチフィルタのオーバーシュートする/しない」が問題となります。ノッチフィルタを挿入すると時間遅れが発生し、ノッチを深く・幅が太くしていくとより遅れは大きくなります。それを嫌ってノッチを細くするとオーバシュートする指令が生成されることがあるのです。

この数値例は系の振動周波数、すなわち２慣性系の 「反」共振周波数が100[Hz]付近なので、冒頭の(ノッチが100Hzの)フィルタを試してみます。ノッチの深さは変更しないことにして、冒頭の赤/青のオーバーシュートする/しないの２通りを試します。

先ほどの系に、オーバーシュートする/しないノッチフィルタを挿入して得られた回転速度のステップ応答を示すのが下図です。

オーバーシュートするノッチフィルタを挿入した場合、確かにフィルタなしに比べて振動は抑えられるものの、やはりオーバーシュートの大きさが気になります。特に、モータ角速度の方はフィルタなしに比べてもオーバーシュートの最大値が増えるようです(これは数値例で変わります)。冒頭で同じ波形を見せたので省略しますが、この違いは「整形された回転速度指令」の段階でオーバーシュートする波形になっているためです。

また、こうした「指令波形がオーバーシュートしている」状況はソフト実装上いろいろな不具合の原因になります。通常、目標指令の生成は完全にソフトウェアで完結した話なので、時々刻々の目標値が、最終的な目標値を超えることはないだろう…と思い込みがちです。特に位置指令を考えた場合、こうしたオーバーシュートする位置指令を使ってもろもろの判定・処理をすると意図しない挙動になり得るため、指令のオーバーシュートには注意が必要です。

長くなりましたが、本記事でノッチフィルタのオーバーシュートを問題視するのはこうした背景からです。もちろん、「ある周波数成分を取り除きたい」という用途でノッチフィルタを使う場合、こうしたオーバーシュートの話はあまり気にしなくて良いと思われます。

ここで使ったMATLABコードはこちら

%% 
% 2慣性系の振動抑制にノッチフィルタを使う、の例

clear;
s = tf('s');

omega_n = 100*2*pi; % [rad/s]
zeta1 = 1; % ノッチの幅 1未満だと、オーバーシュートする
d1 = 0.1; % ノッチの深さ

% オーバーシュートしないノッチフィルタ
Notch1 = tf([1 2*zeta1*omega_n*d1 omega_n^2],[1 2*zeta1*omega_n omega_n^2]);

zeta2 = 0.5;
d2 = 0.1;

% オーバーシュートするノッチフィルタ
Notch2 = tf([1 2*zeta2*omega_n*d2 omega_n^2],[1 2*zeta2*omega_n omega_n^2]);
bode(Notch1,Notch2);


%% 
% 2慣性系の記述 (ACモータ可変速制御システムの理論と設計)を参考

ratioJ = 1;% JL/JM
JM = 0.00255;% [kg m^2]
JL = JM*ratioJ; % [kg m^2]

J = JL+JM;
KF = 1000; %[K g cm/rad] 
CF = 0.1000; %[K g cm]

omega_a = sqrt(KF*(1/JM + 1/JL));
omega_z = omega_a / sqrt(1+JL/JM);
zeta_a = CF*(1+JL/JM) / (2*omega_a*JL);
zeta_z = zeta_a / sqrt(1+JL/JM);

omega_a_Hz = omega_a/(2*pi);
omega_z_Hz = omega_z/(2*pi);

%% 
% 制御対象の伝達関数

G_MM = tf([(1/omega_z^2) 2*zeta_z/omega_z 1], [(1/omega_a^2) 2*zeta_a/omega_a 1])/s /J ; 
G_ML = tf([2*zeta_z/omega_z 1], [(1/omega_a^2) 2*zeta_a/omega_a 1])/s /J ; 

% 周波数特性の描画
bode(G_MM,G_ML);
xlim([10 1000]);
grid on;
legend({'G_{MM}(s)','G_{ML}(s)'})
text(100,-40,['\omega_z : ' num2str(omega_z_Hz) ' [Hz]'])
text(100,-50,['\omega_a : ' num2str(omega_a_Hz) ' [Hz]'])
% 速度ループの解析
% 速度PI制御器Cv_pi

Ksp = 20;
Ksi = 500;

Cv_pi = tf([Ksp Ksi],[1 0]);

Gv_motor = feedback(G_MM*Cv_pi,1);
opts = bodeoptions('cstprefs');
opts.PhaseWrapping = 'on'; %位相線図のY範囲を[-180°,180°]にする
opts.Title.String = 'モータ角速度指令→モータ角速度 のボード線図';
bodeplot(Gv_motor,opts);
xlim([10 1000]);
grid on;
legend({'G_v(s)'})

%% 
% モータ角速度→機械端角速度の間に存在する特性W_Lを使って、
% 
% モータ角速度指令→機械端角速度 の閉ループ系の伝達関数Gv_loadを求める。

% モータ角速度→機械端角速度 の伝達関数W_L
W_L = tf([2*zeta_z/omega_z 1], [(1/omega_z^2) 2*zeta_z/omega_z 1]); 
opts.Title.String = 'モータ角速度→負荷角速度 のボード線図';
bodeplot(W_L,opts);

Gv_load = W_L*Gv_motor;
opts.Title.String = 'モータ角速度指令→負荷角速度 のボード線図';
bodeplot(Gv_load,opts);

% まとめて描画する場合
% bode(Gv_motor,W_L,Gv_load)
% xlim([10 1000]);
% grid on;
% legend({'G_{vm}(s)','W_L(s)','G_{vl}(s)'})
%% 
% 

step(Gv_motor,0.1);
title('モータ角速度にステップ指令を与えたときのモータ角速度の応答')
grid on;
step(Gv_load,0.1);
title('モータ角速度にステップ指令を与えたときの負荷角速度の応答')
grid on;
%% 
% ノッチフィルタによる振動抑制効果確認

step(Gv_motor,Gv_motor*Notch1,Gv_motor*Notch2,0.1);
title('モータ角速度にステップ指令を与えたときのモータ角速度の応答')
legend({'ノッチフィルタなし','オーバーシュートしないノッチフィルタの適用','オーバーシュートするノッチフィルタの適用'},'Location',"southeast")

step(Gv_load,Gv_load*Notch1,Gv_load*Notch2,0.1);
title('モータ角速度にステップ指令を与えたときの負荷角速度の応答')
legend({'ノッチフィルタなし','オーバーシュートしないノッチフィルタの適用','オーバーシュートするノッチフィルタの適用'},'Location',"southeast")

オーバーシュートする/しないノッチフィルタの違い

だいぶ横道に入りましたが、本題に戻ります。結論から言うと、冒頭の画像のとおり極の位置が違います。

ボード線図上ではほぼ違いがわからないピンクのものも含めて、それぞれの伝達関数の極・零点配置を示すと、下図となります(pzmapの色を変えるのが大変そうだったので、マゼンタ→ピンクだと思って見てください)。

ゼロ点の位置がほぼ変わらないのに対し、極の位置はオーバーシュートする/しないで明確に異なります。オーバーシュートするものは、極が複素成分を持っていることがわかります。一方、オーバーシュートしないものは実軸上に(画像だと見えませんが)2つ重なった極となっています。2次遅れ系の応答でよく知られている通り、極の複素成分を持つ系は振動を引き起こします。

...なるほど、極の複素成分があるかないかで決まるのかあ…完！ と終わってもいいんですが、まだもう少し掘り下げます。

オーバーシュートする/しないノッチフィルタのボード線図上での違い

なぜこれらの特性がボード線上には見えてこないのかを確認してみます。

上記の特性は、ボード線図上では隠れてしまっています。極零が影響するとわかったので、ノッチフィルタの分母分子で伝達関数を分けると下記となります。

$$G(s) = \frac{s^2+d\cdot2\zeta\omega_ns+\omega_n^2}{s^2+2\zeta\omega_ns+\omega_n^2} =: G_{\rm num}(s)G_{\rm den}(s)$$

ここで、$$G_{\rm num}(s) = \frac{s^2+2(d\zeta)\omega_ns+\omega_n^2}{\omega_n^2}, $$

$$G_{\rm den}(s) = \frac{\omega_n^2}{s^2+2\zeta\omega_ns+\omega_n^2} $$
と置きました。分母分子の$\omega_n^2$は定常ゲインを1にそろえるためのものです。また、違いを分かりやすくするために$d \zeta$をかっこでくくりました。

この分解をしたうえで、冒頭のボード線図を再度描画すると下図となります。ボード線図は対数グラフなので、掛け算の形で分解した伝達関数はグラフとしては足し合わせで考えることができます。

…ちょっと違いがわかりにくいので、$\zeta=0.05$とした、より極端にオーバーシュートする例も見てみます。

この図で大事なのは、分母(den)の周波数特性です。オーバーシュートしない例ではゲインピークがあらわれていませんでしたが、ノッチの周波数帯でピークを持つ特性が隠れている、と言えます。

これは元の伝達関数の式を眺めてみると明らかで、ノッチフィルタは

• 分母(den)の二次遅れ系$G_{\rm den}$
• 分子(num)の二次遅れ系(ただし減衰係数が$\zeta$ではなく、$d \zeta$ )の逆特性$G_{\rm num}$

から構成されています。通常、$d<1$のため$d \zeta < \zeta$となる点も注意ください。

上のボード線図をみると、より細いノッチを表現するために分子(num)だけでは足りず、分子(num)と分母(den)の波形の足し合わせが必要な様子がわかります。このより細いノッチを表すための分母のピーク特性が、複素成分を含む極として表れてきている、と解釈できそうです。分母の２次遅れ系はよく知られている通り、$\zeta \geq 1$なら周波数特性がピークを持たず、ステップ応答はオーバーシュートしません。(実際、先の赤の②オーバシュートする例は$\zeta < 1$で、ボード線図を拡大するとdenのゲインが少しだけ持ち上がっています)

これはなかなか興味深い(はず)です。細い幅で急峻なノッチを非プロパーな2次のゼロ点だけで表現するのには限界があり、極を複素極にしてゼロ点と近づけることで、初めて実現できる…という風に見えます。 理想的な遮断特性のフィルタは実現できず現実にはリプルが生じる…などの性質とも関連してそうに思えますが、あまり深いことは言えなさそうなので詳しい方のツッコミを待つことにします。

オーバーシュートする/しないノッチフィルタのステップ応答の逆ラプラス変換

さらに別角度から理解を深めるために、ノッチフィルタのステップ応答を逆ラプラス変換から求めてみましょう。
MATLAB Symbolic Math Toolboxを使います(これもMATLAB Onlineなら無料で使えます)。

syms s t omega_n zeta d

% ノッチフィルタの伝達関数
G = (s^2 + d * 2 * zeta * omega_n * s + omega_n^2) / (s^2 + 2 * zeta * omega_n * s + omega_n^2);
% ノッチフィルタのステップ入力印加時の伝達関数
G_u = G/s;

% 逆ラプラス変換によるステップ応答
g_t = ilaplace(G_u, s, t);
disp(g_t);

表示された式を整理すると下記になります。

\begin{align} 
g(t) = 1-\frac{2\zeta \left(1-d\right)}{\sqrt{1-\zeta^2} }e^{- \zeta \omega_n t } \sin \left( \sqrt{1-\zeta^2 } \omega_n t\right) 
\end{align}

時間応答を確認する場合は下記のコードを活用ください。

zeta_val = 0.8; % 1ぴったりは分母ゼロになるので不可。1未満で、オーバーシュートするパターン。
% zeta_val = 1.1; % 1ぴったりは分母ゼロになるので不可。1以上で、オーバーシュートしないパターン。

omega_n_val = 100 * 2 * pi;
d_val = 0.1;

% シンボリック関数の数値化
g_sim = subs(g_t, [zeta, omega_n, d], [zeta_val, omega_n_val, d_val]);
disp(g_sim);

% 数値関数に変換
g_sim_func = matlabFunction(g_sim);

% 時間範囲の設定
% t_values = linspace(0, 0.001, 10); 
t_values = 0:0.0001:1; 

% g_sim(t)の値を計算
g_values = g_sim_func(t_values);

% グラフ描画
figure;
plot(t_values, g_values);
xlabel('Time (s)');
ylabel('g(t)');
title('Time Response of the Notch Filter');
grid on;

ちなみに、このノッチフィルタのステップ応答は手計算でも求められます。驚くほどスッキリ分解できたため、せっかくなので載せておきます。

ノッチフィルタのステップ応答は、伝達関数$G(s)$に$1/s$をかけて逆ラプラス変換することで得られます。

\begin{align} 
\frac{G(s)}{s} &= \frac{1}{s}\frac{s^2 + 2 d \zeta \omega_n s+ \omega_n^2}{s^2 + 2 \zeta \omega_n + \omega_n^2} \\
&= \frac{1}{s}\left(1-\frac{ 2 \zeta \omega_n (1-d)s }{s^2 + 2 \zeta \omega_n + \omega_n^2} \right) \\
&= \frac{1}{s}-\frac{ 2 \zeta \omega_n (1-d) }{s^2 + 2 \zeta \omega_n + \omega_n^2} \\
&= \frac{1}{s} - \frac{ 2 \zeta \omega_n (1-d) }{(s+\zeta \omega_n)^2 + \omega_n^2 (1-\zeta^2)} \\
&=  \frac{1}{s} - \frac{2 \zeta (1-d) }{\sqrt{1-\zeta^2}}\frac{ \omega_d}{(s+\zeta \omega_n)^2 + \omega_d^2 }
\end{align}

最後の式で、$\omega_d^2 = \omega_n^2 (1-\zeta^2)$ と置き換えました。

逆ラプラス変換$\mathcal{L}^{-1} \left[ \frac{1}{s}\right] = 1$、$\mathcal{L}^{-1} \left[ \frac{\omega}{(s+a)^2 + \omega^2}\right] = e^{-a t} \sin \omega t$ と、$\omega_d$の置き換えを戻すことで下記の時間応答が得られます。

\begin{align} 
g(t) = 1-\frac{2\zeta \left(1-d\right)}{\sqrt{1-\zeta^2} }e^{- \zeta \omega_n t } \sin \left( \sqrt{1-\zeta^2 } \omega_n t\right) 
\end{align}

この応答の式をみると、やはり$\sqrt{1-\zeta^2}$の部分で応答が大きく変わります。$0<\zeta<1$なら、三角関数と指数減衰の組み合わせで振動し、1を超えることがありそうです。$\zeta>1$なら$\sqrt{1-\zeta^2}$が複素数になり、少々変わった波形になりそうです。一方、$d$の方は指数関数・三角関数の係数にだけ効いてきて、$0<d<1$ (0に近い方がノッチが深い)なので波形はあまり大きく変わらないことがわかります。

ノッチフィルタで減衰係数ζ>1となる場合の時間応答

$\sqrt{1-\zeta^2}$が複素数になる場合について、上述の式を$j \alpha$という純虚数で置き換えてみます。$\alpha$自体は実数で、$\alpha=\sqrt{\zeta^2-1}>0$となります。

\begin{align} 
g(t) &= 1-\frac{2\zeta \left(1-d\right)}{j \alpha}e^{- \zeta \omega_n t } \sin \left( j \alpha \  \omega_n t\right)  \\
&= 1-\frac{2\zeta \left(1-d\right)}{j \alpha}e^{- \zeta \omega_n t } \frac{ e^{ -\alpha \omega_n t} - e^{\alpha \omega_n t}}{2j} \\
&= 1+\frac{\zeta \left(1-d\right)}{\alpha}e^{- \zeta \omega_n t } \left( e^{ -\alpha \omega_n t} - e^{ \alpha \omega_n t} \right) \\
&= 1+\frac{\zeta \left(1-d\right)}{\alpha} \left(   e^{-(\alpha+\zeta)\omega_nt } -e^{(\alpha-\zeta)\omega_nt}  \right)
\end{align}

途中、複素数がキレイに消えて指数関数だけになるのが面白いですね。$\alpha$の置き換えを戻すと下記となります。

\begin{align} 
g(t) &=  1+\frac{\zeta \left(1-d\right)}{\sqrt{\zeta^2 -1 }} \left(  e^{-(\sqrt{\zeta^2 -1 }+\zeta)\omega_nt } - e^{(\sqrt{\zeta^2 -1 }-\zeta)\omega_nt} \right)
\end{align}

3項目の指数関数の指数部は、$\sqrt{\zeta^2-1}<\zeta$なので指数部が負の指数関数となります。つまり、２項目も３項目も0へ減少する指数関数で、２項目の方が0への収束が早くなります。振る舞いを確認するため簡単なMATLABコードを実行すると、下記のようになります。

omega_n = 2*pi*100;
zeta = 1.2;
t = 0:0.0001:0.05; 

SecondTerm = exp( -(sqrt(zeta^2-1)+zeta)*omega_n*t );
ThirdTerm = exp( (sqrt(zeta^2-1)-zeta)*omega_n*t );

plot(t,SecondTerm,t,ThirdTerm,t,SecondTerm-ThirdTerm)
legend({'$e^{-(\sqrt{\zeta^2 -1 }+\zeta)\omega_nt}$','$e^{(\sqrt{\zeta^2 -1 }-\zeta)\omega_nt}$','$e^{-(\sqrt{\zeta^2 -1 }+\zeta)\omega_nt} - e^{(\sqrt{\zeta^2 -1 }-\zeta)\omega_nt}$'},'Interpreter','latex')

これを見ると、ノッチフィルタのステップ応答の特徴的な挙動があらわれています。↑の黄色の波形は大小関係的に、0を上回ることも-1を下回ることもありません。指数関数にかかっている係数は0以上1未満なので、１＋↑の黄色の波形を0～1の範囲で拡大縮小したもの、がノッチフィルタのステップ応答となっていたわけです。

ちなみに、途中で双曲線関数sinhで置き換えるともう少しスッキリした式になります。こちらの方が解釈しやすい…という方はご利用ください。

\begin{align} 
g(t) = 1-\frac{2\zeta \left(1-d\right)}{ \sqrt{\zeta^2-1} }e^{-\zeta \omega_n t } \sinh( \sqrt{\zeta^2-1} \ \omega_n t) 
\end{align}

オーバシュートする/しない伝達関数の見極め方

結論から言うと、より一般的な条件で見極める方法はないと思われます。(あったら教えてほしい…)

これまで述べてきた通り、「ノッチフィルタの極に複素成分が含まれている場合はオーバーシュートする」ということがわかりました。ただし、これを拡大解釈して「極に複素成分がある伝達関数は、ステップ応答がオーバーシュートする」という結論にはなりません。例をあげてみましょう。

下記のように、$-2, -3\pm 4 i$の3つの極を持つ伝達関数を考えます。このステップ応答は、オーバーシュートしません。

Q = zpk([],[-2 -3+4i -3-4i],1);
pzmap(Q,'g')
step(Q,'g',5)

この伝達関数はゼロ点がなく、複素成分の影響がキャンセルされている…というわけでもありません。(慣れている方は、虚軸に近い方の極が支配的になるから…と推察できそうではありますね）

こうした伝達関数に、(たとえ安定ゼロに限ったとしても)ゼロ点が加わってくると解釈はより難しくなりそうです。

インパルス応答が負になること、に着目したオーバーシュートの判定条件

オーバーシュート有無の判定条件に関しては多数文献があり、例えばこれら⁴⁵などがあります。そのいくつか(多く？)は「インパルス応答が負になる → ステップ応答が極大値を持つ」という性質に着目しています。

例えば先の3つの極を持つ伝達関数$Q$について、インパルス応答を求めると下記になります。インパルス応答とステップ応答との関係性から明らかに、インパルス応答は「ステップ応答の時間変化率」を表しています。したがって、インパルス応答の符号が変化しない場合はステップ応答の波形は単調増加となり、極大値を持たない、と言えます。

しかし、「極大値を持つ」のと「オーバーシュートする」は似て非なるもので、「オーバーシュートする伝達関数は必ずインパルス応答が負になっているか」というとそうではありません。

例えば、先のオーバーシュートしないノッチフィルタの伝達関数について、インパルス応答をプロットすると下記です。これはステップ応答の時間変化を振り返ると明らかで、初期値が1の後急激に一度0へ向かい、そこから再度1へ向かうという挙動からくるものです。

したがってインパルス応答が負になること、に着目した場合、先に述べた「オーバーシュートと呼ばない」ものも弾いてしまうことになり、(本記事にとっては)保守的な条件となります。もちろん、先の文献に誇張表現があるわけではありません。ただ、このように「オーバーシュートするか否か」と「ステップ応答が極大値をとるか否か」を明確に分けている文献は少ない印象です。

これらが難しいポイントで、より一般的に伝達関数からオーバーシュートを見分ける術はどうしたらいいのか…の試みについては後日(できたら)に回したいと思います。

おわりに

本記事では、ノッチフィルタのステップ応答がオーバシュートする/しないがどこで変わってくるのかを述べました。今回扱ったノッチフィルタは2次・不安定ゼロを持たないという比較的身近な存在な気がしますが、その時間応答はけっこう不思議な点がいっぱいありますね。(古典制御をもっとまじめに学んでいれば…)

後日もう少しだけ一般的なことが言えればと思いますが、定量的なことは言えず、試行錯誤する見込みです。これらに関わる文献をご存じの方はぜひぜひ教えてください…。

1. 「ACモータ可変速制御システムの理論と設計」, 杉本英彦, 小山正人, 玉井伸三, 森北出版, https://www.morikita.co.jp/books/mid/074401 ↩

2. 「精密位置決め・送り系設計のための制御工学」, 松原 厚, 森北出版, https://www.morikita.co.jp/books/mid/091981 ↩

3. 「よくわかる機械の制振設計 防振メカニズムとフィードフォワード制御による対策法」, 三好孝典, 日刊工業新聞社, https://pub.nikkan.co.jp/book/b10021429.html ↩ ↩²

4. "Some conditions on zeros to avoid step-response extrema,", A. Rachid, in IEEE Transactions on Automatic Control, vol. 40, no. 8, pp. 1501-1503, Aug. 1995, https://ieeexplore.ieee.org/document/402253/ ↩

5. "Influence of Zero Locations on the Number of Extrema in the Step Response," M. El-Khoury and R. Longchamp, American Control Conference, Boston, MA, USA, 1991, pp. 1192-1194, https://ieeexplore.ieee.org/document/4791567 ↩