この記事では、離散時間の伝達関数について説明します。

はじめに

前回の記事では、ラプラス変換を用いて導出された伝達関数の概念について説明しました。

しかし、ラプラス変換から得られる伝達関数は連続時間システムで使用するものであり、離散時間システムで使用する場合はZ変換を用います。Z変換に関する内容も以前の記事で説明しました。

今回は、上記の伝達関数の記事とZ変換の記事で理解した知識を基礎として、離散時間の伝達関数について説明します。

運動方程式の伝達関数の例

まず、最も簡単な運動方程式を見てみましょう。時間$t$の関数である入力の力$f(t)$を受けながら、空気摩擦抵抗係数$c_v$の空気中を走る質量$m$の薄い空中円筒の速度$v(t)$の運動方程式は、次のように書けます（詳しい説明は前回の記事に記載したため省略します）。

\begin{align}
\dot{v}(t) &= \frac{f(t)-c_v v(t)}{m}
\end{align}

ここでは、まずオイラー法を用いて微分方程式から漸化式を導きます。

\begin{align}
\dot{v}_{n} = \frac{v_{n+1} - v_n}{\Delta_T}
\end{align}

$\dot{v}_{n}$を代入し、Z変換を行い、線形性とシフト性を用いて、$z$空間における$v_n \to V(z)$と$f_n \to F(z)$の関係を求めます。

\begin{align}
\frac{v_{n+1} - v_{n}}{\Delta_T} &= \frac{f_n - cv_n}{m} \\
\frac{\mathcal{Z}[v_{n+1}] - \mathcal{Z}[v_{n}]}{\Delta_T} &= \frac{\mathcal{Z}[f_n] - c_v\mathcal{Z}[v_n]}{m} \\
\frac{zV(z) - V(z)}{\Delta_T} &= \frac{F(z) - c_vV(z)}{m} \\
(mz - m + c_v\Delta_T)V(z) &= \Delta_TF(z) \\
V(z) &= \frac{\Delta_T}{mz - m + c_v\Delta_T}F(z) \\
\end{align}

これにより、ラプラス変換と同様に、$z$で表される伝達関数が得られます。

G(z) = \frac{\Delta_T}{mz - m + c_v\Delta_T}

尚、ここで用いたシフト性は、初期値を0としたものです。つまり

\mathcal{Z}[x_{n+k}] = z^kX(z)

ラプラス変換の時と同様に、伝達関数を用いる問題では通常、初期値は無視されます。

離散伝達関数の一般的な形

次に、より一般的な形式について見てみましょう。

通常、伝達関数で分析しようとするシステムは、入力$u$と出力$y$の線形関係を表す漸化式で次のように表現できます。

\begin{align}
& y_{n+i} + a_{i-1}y_{n+i-1} + \cdots + a_1y_{n+1} + a_0y_{n} \\
&= b_{j}u_{n+j} + b_{j-1}u_{n+j-1} + \cdots + b_1u_{n+1} + b_0u_{n}
\end{align}

ここでZ変換を適用し、シフト性$\mathcal{Z}[x_{n+k}] = z^kX(z)$を用います。

\begin{align}
      & z^{i} Y(z) + a_{i-1}z^{i-1}Y(z) + \cdots + a_1 zY(z) + a_0 Y(z) \\
      &= b_j z^{j} U(z) + b_{j-1}z^{j-1}U(z) + \cdots + b_1zU(z) + b_0U(z)
\end{align}

これにより、伝達関数の一般的な形が得られます。

\begin{align}
      G(z) = \frac{Y(z)}{U(z)} &= \frac{b_j z^{j} + b_{j-1}z^{j-1} + \cdots + b_1z + b_0}{z^{i} + a_{i-1}z^{i-1} + \cdots + a_1 z + a_0}
\end{align}

連続時間システムが$s$を用いて伝達関数を表現するのと同様に、離散時間システムでは$z$を用いてこのように伝達関数を表現できます。

Pythonのcontrolモジュールによる離散時間伝達関数の実装

離散時間の伝達関数の実装は、連続時間の場合と同じtf関数を使用しますが、分子と分母の後に時間間隔dt（方程式では$\Delta_T$）を指定します。

