プラントモデルの離散化:Z変換表

Posted on 3/28/2016
このエントリーをはてなブックマークに追加
マイコンで制御される被制御系(プラント)の安定を評価するには、プラントの標本化(離散化)が必要でした。
⇒ <参考>デジタル制御:最初の最初 ③ 

このデジタル制御におけるプラントの離散化は一般化して書くとShalom D. Ruben先生の資料をお借りすると下記のように表せます。
Shalom D. Ruben
UCLA -Winter 2011 - MAE 171B Digital Control of Physical Systems
Chapter 4 - Discrete-Time System Representation
http://ecee.colorado.edu/shalom/DTRep.pdf
連続出力$ y(t) $のラプラス変換 $ Y(s) $は\[
Y(s) = G_{p}(s) \cdot G_{ZOH}(s) \cdot U^*(s)
\]
したがって、この出力をサンプリング(A/D変換)した結果をスター変換で表すと
(サンプリング周期 T[sec])\[\begin{align*}
Y^*(s) &= \frac{1}{2 \pi i} \lim_{T \to \infty}\int_{c - iT}^{c+iT}Y(\sigma)\frac{1}{1-e^{-T(s-\sigma)}} d\sigma \\
 &= \frac{1}{2 \pi i} \lim_{T \to \infty}\int_{c - iT}^{c+iT}G_{p}(\sigma) \cdot G_{ZOH}(\sigma) \cdot U^*(\sigma)\frac{1}{1-e^{-T(s-\sigma)}} d\sigma \\
 &= \Bigl(\frac{1}{2 \pi i} \lim_{T \to \infty}\int_{c - iT}^{c+iT}G_{p}(\sigma) \cdot G_{ZOH}(\sigma) \frac{1}{1-e^{-T(s-\sigma)}} d\sigma \Bigr) U^*(s)
\end{align*}\]
Note: 文献[1]の式(4.3)
となるから離散系にも伝達関数を定義できて、$ z = e^{sT} $ としてZ変換と紐づければ\[\begin{align*}
\mathcal{Z}[G(z)] \equiv \frac{Y(z)}{U(z)} &= \mathcal{Z}[G_{p}(s) \cdot G_{ZOH}(s)] \\
&= (1-z^{-1}) \mathcal{Z}[\frac{G_{p}(s)}{s}]
\end{align*}\]
Note: 文献[1]の式(4.10)

以下では、このD/A変換(ZOH)とプラントを直列にしたものを $ G(z) $の列とします。

Z変換表 : デジタル制御変換表
プラント伝達関数 $ G_p(s) $デジタル制御伝達関数 $ G(z) $
(サンプリング周期 T[sec])
積分系\[\frac{1}{s} \]\[\frac{T}{1-z^{-1}}\]
1次遅れ系
\[\frac{1}{\tau s+1} \]\[\frac{(1-e^{-\frac{T}{\tau}})z^{-1}}{1-e^{-\frac{T}{\tau}}z^{-1}} \]
\[\frac{\omega}{s+\omega} \]\[\frac{(1-e^{-\omega T})z^{-1}}{1-e^{-\omega T}z^{-1}} \]
2次遅れ系
$ 0 < \zeta < 1 $
\[\frac{\omega^2}{s^2+2\zeta \omega s + \omega^2}\]\[\frac{[1-\frac{e^{-\zeta \omega T}}{\sqrt{1-\zeta^2}}sin(\omega*\sqrt{1-\zeta^2}T+\phi)]z^{-1}+[e^{-2\zeta \omega T}+\frac{e^{-\zeta \omega T}}{\sqrt{1-\zeta^2}}sin(\omega \sqrt{1-\zeta^2}T-\phi)]z^{-2}}{1-2e^{-\zeta \omega T}cos(\omega \sqrt{1-\zeta^2} T)z^{-1}+e^{-2 \zeta \omega T}z^{-2}}\]ただし、$ \phi = acos(\zeta) $

