プラントモデルの離散化: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{aligned} 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{aligned}\]
Note: 文献[1]の式(4.3)
となるから離散系にも伝達関数を定義できて、$ z = e^{sT} $ としてZ変換と紐づければ\[\begin{aligned} \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{aligned}\]
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でグラフを表示する方法

制御シミュレーションでPythonやC++を使うことが増え、その結果を手軽にグラフ表示したくなることが多くなってきました.

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

この方法はAtsushiさんの[C++のコードから簡単にmatplotlibを使ってグラフを作成する方法]を参考に会社などWindowsしか選べない環境でC++からmatplotlibを使う方法です.

環境

  • Windows 10 - 64bit (Windows 7など)
  • Anaconda (Python 2.7, 3.5, 3.6いずれでも)
  • Visual Studio 2013, 2015, 2017 いずれでも

1. Visual Studio 2017の場合

 環境変数VS150COMNTOOLSを設定。Visual Studio Build Toolsとしてデフォルトのパスにインストールしている場合は、"C:\Program Files (x86)\Microsoft Visual Studio\2017\BuildTools\Common7\Tools\"

その他はご自分の環境に合わせて. Visual Studio 2013, 2015の場合はデフォで設定済みになっています。

2. Matplotlib-cppを取得し、設定

Chachay/matplotlib-cppからクローンし、Contrib/WinBuild.cmdをエディタで開き、環境に合わせて設定。


if NOT DEFINED MSVC_VERSION set MSVC_VERSION=12または14,15
if NOT DEFINED CMAKE_CONFIG set CMAKE_CONFIG=Release
if NOT DEFINED PYTHONHOME set PYTHONHOME=Pythonフォルダへのリンク

3. ビルドして実行

コマンドプロンプトやPowerShellでcontribフォルダを開き、Contrib/WinBuild.cmdを実行. examples/build/releaseにサンプルファイルがリリースされているはず. 


4. 自分のファイルは?

cmakelist.txtを参考にどうぞ!

環境


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 

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

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

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

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

デジタル制御:最初の最初①で、デジタル制御を連続系で近似する方法があることを知りました.
今回は近似の効能を確認するところをスタートにしてみます.

制御系のモデルとしては、ここまでで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という数字、どこから来たのでしょう.

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