線形制御系設計

    線形システムの表現
    線形システムの特性解析
    制御器設計
    観測器設計
    離散時間システム

    モーションコントロール

    加速度制御
    剛体の姿勢制御

    モータドライブ

    駆動制御法
    磁極位置推定法

    システム同定

    ホワイトボックスモデル

伝達関数の離散化

制御器やフィルタに所望の周波数特性を持たせることを検討した場合,その設計は周波数領域行われる。 デジタルプロセッサを用いてこれらを実装するためには,離散的にサンプルされたデータを用いて所望の周波数特性を実現する入出力写像を求める必要がある。

問題定義

以下では所望の周波数特性を有する系を ff, この系に対する連続時間入力信号および出力信号を x, yx,\ y とする。ただし,x, yx,\ y は時刻 t0t\ge0において有限の値を持つものとする。サンプラを有する離散システムにおいて利用可能な状態量は,現時刻のサンプル値と過去のサンプル値群である。サンプル周期を TsT_{\rm s}とし,最新のサンプル値が時刻 tk, θ=kTs+θ (k0, 0<θ<Ts)t_{k,\ \theta}=kT_{\rm s} + \theta\ (k\ge0,\ 0<\theta < T_{\rm s}) にて取得されたとすると,利用可能なデータの集合は以下のように記述される。
X(k, θ)={x((kn)Ts+θ) n[0, k]}Y(k, θ)={y((kn)Ts+θ) n[1, k]}\begin{align} X(k,\ \theta)&=\left\{x\left((k - n)T_{\rm s} + \theta\right) |\ n \in [0,\ k]\right\}\\ Y(k,\ \theta)&=\left\{y\left((k - n)T_{\rm s} + \theta\right) |\ n \in [1,\ k]\right\} \end{align}
ここで,連続系の写像 ff を離散系で実現する場合の写像を fdf_{\rm d} とすると,この写像は次の形式で記述される。
fd: (X(k, θ), Y(k, θ)) y(tk, θ)\begin{align} f_{\rm d}:\ (X(k,\ \theta),\ Y(k,\ \theta)) \mapsto \ y(t_{k,\ \theta}) \end{align}

時間シフトを用いたシステムの実現

システム fdf_{\rm d} の実装について,次の線形システムを実現することを考える。
y(t)=k=0pakx(tkTs)k=1qbky(tkTs)\begin{align} y(t) &= \sum_{\rm k=0}^{\rm p} a_{\rm k}x(t - kT_{\rm s}) - \sum_{\rm k=1}^{\rm q} b_{\rm k}y(t - kT_{\rm s}) \end{align}
ただし, ai,bj (i[0, p], j[1, q])a_{\rm i}, b_{\rm j}\ ({\rm i}\in[0,\ {\rm p}],\ {\rm j}\in[1,\ {\rm q}]) は定数とする。このシステムでは,任意の時刻における出力が一定時刻前の入出力信号から生成される。すなわち,サンプラが導入されるシステムにおいてもサンプル値の関係性がそのまま説明可能であり,x, yx,\ y のサンプル値について次の関係が成立する。
y[n]=k=0pakx[nk]+k=1qbky[nk]\begin{align} y[{\rm n}] &= \sum_{\rm k=0}^{\rm p} a_{\rm k}x[{\rm n} - k] + \sum_{\rm k=1}^{\rm q} b_{\rm k}y[{\rm n} - k] \end{align}
ここで,システムの入出力をラプラス変換することにより,次を得る。
Y(s)=k=0pakeskTsX(s)+k=1qbkeskTsY(s)(1+k=1qbkeskTs)Y(s)=k=0pakeskTsX(s)Y(s)=k=0pakeskTs1+k=1qbkeskTsX(s)\begin{align} Y(s) &= \sum_{\rm k=0}^{\rm p} a_{\rm k}e^{-s{\rm k}T_{\rm s}}X(s) + \sum_{\rm k=1}^{\rm q} b_{\rm k}e^{-s{\rm k}T_{\rm s}}Y(s) \\ \Leftrightarrow (1 + \sum_{\rm k=1}^{\rm q} b_{\rm k}e^{-s{\rm k}T_{\rm s}})Y(s) &= \sum_{\rm k=0}^{\rm p} a_{\rm k}e^{-s{\rm k}T_{\rm s}}X(s) \\ \therefore Y(s) &= \frac{\sum_{\rm k=0}^{\rm p} a_{\rm k}e^{-s{\rm k}T_{\rm s}}}{1 + \sum_{\rm k=1}^{\rm q} b_{\rm k}e^{-s{\rm k}T_{\rm s}}}X(s) \end{align}
すなわち,システム fdf_{\rm d} の伝達関数 GdG_{\rm d}は次のように記述される。
Gd(s)=k=0pakeskTs1+k=1qbkeskTs\begin{align} G_{\rm d}(s) &= \frac{\sum_{\rm k=0}^{\rm p} a_{\rm k}e^{-s{\rm k}T_{\rm s}}}{1 + \sum_{\rm k=1}^{\rm q} b_{\rm k}e^{-s{\rm k}T_{\rm s}}} \end{align}
離散時間系の実装問題は,GdG_{\rm d} が所望の伝達関数を持つように係数 ai,bja_{\rm i}, b_{\rm j} を決定することに帰着する。