Note 1:  高周波に対する安定の確認については別途検討が必要 http://ysserve.wakasato.jp/Lecture/ControlMecha3/node23.html

Note 2: 時間おくれについては拡張Z変換を使います
http://ysserve.wakasato.jp/Lecture/ControlMecha3/node16.html

【参考】

1. UCLA -Winter 2011 - MAE 171B Digital Control of Physical Systems
    Chapter 4 - Discrete-Time System Representation
    http://ecee.colorado.edu/shalom/DTRep.pdf

2. Neuman, C. P., and C. S. Baradello. "DIGITAL TRANSFER-FUNCTIONS FOR MICROCOMPUTER CONTROL." IEEE TRANSACTIONS ON SYSTEMS MAN AND CYBERNETICS 9.12 (1979): 856-860.

3. 株式会社デンソー, 制御装置. 特開2010-19105. 2010-1-28.

[Windows] C++からmatplotlibでグラフを表示する方法

Posted on 3/06/2016
このエントリーをはてなブックマークに追加
制御の勉強メモ…というブログの趣旨からは外れてしまいますが、
制御シミュレーションでPythonやC++を使うことが増え、
その結果を手軽にグラフ表示したくなることが多くなってきました.

シミュレーションの結果はCSVで書き出してMicrosoft Excelでグラフ化するなどしてますが、
もっと手軽にする方法を見つけたのでメモします.

会社などWindows以外を選べない環境でC++からmatplotlibを使う方法です.


環境


1. [準備] Libpython27.aの入れ替え

後述の工程でg++を使ってのコンパイル時に
libpython27.a: error adding symbols: File format not recognized
などというエラーが出てくるのでその対策をしておきます.
(最初はPython2.7のインストールで何かを間違えたかと思いましたが、
 良くあることらしくいくつか事例を見かけました.→ 事例① 事例② 事例③ )

以下は元はと言えば64bit版環境を構築する説明ですが、これを参考にlibpython27.aを置き換えました.
You can try to make libpython27.a or download from here (tip: you can unpack installer with 7-Zip). Make sure that you put libpython27.a in C:\Python27\libs directory.
Windows Installation Guide · Valloric/YouCompleteMe Wiki · GitHub 
ものぐさなので自分ではビルトしないことにしました
    libpython27.aをダウンロードしてしまう方法
  1. ここからlibpython-2.7.10-cp27-none-win32.whlをダウンロード
  2. 7Zipとして解凍
  3. 中に入っているlibpython27.aをC:\python27\libsへコピー (上書き)


2. [準備] matplotlib-cppの準備

Githubのリポジトリからmatplotlib-cpp-starterをダウンロードしてきます.
これを(Windows用に?)、少し改変します.
  • matplotlibcpp.h
    • 13行目
      #include <python2.7/Python.h>を#include <Python.h>に変更
  • main.cpp
    • #define M_PI 3.14159265などと追加

3. matplotlib-cppのサンプル実行

コマンドプロンプトでの実行用にバッチファイルを書きます.
(build_and_run.shの代わりになるファイル)

main.cppと同じフォルダにbuild_and_run.batファイルを作り、中に、
g++ ./main.cpp -lstdc++  -std=c++11 -I"C:\python27\include" -L"C:\Python27\libs" -lpython27 -o test && test
Pause
と書きます.

これを実行するとサンプルがコンパイルされて、以下が表示されるはずです.
(Atsuhiさんの見本と違う結果ですが2016/03/16現在は以下の結果になるmain.cppが添付されています)

4. 使い方


Gnuplotをパイプするより簡単で助かります!!


参考

[1] C++のコードから簡単にmatplotlibを使ってグラフを作成する方法 - MyEnigma 

[2] MinGWのインストールと使い方(2014年版) 

[3] Matplotlibの使い方 : matplotlib入門 - りんごがでている  http://bicycle1885.hatenablog.com/entry/2014/02/14/023734