使用例：

import control

bunshi = [2,3]
bunbo = [4,5,6]
dt = 0.1
Gz = control.tf(bunshi,bunbo,dt)
print(Gz)

結果：

dt = 0.1

      2 z + 3
  ---------------
  4 z^2 + 5 z + 6

これにより、表記が$s$ではなく$z$になります。

連続時間と同様に、forced_response関数を使用して入力から出力を求めることができます。

import numpy as np
import matplotlib.pyplot as plt
import control

m = 500
c_v = 60
Delta_T = 0.5
bunbo = [m,c_v*Delta_T-m]
Gz = control.tf([Delta_T],bunbo,Delta_T)

n = np.arange(101)
t = n*Delta_T
f = np.zeros(len(n))
f[(t>2)&(t<40)] = 800
v = control.forced_response(Gz,t,f).outputs

plt.figure(figsize=[6,6],dpi=100,facecolor='#eeeeff')
plt.subplot(211)
plt.scatter(t,f,5,'m')
plt.ylabel(r'$f_n$',size=14)
plt.grid(ls=':')
plt.subplot(212)
plt.scatter(t,v,5,'r')
plt.ylabel(r'$v_n$',size=14)
plt.xlabel(r'$t = \Delta_Tn$',size=14)
plt.grid(ls=':')
plt.show()

ここで、forced_responseの入力はnではなくtを使用することに注意してください。また、その時間間隔は、この伝達関数の定義時に指定した時間間隔と一致するか、その倍数である必要があります。

連続時間の伝達関数から離散時間の伝達関数への変換

以上では、漸化式からZ変換を行って離散時間の伝達関数を求める方法を紹介しましたが、実際には連続時間の伝達関数から変換する方法の方がよく用いられます。

しかし、面倒なことに、伝達関数を変換する方法は一つではなく多数存在します。それぞれの方法は異なる状況で有効であるため、各方法を理解した上で使い分ける必要があります。

ここでは、代表的な4つの方法を紹介します。これらの方法は、PythonのcontrolモジュールやMATLABでも実装できます。

ゼロ次ホールド

ゼロ次ホールドの変換方法

これは最も基本的で、最もよく使われる方法です。

まずは使用方法から紹介し、理論や理由については後で説明します。

連続時間の伝達関数$G(s)$があるとき、離散時間の伝達関数$G(z)$は次のように求められます。

G(z) = \left(\frac{z - 1}{z}\right)\mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right]\right]

ゼロ次ホールドの使用例

前述と同じ、空気摩擦抵抗を受ける空中円筒の制御を考えます。前回の記事で述べたように、連続時間の伝達関数は次のようになります。

G(s) = \frac{1}{ms + c_v}

この$G(s)$から出発して、順を追って$G(z)$を導出します。

\begin{align}
\frac{1}{s}G(s) &= \frac{1}{s(ms + c_v)} = \frac{1}{c_v}\left(\frac{1}{s}-\frac{1}{s + c_v/m}\right) \\
\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right] &= \frac{1}{c_v} - \frac{1}{c_v}e^{-c_vt/m} = \frac{1}{c_v}\left(1 - \left(e^{-c_v\Delta_T/m}\right)^n\right) \\
\mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right]\right] &= \frac{1}{c_v}\left(\frac{z}{z-1} - \frac{z}{z-e^{-c_v\Delta_T/m}}\right) \\
\left(\frac{z - 1}{z}\right)\mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right]\right] &= \frac{1}{c_v}\left(1 - \frac{z-1}{z-e^{-c_v\Delta_T/m}}\right) \\
G(z) &= \frac{1-e^{-c_v\Delta_T/m}}{c_v(z-e^{-c_v\Delta_T/m})}
\end{align}

ゼロ次ホールドの考え方

ここからは少し複雑な話になります。

ゼロ次ホールドの考え方は、離散時間の入力から、できるだけ元の連続時間での計算と同じ出力を得たいというものです。

前回の記事でも述べたように、ある連続時間関数があり、それを離散化してラプラス変換すると、それはZ変換になります。