シフトオペレータの導入

サンプル値の関係性を簡潔に記述するために,シフトオペレータ zz を導入する。
zesTs\begin{align} z &\equiv e^{\rm sT_{\rm s}} \end{align}
これにより,GdG_{\rm d}zz の関数として記述することができる。
Gd(z)=k=0pakzk1+k=1qbkzk\begin{align} G_{\rm d}(z) &= \frac{\sum_{\rm k=0}^{\rm p} a_{\rm k}z^{-{\rm k}}}{1 + \sum_{\rm k=1}^{\rm q} b_{\rm k}z^{-{\rm k}}} \end{align}

離散化手法

以下ではシステム ff の伝達関数を GG とする。一般に GG は周波数領域で設計された ss の関数であるが,これを zz の関数に変換することができればサンプル値を用いて実装可能となる。これには複数の方法があり,ここでは一部のみを紹介する。

後退差分法

シフトオペレータの定義より,次の関係が成立する。
s=1Tsln(z)\begin{align} s = \frac{1}{T_{\rm s}}\ln(z) \end{align}
右辺について,z=1z=1 を中心としたテイラー近似多項式を求めることを考える。この中心点の選択については,s=0s=0 近傍(低周波数帯)での近似精度を確保するために必要となる。まず,自然対数関数の n\rm n 次導関数は次のように記述される。
dndznln(z)=(1)n1(n1)!xn\begin{align} \frac{d^{\rm n}}{dz^{\rm n}}\ln(z) = (-1)^{{\rm n} -1}\frac{({\rm n} - 1)!}{x^{\rm n}} \end{align}
したがって,ss の近似多項式は次のように記述される。
s1Tsn=0(1)n11n(z1)n1\begin{align} s \approx \frac{1}{T_{\rm s}}\sum_{\rm n=0}^{\infty} (-1)^{{\rm n} -1}\frac{1}{{\rm n}}(z-1)^{{\rm n} - 1} \end{align}
この近似を1次で打ち切ることにより,ss を前進差分によって近似した式が現れる。
sz1Ts\begin{align} s \approx \frac{z-1}{T_{\rm s}} \end{align}
ここで,制御器実装に際しては未来値を使用することはできないため,zz による位相遅れを許容して後退差分とする。
sz1z1Ts1z1Ts\begin{align} s \approx z\frac{1-z^{-1}}{T_{\rm s}} \approx \frac{1-z^{-1}}{T_{\rm s}} \end{align}
以上を用いて,zz の関数として記述される近似伝達関数 GdG_{\rm d} を得る。
Gd(z)=G(s)s=1z1Ts\begin{align} G_{\rm d}(z) = G(s)|_{s=\frac{1-z^{-1}}{T_{\rm s}}} \end{align}

双一次変換

シフトオペレータの定義より,ss を逆双曲正接関数を用いて表現する。
s=1Tsln(z)=2Tstanh1(1z11+z1)\begin{align} s = \frac{1}{T_{\rm s}}\ln(z) = \frac{2}{T_{\rm s}}\tanh^{-1}\left(\frac{1-z^{-1}}{1+z^{-1}} \right) \end{align}
ここで,逆双曲正接関数の近似多項式は次のように記述される。
tanh1(w)=n=0w2n+12n+1, w<1w=1z11+z1\begin{align} \tanh^{-1}(w) &= \sum_{\rm n=0}^{\infty} \frac{w^{2{\rm n}+1}}{2{\rm n}+1},\ |w|<1 \\ w &=\frac{1-z^{-1}}{1+z^{-1}} \end{align}
この近似は z=1z=1 近傍において成立する。したがって,ss の近似多項式は次のように記述される。
s2Tsn=012n+1(1z11+z1)2n+1\begin{align} s \approx \frac{2}{T_{\rm s}}\sum_{\rm n=0}^{\infty} \frac{1}{2{\rm n}+1} \left(\frac{1-z^{-1}}{1+z^{-1}} \right)^{2{\rm n}+1} \end{align}
この近似を1次で打ち切ったものが双一次変換と呼ばれる。
s2Ts1z11+z1\begin{align} s \approx \frac{2}{T_{\rm s}}\frac{1-z^{-1}}{1+z^{-1}} \end{align}
以上を用いて,zz の関数として記述される近似伝達関数 GdG_{\rm d} を得る。
Gd(z)=G(s)s=2Ts1z11+z1\begin{align} G_{\rm d}(z) = G(s)|_{s=\frac{2}{T_{\rm s}}\frac{1-z^{-1}}{1+z^{-1}}} \end{align}