[4] Windows Installation Guide · Valloric/YouCompleteMe Wiki · GitHub 

デジタル制御:最初の最初 ③

Posted on 3/05/2016
このエントリーをはてなブックマークに追加
制御周期の影響を変更確認するためのシミュレーションをしました。

下記パラメータでは制御周期70Hz前後が発散と持続振動の境目でしたが、この理由について確認してみます。

$ K_E $1 [V/rad/s]
$ \tau_m $1 [sec]
$ K_p $112 [V/rad/s]
$ K_i $3947 [V/rad]
(連続系であれば閉ループ伝達関数 $ \frac{2 \omega \zeta s + \omega ^2}{s^2 +2 \omega \zeta s + \omega ^2} $に対して $ \omega = 10[Hz] $, $ \zeta = 0.9 $となる条件)

デジタル制御
については、A/D変換およびD/A変換をより正確に書くと以下の通りになります。


答えからいうと、デジタル制御のお勉強: デジタル制御:最初の最初②では、
離散系を連続系に近似して制御系の評価をしようとしましたが、
サンプリング時間が長くなってくると(周期が下がると)、
離散系の特性が無視できなくなってきます。

これが、これまでの連続系を使った評価で「70Hz前後」を説明できない理由です。
今回は、逆に、連続系を離散系にサンプリングして制御系の評価をします。

まず、プラントであるDCモーターを離散表現するためにサンプリングについて考えます。

理想的な「サンプリング(A/D変換)」は入力に対し、くし型関数との掛け算で表されます。
入力を$ x(t) $ [連続値]とすると、くし型関数 $ \delta _T (t) $を使って、出力値$ x^*(t) $が $ x^*(t) = x(t) \cdot \delta_T(t) $と表されます。
関数の積に対するラプラス変換の関係は、
\[
\mathcal{L}[f(t) \cdot g(t)] = \frac{1}{2 \pi i} \lim_{T \to \infty}\int_{c - iT}^{c+iT}F(\sigma)G(s-\sigma) d\sigma
\]ですから、$ x(t) $のラプラス変換が $ X(s) $で、くし関数のラプラス変換が $ \mathcal{L}[\delta_T(t)]=\frac{1}{1-e^{-Ts}} $から
\[
X^*(s) = \frac{1}{2 \pi i} \lim_{T \to \infty}\int_{c - iT}^{c+iT}X(\sigma)\frac{1}{1-e^{-T(s-\sigma)}} d\sigma
\]となり、これがスター変換になります。これが離散系のラプラス変換にあたるZ変換につながるのですが
Yasunari SHIMADA先生のサイトがここを解説しています。
⇒ http://ysserve.wakasato.jp/Lecture/ControlMecha3/node5.html#SECTION00322300000000000000

また、D/A変換はデジタル制御:最初の最初①の通り0次ホールド(ZOH)で表現できます。

これらを合わせて制御系の離散系を求めるには、DCモータのZ変換(離散表現)を求める必要があり、
サンプリング周期 $ T_c $[sec]を使って ( $z = e^{s T_c} $)
\[\begin{align*}
\mathcal{Z}[ZOH(s) \cdot \frac{1/K_E}{s\tau_m+1}]&=(1-z^{-1})\mathcal{Z}[\frac{1}{s}\frac{1/K_E}{s\tau_m+1}] \\
&=\frac{1}{K_E}\frac{(1-e^{-\frac{T_c}{\tau_m}})z^{-1}}{1-e^{-\frac{T_c}{\tau_m}}z^{-1}}
\end{align*}
\]で求められます。
(幸い1979年のC. P. NEUMAN氏によるDigital Transfer Functions for Microcomputer Controlがオープンアクセスになっているので、こちらを参照すればデジタル制御で使うプラントの離散化には困らないでしょう)