具体的には、まず元の入力の連続時間関数を$u_a(t)$とし、サンプリング周期$n\Delta_T$の時点でのみ値を持つ関数を$u_n = u_a(n\Delta_T)$とします。この$u_n$は、くし型関数を用いて「インパルス列」としての連続時間関数$u_d(t)$で次のように表現できます。

\begin{align}
u_d(t) &= u_a(t)Ш(t) \\
&= u_a(t)\sum_{n=-\infty}^{\infty}\delta(t-n\Delta_T)
\end{align}

これにより、

\begin{align}
\mathcal{L}[u_d(t)] &= \mathcal{Z}[u_n] \\
U_d(s) &= U(z)
\end{align}

と書くことができます。

ここで、連続時間で入力$U_a(s)=\mathcal{L}[u_a(t)]$に対する出力$Y_a(s)$を求める伝達関数$G(s)$があるとします。

\begin{align}
Y_a(s) &= G(s)U_a(s)
\end{align}

この同じ伝達関数を用いて、$U_d(s)$に対する出力$Y_d(s)$を求めることもできます。

\begin{align}
Y_d(s) &= G(s)U_d(s)
\end{align}

ここで、$U_d(s)$はサンプリング時点$n\Delta_T$のみで値を持つインパルス入力ですが、出力$Y_d(s)$はそうではないことに注意してください。なぜなら、たとえインパルス入力であっても積分が行われ、連続値の出力が得られるからです。

このような$Y_d(s)$でも、$Y_a(s)$の近似としてはある程度使用できますが、誤差が大きくなります。

私達が求めたいのは、できるだけ本来の連続値$u_a(t)$から計算した出力$y_a(t)$に近い出力です。そのためには、できるだけ$u_a(t)$に近い入力が望ましいですが、実際に利用できるのはサンプリング時点でのみ値を持つ$u_d(t)$（または$u_n$）です。そこで、この$u_d(t)$から、$u_a(t)$を近似する新しい連続値$u_c(t)$を生成する必要があります。

そして、最も一般的で簡単な近似方法は、直前のサンプリング値を周期が終わるまで保持し続けるという方法です。つまり、次の図のようになります。

この方法はゼロ次ホールド（Zero-Order Hold、略してZOH）と呼ばれます。

これは矩形関数と畳み込みを行うことで実現できます。

矩形関数は、2つの単位ステップ関数を用いて次のように定義できます。

h_{\text{zoh}}(t) = \theta(t)-\theta(t-\Delta_T)

ここでステップ関数の定義は：

\begin{align}
      \theta(t-\alpha) = \left\{ \begin{matrix} 1 & (t \ge \alpha) \\
       0 & (t < \alpha) \end{matrix} \right.
\end{align}

$u_c(t)$は、この$h_{\text{zoh}}(t)$との畳み込み、つまり次のように計算されます。

\begin{align}
u_c(t) &= u_d(t)*h_{\text{zoh}}(t) \\
&= \int_{-\infty}^{\infty} u_d(\tau)h_{\text{zoh}}(t-\tau)d\tau
\end{align}

畳み込みに関する詳細な解説は省略しますが、簡単にまとめると、次のGIF画像のように説明できます。

上側が$u_d(t)$で、移動する矩形が$h_{\text{zoh}}(t)$を表し、下側のような$u_c(t)$が生成されます。

畳み込みは一見複雑に見えますが、ラプラス変換すると、畳み込みは単なる乗算になります。

\begin{align}
\mathcal{L}[a(t)*b(t)] &= \mathcal{L}[a(t)]\mathcal{L}[b(t)] \\
&= A(s)B(s) \\
\end{align}

これもラプラス変換の重要な特性の一つです。

したがって、

\begin{align}
U_c(s) &= \mathcal{L}[h_{\text{zoh}}(t)*u_d(t)] \\
&= \mathcal{L}[h_{\text{zoh}}(t)]\mathcal{L}[u_d(t)] \\
&= H_{\text{zoh}}(s)U_d(s)
\end{align}

ここで、

\begin{align}
H_{\text{zoh}}(s) &= \mathcal{L}[h_{\text{zoh}}(t)]\\
&= \mathcal{L}\left[\theta(t)-\theta(t-\Delta_T)\right] \\
&= \frac{1}{s} - \frac{e^{-\Delta_T s}}{s} \\
&= \frac{1-e^{-\Delta_T s}}{s}
\end{align}

となります。したがって、出力$Y_c(s)$は、入力$U_d(s)$から次のように計算できます。

\begin{align}
Y_c(s) &= G(s)U_c(s)\\
&= G(s)H_{\text{zoh}}(s)U_d(s)
\end{align}

そして、新しい伝達関数として次のように書くことができます。

\begin{align}
Y_c(s) &= G_{\text{zoh}}(s)U_d(s)
\end{align}

\begin{align}
G_{\text{zoh}}(s) = G(s)H_{\text{zoh}}(s)
\end{align}

さて、ここまではラプラス変換の範囲でしたが、いよいよZ変換の話に戻ります。私たちが求めたいのは、この$G_{\text{zoh}}(s)$に対応するz領域の伝達関数です。

これは次のように求められます。

\begin{align}
G(z) &= \mathcal{Z}[\mathcal{L}^{-1}[H_{\text{zoh}}(s)G(s)]] \\
&= \mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{1 - e^{-\Delta_T s}}{s}G(s)\right]\right]
\end{align}