周波数ワーピング

離散システム GdG_{\rm d} は近似によって得られたシステムであるため,その特性は GG と一致しない。特に周波数軸の歪み(ワープ)が生じるため,プリワーピングと呼ばれる処理が必要となる。これは G,GdG, G_{\rm d} のそれぞれの周波数伝達関数から確認することができる。
周波数伝達関数を確認するため,以下では各システムが jωj\omega の関数であるとして考える。GdG_{\rm d}GG によって表すと以下のようになる。
Gd(jω)=G(2Ts1ejωTs1+ejωTs)=G(2TsejωTs/2ejωTs/2ejωTs/2+ejωTs/2)=G(2Ts2jsin(ωTs/2)2cos(ωTs/2))=G(j2Tstan(ωTs/2))=G(jωa)ωa2Tstan(ωTs/2)\begin{align} G_{\rm d}(j\omega) &= G\left(\frac{2}{T_{\rm s}}\frac{1-e^{-j\omega T_{\rm s}}}{1+e^{-j\omega T_{\rm s}}}\right) \\ &= G\left(\frac{2}{T_{\rm s}}\frac{e^{j\omega T_{\rm s}/2}-e^{-j\omega T_{\rm s}/2}}{e^{j\omega T_{\rm s}/2}+e^{-j\omega T_{\rm s}/2}}\right) \\ &= G\left(\frac{2}{T_{\rm s}}\frac{2j\sin(\omega T_{\rm s}/2)}{2\cos(\omega T_{\rm s}/2)}\right) \\ &= G\left(j\frac{2}{T_{\rm s}}\tan(\omega T_{\rm s}/2)\right) \\ &= G\left(j\omega_{\rm a}\right) \\ \omega_{\rm a} &\equiv \frac{2}{T_{\rm s}}\tan(\omega T_{\rm s}/2) \end{align}
これは GdG_{\rm d} の周波数伝達関数が「周波数がワーピングされた空間における GG の周波数伝達関数」と一致することを意味している。低周波帯域ではこのワーピングは問題にならないが,高周波帯域では振幅・位相特性が大きくずれるため,GdG_{\rm d} は設計段階の想定通りに動作しない。そのため,予め周波数ワーピングを考慮して GG を設計する必要があり,この処理をプリワーピングと呼ぶ。
また,この周波数歪み歪みのため GGGdG_{\rm d}の極と零点は一致しないことに留意する。この変換では直流から中周波帯域までの周波数特性を精度良く変換しており,内部変数としての極と零点の位置には焦点を当てていない。

エイリアシングフリー

周波数ワーピングは高周波帯域の周波数を歪めるため伝達関数の特性を歪める一方で,エイリアシングを回避するといった利点も存在する。本変換によって ss 平面の虚軸上の範囲 (j, j)(-j\infty,\ j\infty) の点は zz 平面上の単位円上に移り,その際の位相は(π, π)(-\pi,\ \pi)に収まる。
ss 平面の虚軸上の任意の点 jωj\omegazz 平面上の点 ejΩe^{j\Omega} に変換する作用を考える。双一次変換の定義式より,これらの点は次のような関係にある。
jω=2Ts1ejΩ1+ejΩ=j2Tstan(Ω/2)Ω=2arctan(ωTs/2)\begin{align} j\omega &= \frac{2}{T_{\rm s}}\frac{1-e^{-j\Omega}}{1+e^{-j\Omega}} = j\frac{2}{T_{\rm s}}\tan(\Omega/2) \\ \therefore \Omega &= 2\arctan \left(\omega T_{\rm s}/2\right) \end{align}
ss 平面虚軸上の点は単射かつ単調増加する非周期関数である逆正接関数によって一意に変換され,その範囲は Ω(π, π)\Omega\in(-\pi,\ \pi) となるため,離散化に伴うエイリアシングは発生しない。すなわち,任意の極零配置を持つシステムについて,ナイキスト周波数の制約なしに離散化して問題ない。

整合z変換