この結果を用いて、離散系の一巡伝達関数を求めると
\[
G(z) \equiv (K_p + K_i T_c \frac{1}{1-z^{-1}})(\frac{1}{K_E}\frac{(1-e^{-\frac{T_c}{\tau_m}})z^{-1}}{1-e^{-\frac{T_c}{\tau_m}}z^{-1}})
\]
よって、この離散系の特性方程式は、
\[
1+G(z) = 0
\]となり、特性根を $ \gamma_1 \cdots \gamma_n $とすると、制御系の閉ループ伝達関数は定数$ C_1 \cdots C_n $を用いて
\[
\frac{G(z)}{1+G(z)}  = \frac{C_1}{z-\gamma_1}+\cdots+\frac{C_n}{z-\gamma_n}
\]と表現できます。

したがって、
全ての極に対して$ | \gamma | < 1 $であるとき、系は収束
どれか一つでも1に等しくなると、振動系
1より大きくなると発散してしまいます。

先にシミュレーションした条件では70Hz前後が、この極が1より大きくなる閾値だったということです。

確かめてみましょう。
\[
1+G(z) = 0 \leftrightarrow z^2+((1-e^{-\frac{T_c}{\tau_m}})(\frac{K_p}{K_E}+\frac{K_i T_c}{K_E})-(e^{-\frac{T_c}{\tau_m}}+1))z+(-(1-e^{-\frac{T_c}{\tau_m}})\frac{K_p}{K_E}+e^{-\frac{T_c}{\tau_m}})=0
\]
この特性方程式に不安定根が含まれるかどうかは$ z = \frac{1+w}{1-w} $という変換をしたうえで、
Routh Hurwitzの判別法が使えます。
2次のときは後述資料の通り下記の不等式を満たすことが条件です。

  • $ 2(1+e^{-\frac{T_c}{\tau_m}})-(1-e^{-\frac{T_c}{\tau_m}})(2\frac{K_p}{K_E}+\frac{K_i T_c}{K_E}) > 0$
  • $ (1-e^{-\frac{T_c}{\tau_m}})(\frac{K_p}{K_E}+1) > 0$
  • $ 2+(1-e^{-\frac{T_c}{\tau_m}})\frac{K_i T_c}{K_E} > 0$
  • $ 1+e^{-\frac{T_c}{\tau_m}}-(1-e^{-\frac{T_c}{\tau_m}})\frac{K_p}{K_E}>0$

http://wolfweb.unr.edu/~fadali/ee472/Stability.pdf

第2式、第3式は $ K_p, K_i, T_c > 0 $から恒等的に成立します。
第1式、第4式については前記のパラメータを使って、$ 2(1+e^{-T_c})-(1-e^{-T_c})(224+3947 T_c) > 0$および$ 1+e^{-T_c}-(1-e^{-T_c})112>0$から$ 0 < T_c < 0.0142696 $が求まります。
(数値解: Wolfram Alpha)

使用したパラメータ(再掲)
$ K_E $1 [V/rad/s]
$ \tau_m $1 [sec]
$ K_p $112 [V/rad/s]
$ K_i $3947 [V/rad]

制御周期 0.0142696 secのとき、制御周波数は70.079..[Hz]ですから、
持続振動と発散(不安定)の境界が約70Hzにあることが確認できました。

デジタル制御:最初の最初 ②

Posted on 3/04/2016
このエントリーをはてなブックマークに追加
デジタル制御:最初の最初①で、デジタル制御を連続系で近似する方法があることを知りました.
今回は近似の効能を確認するところをスタートにしてみます.

制御系のモデルとしては、ここまでで3つあります.

  1. 連続系:オリジナル
  2. デジタル制御
  3. ②を連続系に近似したもの


ここで、$ K_E = 1 $, $ \tau_m = 1 $ とします.
さらに、PIコントローラの係数として$ K_p = 112 $, $ K_i = 3947 $を採用して、
①の閉ループ伝達関数 $ \frac{2\omega \zeta s + \omega^2}{s^2+2\omega \zeta s + \omega^2} $が、$ \omega = 10Hz $, $ \zeta = 0.9 $程度になるように調整しました.
センサのサンプリング周期 $ T_s $は制御周期 $ T_c $の1/100未満として、
サンプリングによるむだ時間などは無視できると仮定しましょう.