ただし、$z=e^{\Delta_T s}$の関係が知られているので、これを代入し、$1 - e^{-\Delta_T s} = 1-z^{-1}$を外に出すことで、次のように書き換えられます。

\begin{align}
G(z)  &= \left(1-z^{-1}\right)\mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right]\right] \\
&= \left(\frac{z - 1}{z}\right)\mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right]\right]
\end{align}

以上が、ゼロ次ホールド法の説明でした。

ステップ関数入力には最強

ゼロ次ホールド法は、「ステップ応答を完璧に再現できる方法」と説明されることがあります。これもゼロ次ホールド法の考え方の一つです。

ステップ応答とは、ステップ関数入力に対する応答です。つまり、$u(t)$は$t \ge 0$で1（または時間によらない何らかの定数）である入力です。

この場合、ラプラス変換は$\begin{align} U(s) = \frac{1}{s} \end{align}$、Z変換は$\begin{align} U(z) = \frac{z}{z-1} \end{align}$となります。そして、それぞれ伝達関数を乗じることで得られる出力は、$\begin{align} Y(s) = \frac{1}{s}G(s) \end{align}$、$\begin{align} Y(z) = \frac{z}{z-1}G(z) \end{align}$となります。

一方、ゼロ次ホールドにおける$G(s)$と$G(z)$の関係から、次のように書き換えられます。

\begin{align}
G(z)  &= \left(\frac{z - 1}{z}\right)\mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right]\right] \\
\frac{z}{z - 1} G(z) &= \mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right]\right] \\
\mathcal{Z}^{-1}\left[\frac{z}{z - 1} G(z)\right] &= \mathcal{L}^{-1}\left[\frac{1}{s}G(s)\right] \\
\end{align}

これはステップ関数入力に対する出力と一致します。したがって、ゼロ次ホールドから得られた$G(z)$は、ステップ関数入力の場合を完璧に再現できるわけです。

一次ホールド

ゼロ次ホールド法の改良版として、一次ホールド（First-Order Hold、略してFOH）という方法があります。三角形近似とも呼ばれます。

ゼロ次ホールドと似ていますが、連続入力$u_c(t)$を生成するための近似方法が異なります。ゼロ次ホールドでは各周期で一定値を使用するのに対し、一次ホールドでは前後の値の両方を用いて線形計算を行います。つまり、次の図のように、各サンプリング点を直線で結ぶのです。

そのため、ゼロ次ホールドでは$H_{\text{zoh}}(s)$を使用するのに対し、一次ホールドではより複雑な$H_{\text{foh}}(s)$を使用します。

H_{\text{foh}}(s) = \frac{\Delta_Ts+1}{\Delta_T}\left(\frac{1-e^{-\Delta_Ts}}{s}\right)^2

