guiメモ

技術寄りの雑記(予定)

パルス伝達関数をプログラムに落とし込む


\def\lr#1{{\left( #1 \right)}}
\def\slr#1{{\left[ #1 \right]}}
\def\eq#1{{\begin{align} #1 \end{align}}}
\def\mat#1{{\left[ \begin{array}{cccccccccc} #1 \end{array} \right}} ]

概要

z領域における伝達関数(パルス伝達関数)をプログラムに落とし込めるよう漸化式の形式に変形する方法についてメモしておきます。 パルス伝達関数をプログラムに落とし込むことに焦点を当てるため、z変換の詳細な性質等には触れません。

本稿はesaに投稿した記事を加筆・修正したものです。

サンプル遅延の表現

唐突かつ天下り的ですが、単位サンプル遅延はz^{-1}で表現されます。(詳しい説明は数ある参考書や他のサイト様にお譲りします。) ここでは例を見た方が理解が早いので簡単な例を示します。あるkサンプル目のデータa\slr{k}に対して単位サンプル遅延z^{-1}を作用させると次のようになります。

\displaystyle{
z^{-1}a\slr{k}=a\slr{k-1}
}

対照的に、zを作用させると単位サンプル進みます。

\displaystyle{
za\slr{k}=a\slr{k+1}
}

パルス伝達関数はこれらの関係を利用することで簡単に漸化式へ変形することができます。

漸化式の一般形と実装

漸化式を求めたいパルス伝達関数H\slr{z}とします。ただしH\slr{z}はプロパーな伝達関数(m\leq n)とします。

\displaystyle{
H\slr{z}=\frac{b_{m}z^m+\cdots+b_1z+b_0}{z^n+a_{n-1}z^{n-1}+\cdots+a_1z+a_0}
}

入力をu\slr{k}、出力をy\slr{k}とすると、伝達関数の定義から次のようになります。

\displaystyle{
\frac{y\slr{k}}{u\slr{k}}
=\frac{b_{m}z^m+\cdots+b_1z+b_0}{z^n+a_{n-1}z^{n-1}+\cdots+a_1z+a_0}
}

分母を払うと、次のようになります。

\displaystyle{
y\slr{k}\left(z^n+a_{n-1}z^{n-1}+\cdots+a_1z+a_0\right)
=u\slr{k}\left(b_{m}z^m+\cdots+b_1z+b_0\right)
}

両辺にz^{-n}を作用させて整理すると、漸化式の一般形が得られます。

\displaystyle{
\begin{align}
y\slr{k}=
&b_mu\slr{k-\left(n-m\right)}+\cdots+b_1u\slr{k-\left(n-1\right)}+b_0u\slr{k-n}\\
&-a_{n-1}y\slr{k-1}-\cdots-a_1y\slr{k-\left(n-1\right)}-a_0y\slr{k-n}
\end{align}
}

これをC++風に落とし込んだ例を以下に示します。

float dtf(float input)
{
    constexpr int n = 1;            // 分母多項式の次数(max(1,m)以上)
    constexpr float a[n] = {};      // 分母多項式の係数(最高次数は含まないことに注意)
    constexpr int m = 0;            // 分子多項式の次数
    constexpr float b[m+1] = {};    // 分子多項式の係数
    
    static float u[m+1] = {};
    static float y[n] = {};
    float output = 0.0f;
    u[m] = input;
    
    // 出力値の計算
    for(int i = 0; i <= m; i++)
    {
        output = output + b[i]*u[i];
    }
    for(int i = 0; i <= n-1; i++)
    {
        output = output - a[i]*y[i];
    }
    
    // 保持する値の更新
    for(int i = 0; i < m; i++)
    {
        u[i] = u[i+1];
    }
    for(int i = 0; i < n-1; i++)
    {
        y[i] = y[i+1];
    }
    y[n-1] = output;
    
    return output;
}

一応これでプロパーなパルス伝達関数をプログラムに落とし込むことができますが、ロボットをちょろっと動かす程度では色々効率が悪いコードなので、実際は以後の節で示すような個別の要素毎に実装することが多いです。

要素毎の実装例

積分器の実装

積分器のパルス伝達関数\frac{1}{z-1}となります。前節と同様に変形すると次のようになります。

\displaystyle{
y\slr{k}=x\slr{k-1}+y\slr{k-1}
}

C++風に落とし込むと次のようになります。

float integrator(float input)
{
    static float x = 0.0f;
    static float y = 0.0f;
    
    float res = x + y;
    
    y = res;
    x = input;
    
    return res;
}

ディジタルPID制御器の実装

積分器の実装で用いた考え方を応用します。簡単のため、制御則を次のようにします。

\displaystyle{
u\slr{k}=K_Pe\slr{k}+K_I\frac{z}{z-1}e\slr{k}+K_D\frac{z-1}{z}e\slr{k}
}

ただしe\slr{k}は偏差信号(目標値と出力値の差)とします。

積分項は前節で示した積分器と異なり、1サンプルの遅れがない形式での実装となっています。漸化式は次のようになります。

\displaystyle{
y_I\slr{k}=K_Ie\slr{k}+y_I\slr{k-1}
}

この形式ではパルス伝達関数\frac{K_Iz}{z-1}として見たときのものとなっていますが、次のようにすることで若干見通しが良くなります。

\displaystyle{
\eq{
y_I'\slr{k}&=e\slr{k}+y_I'\slr{k-1}\\
y_I\slr{k}&=K_Iy_I'\slr{k}
}
}

この形式ではパルス伝達関数\frac{z}{z-1}K_I倍している形となり同じ結果となります。

微分項は1サンプル前の偏差信号との差分となっています。よって、漸化式が次のようになります。

\displaystyle{
y_D\slr{k}=K_D\lr{e\slr{k}-e\slr{k-1}}
}

比例項は偏差信号の定数倍なので漸化式については割愛します。以上をまとめると、ディジタルPID制御器の漸化式が次のようになります。

\displaystyle{
\eq{
y_I'\slr{k}&=e\slr{k}+y_I'\slr{k-1}\\
u\slr{k}&=K_Pe\slr{k}+K_Iy_I'\slr{k}+K_D\lr{e\slr{k}-e\slr{k-1}}
}
}

C++風に落とし込むと次のようになります。

float PID(float input, float target)
{
    // input: 出力値
    // target: 目標値
    constexpr float Kp = 0.0f;   // Pゲイン
    constexpr float Ki = 0.0f;   // Iゲイン
    constexpr float Kd = 0.0f;   // Dゲイン

    static float error = 0.0f;   // 偏差信号 e[k]
    static float ierror = 0.0f;  // 積分項 y_I'[k]
    static float derror = 0.0f;  // 1サンプル前の偏差信号 e[k-1]
    static float output = 0.0f;  // 制御入力 u[k]

    error = target - input;     // 偏差信号=目標値-出力値
    ierror += error;            // 積分項の更新

    output = Kp*error + Ki*ierror + Kd*(error-derror);

    derror = error;

    return output;
}

参考まで、先に示したPID制御則の両辺にz-1を作用させることで速度型のPID制御則が得られ、積分項をキャンセルすることができます。一方で逆応答と呼ばれる望ましくない応答を示すことがあり、私は速度型ではない(位置型の)PID制御器を良く用いています。

まとめ

本稿ではパルス伝達関数を漸化式へ変形する方法について検討しました。また、C++風に実装した例も一緒に示しました。内容が基本的すぎるのか、手元にある制御理論の参考書には本稿のようなプログラムに落とし込む方法に関する記述が少ないように見受けられます。実装時に困っている誰かのお役に立てれば幸いです。