確認のため、連続系(オリジナル)のステップ応答を確認すると以下の通りです.
次に、制御周期を1kHz ( $ T_c = 1msec $)として、計算してみます.
閉ループ系の周波数帯域が10Hz程度であるのに対して制御系が1kHzあるので、
デジタル化による問題は生じてないように見えます.

まだまだ制御周期を落としてみましょう.500Hzと250Hzです.
250Hzで「連続系近似」が「デジタル制御」を近似できなくなってきているように見えます.
さらに制御周期を落として、125Hzと62.5Hzです.62.5Hzに至っては遂に発振してしましました.
この発振するかしないかの閾値は、この系の場合だと70Hz前後にあり、
そこ周辺の制御周波数でシミュレーションしてみると
75Hzで収束し、70.5Hzで持続的な振動が起こる…という現象が見られます.

この辺では「連続系近似」はもはや用をなしていないですね.

この70Hzという数字、どこから来たのでしょう.

つづく⇒ デジタル制御:最初の最初 ③ 

デジタル制御:最初の最初 ①

Posted on 2/25/2016
このエントリーをはてなブックマークに追加
3部構成の①です。

DCモータの速度制御を例にデジタル制御を考えみたい.

DCモータPI制御については、以前、連続系PID制御 : DCモータの速度制御で少し取り扱ったが、
多摩川精機さんのカタログ P.5にならって、簡略化したDCモータの伝達関数を使うことにしました[1].

DCモータの端子電圧V [V]に対するモータ回転数ω [rad/s]の伝達関数は
\[ G(s) = \frac{\omega (s)}{V(s)}=\frac{1/K_E}{s \tau_m+1} \]
ここで
$ K_E $ : 誘起電圧定数 $ [V/\frac{rad}{s}] $
$ \tau_m $ : 機械的時定数 $ [sec] $ 

このモータを連続系のPI制御で速度制御する場合のブロック図は、
これを例えばArduinoのようなものでデジタル制御する場合、
モータをPWM制御することが前提となって、実装上、だいたい以下のような形になります.


「A/D変換(デジタル変換)」「D/A変換(アナログ変換)」「デジタルPI」が連続系と異なるところで、
センサの値をマイコンなんかに取り込むために必須の処理です.

このデジタル変換で、「量子化」(値の離散化)と「標本化」(時間の離散化)をいっぺんにやってしまうのですが [2]、
中でも「標本化」が制御の安定性に大きな影響を与えます.

1. 標本化 と サンプル&ホールド

「標本化」「量子化」の具体的な仕組みは[Wikipedia : A/D変換回路]を参考にするとして、
標本化(サンプリング)だけに注目すると下図のような仕組みでA/D変換がなされる.

[標本化の流れ]
① 連続時間の信号が入力される
② サンプラ部がサンプリング周期T[sec]毎に値を切り出す
③ ホールド部が切り出した値をT[sec]間維持する→出力へ

この③の「切り出した値をT[sec]間維持する」という部分だけを取り出すと、
0次ホールド(Zero Order Hold : ZOH) $ H_{ZOH}(s) = \frac{1-e^{-sT}}{sT} $という伝達関数で表現できるのですが [3]
Matlabなどのシミュレータは、②と③の機能を合わせて"ZOH"と呼んでいるようです[4].
(便宜上でしょうか?)
このZOHは離散系(デジタル)の制御を、連続系として扱って考えるために便利なようです.

2. デジタル制御設計を連続系で…

ここまでを まとめて、ブロック図を書き直すと