導出過程は省略しますが、ゼロ次ホールドと同様に、一次ホールドでは次のように$G(z)$を求めます。

G(z) = \frac{1}{\Delta_T}\left(\frac{z - 1}{z}\right)^2\mathcal{Z}\left[\mathcal{L}^{-1}\left[\frac{\Delta_Ts+1}{s^2}G(s)\right]\right]

より高度な近似を使用するため、一次ホールドはゼロ次ホールドよりも精度が高いですが、このように複雑で計算が難しいため、実際には計算が容易なゼロ次ホールドの方が一般的に使用されるのが現状です。

双一次変換

双一次変換方法

ゼロ次ホールドとも一次ホールドとも異なるアプローチとして、双一次変換（タスティン変換又は台形差分法とも）という方法があります。

この方法では、逆ラプラス変換もZ変換も必要ありません。単に$G(s)$の中の$s$に

\begin{align}
s &= \frac{2}{\Delta_T} \frac{z - 1}{z + 1} \\
\end{align}

を代入するだけで$G(z)$が得られます。

双一次変換の使用例

上記と同じ、空気摩擦抵抗を受ける空中円筒の制御問題を例にします。

\begin{align}
G(s) &= \frac{1}{ms + c_v}
\end{align}

双一次変換により、次のように$G(z)$が求められます。

\begin{align}
G(z) &= G(s)\left| _{s=\frac{2}{\Delta_T} \frac{z - 1}{z + 1}} \right. \\
&= \frac{1}{m\frac{2}{\Delta_T} \frac{z - 1}{z + 1} + c_v} \\
&= \frac{1}{\frac{2m}{\Delta_T} \frac{z - 1}{z + 1} + \frac{c_v\Delta_T}{\Delta_T}\frac{z + 1}{z + 1}} \\
&= \frac{\Delta_T(z + 1)}{2m(z - 1) + c_v\Delta_T(z + 1)} \\
&= \frac{\Delta_Tz + \Delta_T}{(2m + c_v\Delta_T)z - 2m + c_v\Delta_T}
\end{align}

双一次変換の考え方

双一次変換の考え方は、台形公式に基づいています。

例えば、位置$x_n$の微分である速度$v_n$があるとして、次の時刻の変位$x_{n+1}$は、前回の記事で紹介したオイラー法では$x_{n+1} = x_n + v_n\Delta_T$となりますが、台形公式では次のようになります。

x_{n+1} = x_n + \frac{(v_{n+1} + v_n)\Delta_T}{2}

Z変換を適用し、z空間での関係を考慮します。

\begin{align}
zX(z) &= X(z) + \frac{(zV(z) + V(z))\Delta_T}{2} \\
(z - 1)X(z) &= \frac{(z + 1)\Delta_T}{2}V(z) \\
V(z) &= \frac{2}{\Delta_T}\frac{(z - 1)}{(z + 1)}X(z) \\
\end{align}

速度と位置の関係はs空間では$V(s) = sX(s)$です。

したがって、

\begin{align}
s &= \frac{2}{\Delta_T}\frac{(z - 1)}{(z + 1)} \\
\end{align}

と考えることができます。

ちなみに、もう一つの導出方法があります。まず、本来$s = \frac{1}{\Delta_T} \ln(z)$という関係があります。ここで、$\ln$は$\text{arctanh}$と次の関係があります。

\begin{align}
\text{arctanh}(x) &= \frac{1}{2}\ln\left(\frac{1+x}{1-x}\right) \\
\ln(y) &= \text{2 arctanh}\left( \frac{y-1}{y+1}\right)
\end{align}

$\ln$を$\text{arctanh}$で置き換え、テイラー展開（参考）を行い、第1項のみを採用します。

\begin{align}
s &= \frac{1}{\Delta_T} \ln(z)  \\
  &= \frac{2}{\Delta_T} \text{arctanh}\left( \frac{z-1}{z+1}\right) \\
  &= \frac{2}{\Delta_T} \left[\frac{z-1}{z+1} + \frac{1}{3} \left( \frac{z-1}{z+1} \right)^3  + \frac{1}{5} \left( \frac{z-1}{z+1} \right)^5 + \cdots \right] \\
  &\approx  \frac{2}{\Delta_T} \frac{z - 1}{z + 1} \\
