離散後の伝達関数の実装方法 覚え書き(c2d -> 実装)

連続系伝達関数(s領域)の離散化はMATLABやPython.controlのc2dコマンドで簡単に出来るけれど、
でどうやって組み込みソフトの周期タスク内に落とし込む形にするのかという話のメモ書きです。



(手計算で離散化させれば実装まで持っていけるのに、c2dコマンドで出てきたものをどう扱ってよいか恥ずかしながら迷子に。)

離散化まで

例として単純な積分器$(\frac{1}{s})$、対象のタスクは1秒周期(sample_time=1)とします。
離散化はc2dコマンド一発で

MATLABの場合は

sample_time = 1;
s= tf('s');
integral_c = 1/s;
integral_d = c2d(integral_c, sample_time, "tustin")

Pythonの場合は

import control
from control import tf, c2d

sample_time = 1
s= tf('s')
integral_c = 1/s
integral_d = c2d(integral_c, sample_time, method='tustin')

出力はどちらを使っても

integral_d =
  0.5 z + 0.5
  -----------
     z - 1
Sample time: 1 seconds
Discrete-time transfer function.

こちらの出てきた式は出力$Y$,入力$U$とする離散系伝達関数で下記を指します。
$$ \frac{Y(z)}{U(z)} = \frac{0.5z + 0.5}{z - 1} $$
これをどう実装できる形にするのかが分からなくなるポイントで、自分は迷子になりました。

実装可能な式へ

ここから$z$の次数を下げてから差分方程式にします。
(分母より分子の次数が高い場合は非プロパーで実装不可なので下には進めません)

まず、右辺の分母分子両方を(分母の)$z$の最高次数で割ってzの最高次数を0とし実装可能な形にします。
この場合は$z$で割るので、
$$ \frac{Y(z)}{U(z)} = \frac{\frac{0.5z}{z} + \frac{0.5}{z}}{\frac{z}{z} - \frac{1}{z}} $$
よって
$$ \frac{Y(z)}{U(z)} = \frac{0.5 + \frac{0.5}{z}}{1 - \frac{1}{z}} $$
左辺を$Y(z)$とするように式変形していくと
$$ Y(z) \times (1 - \frac{1}{z}) = (0.5 + \frac{0.5}{z} ) \times U(z)$$

$$ Y(z) = Y(z) \times \frac{1}{z} + 0.5 \times U(z) + 0.5 \times \frac{1}{z} \times U(z) $$
ここまで変わらず出力は$Y$, 入力は$U$です。

次に、時間領域tへの書き換えを実施します。
これも単純で、$z$がかかった項は今回値の$t$となり、$z-1$がかかった項は前回値$t-1$になります。

具体的には

• $Y(z)$ は　$Y(t)$ : 今回の出力
• $Y(z)\times\frac{1}{z}$ 　は　 $Y(t-1)$ : 1回前の出力
• $U(z)$ は $U(t)$ : 今回の入力
  になります。

したがって、元の式 $$ Y(z) = Y(z) \times \frac{1}{z} + 0.5 \times U(z) + 0.5 \times \frac{1}{z} \times U(z) $$ を変換すると $$ Y(t) = Y(t-1) + 0.5 \times U(t) + 0.5 \times U(t-1) $$ になりました。

日本語にすると

今回出力 = 前回出力 + (0.5 * 今回入力) + (0.5 * 前回入力)

C言語ライクに書くと下記です。

float Y = 0f;　    //出力値
float Y_z1 = 0f;  // 出力の前回値
float U = 0f;     //入力値
float U_z1 = 0f;  //入力の前回値

Task_1sec()　//離散化時のサンプル時間ごとに実行される関数
{
    //入力の取得
    U = funcGetU();
    printf("入力U=%f",U);
    
    //Yの算出
    Y = Y_z1 + (0.5f * U) + (0.5f * U_z1);

    //出力
    printf("出力Y=%f",Y);

    //前回値の更新
    U_z1 = U;
    Y_z1 = Y;
}