制御周期$ T_c $とセンサのサンプリング周期 $ T_s $については、大抵 $ T_c > T_s $なので、 
特にデジタルフィルタ(含 移動平均)などの処理をしていなければ、
A/D変換のブロックは 1 に置き換えてもよいと思います.

すると、Yasunari SHIMADA先生のページ「[5]  連続制御系で近似する方法 」にならって、
連続系に近似的に表現できて、以下の制御ブロックとして扱えます.


背景にはスター変換(ラプラス変換の離散版)の考え方がありますが、
まずは、この近似の効果を見てみることにします.


参考文献
[1] DCサーボモータ - サーボモータ|多摩川精機株式会社 - カタログ (PDFへのリンクあり)
http://www.tamagawa-seiki.co.jp/jpn/servo/tredc.html

[2] CH:2 Discretization of continuous signals
http://www.wakayama-u.ac.jp/~kawahara/signalproc/CH2cont2discrete/contsig2discrete.html

[3] Wikipedia : ZOH
https://en.wikipedia.org/wiki/Zero-order_hold

[4] 連続/離散の変換方法 - MATLAB & Simulink - MathWorks 日本
http://jp.mathworks.com/help/control/ug/continuous-discrete-conversion-methods.html#bs78nig-2

[5] 連続制御系で近似する方法
http://ysserve.wakasato.jp/Lecture/ControlMecha3/node27.html

移動平均とローパスフィルタ

Posted on 1/20/2016
このエントリーをはてなブックマークに追加
デジタル制御ではセンサのAD変換値を平滑化するために移動平均を使います.
この移動平均を連続系で考えてみるために考察します.

サンプリング周波数$ f_s $[Hz], 移動平均数 N個の移動平均は、ローパスフィルタと同じような役割をします.

カットオフ周波数$ f_c $ は、おおよそ下記で表されます.
\[
f_c = \frac{0.443}{\sqrt{N^2 - 1}} \cdot f_s
\]

公式の理屈は参考文献[1]を参照していただいて、
まずステップ応答を比べてみます.

計算条件:カットオフ周波数44.3Hzで合わせたフィルタ
① サンプリング周波数 1kHz, 10個の移動平均
② サンプリング周波数 2kHz, 20個の移動平均
③ カットオフ周波数 44.3Hzの1次ローパスフィルタ $ \frac{2 \pi \cdot 44.3}{s+(2  \pi \cdot 44.3)} $
④ カットオフ周波数 44.3Hz、減衰比0.7の2次ローパスフィルタ $ \frac{(2 \pi \cdot 44.3)^2}{s^2+2 \cdot 0.7 \cdot (2 \pi \cdot 44.3)s+(2 \pi \cdot 44.3)^2} $

応答としては2次のローパスフィルタに近いと考えたらよいのでしょうか?


次はボード線図で比べてみます.Scilabを使って①、②、③、④を比較した図です.
(考え方は[2]を参考にしました)

ボード線図の見どころは

a. 先にあげた公式はゲインが-3.0dBになる周波数を求めている
  (4つの曲線はすべて44.3Hz, -3.0dBを通る)

b. 移動平均のフィルタとしての特性は、移動平均数 ÷ サンプリング周波数で決まる.

c. 0~100Hzの範囲での移動平均は1次のローパスフィルタより強く、2次のローパスフィルタに似ている
      (先のステップ応答で見られた通り)

d. 100~500Hzまで見てみると、ところどころ減衰が悪い周波数があり、
 その周波数では1次LPF程度のレベルで考えたほうが良い

というところだと思います.



参考文献
[1] 移動平均の周波数特性 - メモがわり http://seesaawiki.jp/w/weidows95/d/%B0%DC%C6%B0%CA%BF%B6%D1%A4%CE%BC%FE%C7%C8%BF%F4%C6%C3%C0%AD

[2]デジタルフィルタを計算する
http://bluefish.orz.hm/sdoc/digifil2.html#移動平均フィルタ
Copyright © Since 2014 Chachay All Rights Reserved. Powered by Blogger.