\end{align}

インパルス不変法

ゼロ次ホールドを使用する場合、元の伝達関数$G(s)$に$H_{\text{zoh}}(s)$を乗じることで、元の連続時間関数での出力をより良く近似できると説明しました。しかし、$H_{\text{zoh}}(s)$を乗じず、そのままの$G(s)$を使用しても、ある程度は使用できます。つまり、単に

G(z) = \mathcal{Z}\left[\mathcal{L}^{-1}\left[G(s)\right]\right]

とします。これはインパルス不変法と呼ばれます。名前の通り、入力がインパルス信号の場合に最も有効です。なぜなら、入力が本当にインパルスである場合、わざわざ連続値に変換する必要がないからです。

したがって、インパルス入力を扱う場合に限り、このような単純なインパルス不変法を使用することが適切と言えます。

上記と同じ例に適用したら、このように$G(z)$を求めることができます。

\begin{align}
G(s) &= \frac{1}{ms + c_v} \\
\mathcal{L}^{-1}[G(s)] &= \frac{1}{m}e^{-c_vt/m} = \frac{1}{m}\left(e^{-c_v\Delta_T/m}\right)^n \\
G(z) = \mathcal{Z}\left[\mathcal{L}^{-1}\left[G(s)\right]\right] &= \frac{z}{m(z-e^{-c_v\Delta_T/m})}
\end{align}

各方法の比較

以上で紹介した4つの方法を纏めます。

方法	得意な入力	$\begin{align}\frac{1}{s}\end{align}$	$\begin{align}\frac{1}{s^2}\end{align}$
ゼロ次ホールド	ステップ関数	$\begin{align}\frac{\Delta_T}{z-1}\end{align}$	$\begin{align}\frac{\Delta_T^2}{2}\frac{z+1}{z^2-2z+1}\end{align}$
一次ホールド	ランプ関数	$\begin{align}\frac{\Delta_T}{2}\frac{z+1}{z-1}\end{align}$	$\begin{align}\frac{\Delta_T^2}{6}\frac{z^2+4z+1}{z^2-2z+1}\end{align}$
双一次変換	特になし	$\begin{align}\frac{\Delta_T}{2}\frac{z+1}{z-1}\end{align}$	$\begin{align}\frac{\Delta_T^2}{4}\frac{z^2+2z+1}{z^2-2z+1}\end{align}$
インパルス不変法	インパルス関数	$\begin{align}\frac{\Delta_Tz}{z-1}\end{align}$	$\begin{align}\frac{\Delta_T^2z}{z^2-2z+1}\end{align}$

普段よく使われているのはゼロ次ホールドと双一次変換なので、迷った場合はこの2つの方法のみを考慮すれば十分です。ステップ関数に近い（変動が小さい）入力を扱う場合、ゼロ次ホールドが最適ですが、変化が激しい入力の場合には双一次変換の方が優れています。

一次ホールドはゼロ次ホールドと双一次変換より精度が高いですが、複雑すぎる上に、ゼロ次ホールドと双一次変換で十分な場合が多いため、あまり使用されません。ただし、非常に高い精度が要求され、ゼロ次ホールドと双一次変換だけでは不十分な場合には使用されることもあります。

インパルス不変法はゼロ次ホールドより簡単ですが、良い結果を出すのはインパルス入力を扱う場合くらいです。それ以外の場合は通常、使用されません。

Pythonのcontrolモジュールによる離散時間伝達関数への変換

c2d（continuous to discrete、つまり「連続から離散へ」の略）関数を使用すると、既存の連続時間伝達関数から簡単に離散時間伝達関数に変換できます。

使用例：

Gs = control.tf([1],[1,0])
dt = 0.2
Gz = control.c2d(Gs,dt,'foh')
print(Gz)

結果

dt = 0.2

  0.1 z + 0.1
  -----------
     z - 1

ここで必要な入力は、時間間隔dtと変換に使用するメソッドです。メソッドの指定は以下の表の通りです。

zoh	ゼロ次ホールド
foh	一次ホールド
tustin	双一次変換
impulse	インパルス不変