本変換法では GdG_{\rm d} の極と零点が GG の極と零点を持つように設計する。ここでは GG を零点 bk (k=1p)b_k\ (k=1\dots p) と極 ak (k=1q)a_k\ (k=1\dots q) を持つ有理関数伝達関数として表現する。
G(s)=Kk=1p(sbk)k=1q(sak)\begin{align} G(s) &= K\frac{\prod_{\rm k=1}^{p} (s-b_k)}{\prod_{\rm k=1}^{q}(s-a_k)} \end{align}
このシステムと同一の極と零点を持つ離散系システムの最小実現は次のようになる。
Gd(s)=Kdk=1p(1ebkTsesTs)k=1q(1eakTsesTs)Kd=G(0)k=1q(1eakTs)k=1p(1ebkTs)\begin{align} G_{d}(s) &= K_{\rm d}\frac {\prod_{\rm k=1}^{p} (1-e^{b_kT_{\rm s}}e^{-sT_{\rm s}})}{\prod_{\rm k=1}^{q}(1- e^{a_kT_{\rm s}}e^{-sT_{\rm s}})} \\ K_{\rm d} &= |G(0)| \left| \frac{ \prod_{\rm k=1}^{q}(1-e^{a_kT_{\rm s}})} {\prod_{\rm k=1}^{p}(1-e^{b_kT_{\rm s}})} \right| \end{align}
これを zz の関数に変形して,離散時間系で実装可能な形式として記述する。
Gd(z)=Kdk=1p(1ebkTsz1)k=1q(1eakTsz1)\begin{align} G_{d}(z) &= K_{\rm d}\frac {\prod_{\rm k=1}^{p} (1-e^{b_kT_{\rm s}}z^{-1})}{\prod_{\rm k=1}^{q}(1- e^{a_kT_{\rm s}}z^{-1})} \end{align}
本変換法は極零配置が重要なシステムにて使用することが適している。ただし,ナイキスト周波数を超える範囲に存在する極と零点については,離散化に伴ってエイリアシングの影響を受けることに留意する。これは本手法における szs\rightarrow z が被覆写像であるためであり,たとえば ss 平面上の以下の極は全て zz 平面上の同一点に写像される。
ak=ωr+j(ωi+2πTsk), kZzp,k=eakTs=eωrTs+j(ωrTs+2πk)=e(ωr+jωr)Ts=ea0Ts\begin{align} a_{\rm k} &= -\omega_{r} + j\left(\omega_{i} + \frac{2\pi}{T_{\rm s}}{\rm k}\right),\ {\rm k}\in\mathbb{Z} \\ z_{\rm p,k} &= e^{a_{\rm k}T_{\rm s}} \\ &= e^{-\omega_{r}T_{\rm s} + j(\omega_{r}T_{\rm s} + 2\pi{\rm k})} \\ &= e^{(-\omega_{r} + j\omega_{r})T_{\rm s}} \\ &= e^{a_{0}T_{\rm s}} \end{align}

離散化システムの実装例

ここでは双一次変換を用いて簡単な低域通過フィルタを実装する例を示す。 以下の低域通過特性を示す伝達関数を考える。
Gc(s)=gs+g\begin{align} G_{\rm c}(s) &= \frac{g}{s + g} \end{align}
ここで,gR>0g\in \mathbb{R} > 0 は遮断周波数 [rad/s][{\rm rad/s}] を表す。離散化伝達関数は次のように表現される。
Gd(z)=Gc(2Ts1z11+z1)=gTs(1+z1)(2+gTs)+(2+gTs)z1\begin{align} G_{\rm d}(z) &= G_{\rm c}\left (\frac{2}{T_{\rm s}}\frac{1 - z^{-1}}{1 + z^{-1}} \right) \\ &= \frac{gT_{\rm s}(1 + z^{-1})}{(2+gT_{\rm s}) + (-2+gT_{\rm s})z^{-1}} \end{align}
このフィルタの入力を xx,出力を yy とすると,以下の関係が成立する。
((2+gTs)+(2+gTs)z1)Y(s)=gTs(1+z1)X(s) Y(s)=gTs(1+z1)2+gTsX(s)+2gTs2+gTsz1Y(s) y[k]=gTs2+gTs(x[k]+x[k1])+2gTs2+gTsy[k1]\begin{align} &\left((2+gT_{\rm s}) + (-2+gT_{\rm s})z^{-1}\right)Y(s)= gT_{\rm s}(1 + z^{-1})X(s)\\ \Leftrightarrow\ &Y(s) = \frac{gT_{\rm s}(1 + z^{-1})}{2+gT_{\rm s}}X(s) + \frac{2 - gT_{\rm s}}{2+gT_{\rm s}}z^{-1}Y(s)\\ \Leftrightarrow\ &y[k] = \frac{gT_{\rm s}}{2+gT_{\rm s}}\left(x[k] + x[k-1]\right) + \frac{2 - gT_{\rm s}}{2+gT_{\rm s}}y[k-1]\\ \end{align}