DigitalServo

Blog

About

Library

    線形制御系設計

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

    モーションコントロール

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

    モータドライブ

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

    システム同定

    ホワイトボックスモデル

双一次変換を用いたシステムの実装

制御対象を連続時間系でモデル化して制御器を設計した場合には、制御器が連続時間系の表現で記述されるため計算機に直接実装することができません。そのため、離散化処理を介して実装可能な形式で表現する必要があります。双一次変換は実現したい連続時間系の入出力伝達関数から離散時間系の入出力伝達関数を導出する際に使用できます。

双一次変換

ラプラス変換によって時間遅れ要素は以下のように線形作用素として表現されました。
L[x(t)]≡X(s)⇒L[x(t−L)]=e−sLX(s)\begin{align} \mathcal{L}[x(t)] \equiv X(s) \Rightarrow \mathcal{L}[x(t-L)] = e^{-sL}X(s) \end{align}
ここでサンプリング間隔を TsT_{\rm s} 、信号x(t)x(t)の nn ステップ分の時間が遅れた信号を xn(t)x_n(t) とすると、そのラプラス変換は次のように表現されます。
xn(t)≡x(t−nTs)⇔L[xn(t)]=e−snTsX(s)\begin{align} x_n(t) &\equiv x(t-nT_{\rm s})\\ \Leftrightarrow \mathcal{L}[x_{\rm n}(t)] &= e^{-snT_{\rm s}}X(s) \end{align}
計算機で扱える信号はサンプルされた信号のみとなり、離散信号のみから所望の伝達特性を実現しなければなりません。そこで、シフトオペレータ z≡esTsz \equiv e^{sT_{\rm s}} を導入し、伝達関数を zz を用いて表現することを目指します。zz を用いて伝達関数を表現することは時間遅れ要素のみを用いて伝達関数を表現することであり、過去のサンプルを使用して所望の伝達特性を近似的に実現可能であることを示しています。シフトオペレータの定義より、ss は zz を使用して次のように表現されます。
s=1Tslog⁡(z)=2Tstanh−1(z−1z+1)=2Ts∑n=0∞12n+1(z−1z+1)2n+1\begin{align} s &= \frac{1}{T_{\rm s}}\log(z) \\ &= \frac{2}{T_{\rm s}}{\rm tanh}^{-1}\left( \frac{z - 1}{z + 1} \right)\\ &= \frac{2}{T_{\rm s}}\sum ^{\infty}_{n=0}\frac{1}{2n+1}\left( \frac{z - 1}{z + 1} \right)^{2n+1} \end{align}
級数展開を一次で打ち切った場合の結果が双一次変換となります。
s=2Tsz−1z+1+2Ts∑n=1∞12n+1(z−1z+1)2n+1≈2Tsz−1z+1\begin{align} s &= \frac{2}{T_{\rm s}}\frac{z - 1}{z + 1} + \frac{2}{T_{\rm s}}\sum ^{\infty}_{n=1}\frac{1}{2n+1}\left( \frac{z - 1}{z + 1} \right)^{2n+1} \\ &\approx \frac{2}{T_{\rm s}}\frac{z - 1}{z + 1} \end{align}
ここで、実現したい伝達関数 Gc(s)G_{\rm c}(s) に対して離散化伝達関数 Gd(z)G_{\rm d}(z) は以下のように導出されます。
Gd(z)=Gc(2Tsz−1z+1)\begin{align} G_{\rm d}(z) &= G_{\rm c}\left (\frac{2}{T_{\rm s}}\frac{z - 1}{z + 1} \right) \end{align}
計算機で実装する場合には過去の信号のみが使用可能であることを考慮して、次の形式で導出することもできます。
Gd(z)=Gc(2Ts1−z−11+z−1)\begin{align} G_{\rm d}(z) &= G_{\rm c}\left (\frac{2}{T_{\rm s}}\frac{1 - z^{-1}}{1 + z^{-1}} \right) \end{align}

離散化伝達関数の実装方法