ただし、zohは既定値なので、ゼロ次ホールドを使用したい場合は指定する必要はありません。

ここで、各方法と各種入力に対する比較を行いたいと思います。今回は「バネとダンパーが付いた車」のシステムを例にします。詳細は前の記事を参照してください。このシステムの元の伝達関数は

G(s) = \frac{1}{ms^2 + c_vs + k}

です。この伝達関数をそれぞれの方法で離散時間伝達関数に変換し、計算して結果を比較します。

import numpy as np
import matplotlib.pyplot as plt
import control

k = 5000
m = 50
c_v = 50

t = np.arange(0,10,0.01)
td = np.arange(0,10,0.1)
Gs = control.tf([1],[m,c_v,k])

for inp in ['step','ramp','impulse','cos']:
    if(inp=='step'):
        f = np.zeros(len(t)) + 1000
    elif(inp=='ramp'):
        f = t*100
    elif(inp=='cos'):
        f = 1000*np.cos(t*2)
    fd = f[::10]
    
    if(inp!='impulse'):
        x = control.forced_response(Gs,t,f).outputs
    else:
        x = control.impulse_response(Gs,t).outputs
    
    xd = {}
    for method in ['zoh','foh','tustin','impulse']:
        Gz = control.c2d(Gs,0.1,method)
        if(inp!='impulse'):
            xd[method] = control.forced_response(Gz,td,fd).outputs
        else:
            xd[method] = control.impulse_response(Gz,td).outputs
    
    plt.figure(figsize=[6,8],dpi=100,facecolor='#eeffdd')
    
    if(inp!='impulse'):
        plt.subplot(311,xticklabels=[])
        plt.plot(t,f,'b',label='$f(t)$')
        plt.ylabel(r'$f(t)$',size=14)
        plt.grid(ls=':')
    
    plt.subplot(312,xticklabels=[])
    for method in xd:
        plt.scatter(td,xd[method],5,label=method)
    plt.plot(t,x,'m',label='$x(t)$',lw=1)
    plt.ylabel(r'$x$',size=14)
    plt.legend(loc=4)
    plt.grid(ls=':')
    
    plt.subplot(313)
    for method in xd:
        plt.scatter(td,x[::10]-xd[method],5)
    plt.ylabel(r'$x(t)-x_n$',size=14)
    plt.xlabel(r'$t$',size=14)
    plt.grid(ls=':')
    plt.tight_layout()
    plt.show()

実行すると、4つの入力パターンそれぞれに対するグラフが得られます。それぞれを見ながら比較します。

ステップ関数の場合、やはりゼロ次ホールドが最も優れています。一次ホールドと双一次変換も時間が経つと定常値に収束しますが、インパルス不変法は誤差があって、ほぼ駄目です。

ランプ関数の場合、一次ホールドはほぼ完璧に一致し、双一次変換も時間が経つと誤差がなくなります。一方、ゼロ次ホールドは誤差が残ります。さらに、インパルス不変法は誤差が蓄積して最も悪い結果になります。

インパルス関数の場合、今回だけはインパルス不変法が最も優れています。

正弦波関数（周波数応答）の場合、完璧に再現できる方法はありませんが、一次ホールドは誤差が次第に収束していきます。双一次変換もそれに次いで良好です。

どの方法も、常に最良となるわけではありません。そのため、用途に応じて適切な方法を選択することが重要です。

また、全体的に見ると、一次ホールドが最も良い結果を出す傾向がありますが、計算の複雑さのため、より簡単なゼロ次ホールドと双一次変換の方が一般的に使用されます。

終わりに

本記事では、離散時間システムにおける伝達関数について解説しました。

連続時間伝達関数から離散時間伝達関数への変換には、上述の通り色んな方法がありますが、想定する入力信号の種類や求められる精度に応じて、適切な変換方法を選択することが重要です。特に、ゼロ次ホールドと双一次変換はバランスが良く、多くの実用的な場面で活用されています。

離散時間伝達関数を理解し適切に運用することで、デジタル制御システムの設計と解析が効果的に行えるようになります。

参考