2016/03/28

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

マイコンで制御される被制御系(プラント)の安定を評価するには、プラントの標本化(離散化)が必要でした。
⇒ <参考>デジタル制御:最初の最初 ③ 

このデジタル制御におけるプラントの離散化は一般化して書くと
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.

2016/03/06

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

制御の勉強メモ…というブログの趣旨からは外れてしまいますが、
制御シミュレーションで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 

2016/03/05

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

制御周期の影響を変更確認するためのシミュレーションをしました。

下記パラメータでは制御周期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にあることが確認できました。