離散化伝達関数は実現したい伝達関数を時間遅れ要素を用いて表現したものとなります。シフトオペレータの定義より以下の式が成立することを利用して、実装方法を確認します。
L−1[z−nX(s)]=xn(t)\begin{align} \mathcal{L}^{-1}\left[z^{-n}X(s)\right] &= x_{n}(t) \end{align}
入出力伝達特性が次のように表現される場合について考えます。
Y(s)=Gd(z)X(s)Gd(z)=b0+b1z−1+⋯+bpz−p1+a1z−1+⋯+aqz−q\begin{align} Y(s) &= G_{\rm d}(z)X(s)\\ G_{\rm d}(z) &= \frac{b_{0}+b_{1}z^{-1}+\cdots+b_{p}z^{-p}}{1+a_{1}z^{-1}+\cdots+a_{q}z^{-q}} \\ \end{align}
この式を整理することで、次の関係式を得ます。
(1+∑i=1qaiz−i)Y(s)=∑i=0pbiz−iX(s)∴ Y(s)=∑i=0pbiz−iX(s)−∑i=1qaiz−iY(s)\begin{align} &\left( 1 + \sum^{q}_{i=1}a_iz^{-i} \right) Y(s) = \sum^{p}_{i=0}b_iz^{-i} X(s) \\ \therefore\ & Y(s) = \sum^{p}_{i=0}b_iz^{-i} X(s) - \sum^{q}_{i=1}a_iz^{-i} Y(s) \end{align}
ここで、逆ラプラス変換を用いることで時系列の関係が導出されます。
y(t)=∑i=0pbixi(t)−∑i=1qaiyi(t)\begin{align} y(t) = \sum^{p}_{i=0}b_i x_{i}(t) - \sum^{q}_{i=1}a_i y_{i}(t) \end{align}
計算機上では上記の式を離散化して実装する必要があります。任意のサンプル点 kk を基点として、その点におけるサンプル値を x[k]x[k] および y[k]y[k]、nn ステップ前のサンプル値を x[k−n]x[k-n] および y[k−n]y[k-n] と表現すれば、サンプル値には次の関係が成立します。
y[k]=∑i=0pbix[k−i]−∑i=1qaiy[k−i]\begin{align} y[k] = \sum^{p}_{i=0}b_i x[k-i] - \sum^{q}_{i=1}a_i y[k-i] \end{align}

実装例

ここではシンプルな低域通過フィルタの設計例を示します。以下の低域通過特性を示す伝達関数を考えます。
Gc(s)=gs+g\begin{align} G_{\rm c}(s) &= \frac{g}{s + g} \end{align}
ここで、g∈R>0g\in \mathbb{R} > 0 は遮断周波数 [rad/s][{\rm rad/s}] を表します。離散化伝達関数は次のように表現されます。
Gd(z)=Gc(2Ts1−z−11+z−1)=gTs(1+z−1)(2+gTs)+(−2+gTs)z−1\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)z−1)Y(s)=gTs(1+z−1)X(s)⇔ Y(s)=gTs(1+z−1)2+gTsX(s)+2−gTs2+gTsz−1Y(s)⇔ y[k]=gTs2+gTs(x[k]+x[k−1])+2−gTs2+gTsy[k−1]\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}
以下に設計した低域通過フィルタのクラスを例示します。
 lib/LowPassFilter.h 
class LowPassFilter {
  private:
    double Ts;
    double g;
    double xZ1;
    double yZ1;
    double coeff[2];

  public:
    LowPassFilter(double gn, double SmplTime);
    ~LowPassFilter();
    double Output(double x);
    void ClearStateVars(void);
};
 lib/LowPassFilter.cpp 
#include "LowPassFilter.h"

LowPassFilter::LowPassFilter(double gn, double SmplTime): Ts(SmplTime), g(gn) {
  xZ1 = 0;
  yZ1 = 0;
  coeff[0] = g * Ts / (2.0 + g * Ts);
  coeff[1] = (2.0 - g * Ts) / (2.0 + g * Ts);
}

LowPassFilter::~LowPassFilter() {
}

double LowPassFilter::Output(double x) {
  double y = coeff[0] * (x + xZ1) + coeff[1] * yZ1;
  xZ1 = x;
  yZ1 = y;
  return y;
}

void LowPassFilter::ClearStateVars(void) {
  xZ1 = 0;
  yZ1 = 0;
}
このクラスの使用例は以下のようになります。
 main.cpp 
#include "LowPassFilter.h"

int main (void) {

  const double Ts = 1e-3; // Sampling time
  const double gn = 20.0; // Cutoff frequency

  double t = 0.0;
  double x = 0.0;
  double y = 0.0;

  // Instance
  static LowPassFilter* LPF = new LowPassFilter(gn, Ts);

  for (int i = 0; i < 5e4; i++) {
    x = (t > 1.0)? 1.0 : 0.0;
    y = LPF -> Output(x);
    t += Ts;
  }

  delete LPF;

  return 0;
}

CONTENTS

© DigitalServo