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にあることが確認できました。

2016/03/04

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

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

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

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

2016/02/25

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

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

2016/01/20

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

デジタル制御ではセンサの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#移動平均フィルタ

2015/11/29

むだ時間を持つ1次遅れ系のPI制御パラメータの解析的な分析

プラントにむだ時間を持つ系の制御が難しいことを、
Ziegler Nicholsの限界感度法を試すなかで確認しました.

アクチュエータにむだ時間を持つ1次遅れ系をPI制御する場合について、
$ K_p, K_i $の取りうるパラメータの範囲について考えてみます.

この系の特性方程式は、
\[
(K_p + \frac{K_i}{s})(\frac{e^{-L \cdot s}}{\tau \cdot s + 1})+1=0
\]
この特性方程式からラウス・フルビッツの安定判別法 (Routh–Hurwitz stability criterion)などから方程式を解かずに安定性を判定できれば良いのだけれど、この方程式は多項式ではなく、つかえません.

しかし、あらためて閉ループ伝達関数を求め、
\[\begin{align*}
G(s) &= \frac{P \cdot C}{1+P \cdot C} \\ &= \frac{(K_p s + K_i)e^{-L \cdot s}}{\tau \cdot s^2 + s +(K_p s + K_i)e^{-L \cdot s}}
\end{align*}\]
もし、この伝達関数を部分分数分解できて、ヘヴィサイドの展開定理のように、
\[
G(s) = \sum_{i=1}^{\infty} \sum_{j=1}^{\infty}  \frac{A_{ij}}{(s-a_{i})^j}
\]
などと表せれば、有理関数と同じように考えることができそうです [1].
すると極が全て複素平面の左半面にあれば系は安定と言えます.
⇒ すべての極が複素平面の左半面にある条件を求め、それを元に$ K_p, K_i $の範囲を求める段取りとする.


1. エルミート・ビーラー(Hermite-Biehler)の定理

伝達関数 $ G(s) $の極が全て複素平面の左半面にある条件を求める部分だけが難しそうです.
これについてはエルミート・ビーラーの定理を使って考えることができます[2][3].

先の特性方程式の左辺を変形し、以下のように $ \delta(s) $を定義する.
\[
\delta(s) = \tau s^2 + s + (K_p \cdot s + K_i) e^{-Ls}
\]
ここで、
\[\begin{align*}
\delta^*(s) &= e^{Ls} \delta(s) \\ &= (K_p \cdot s + K_i) + (\tau s^2 + s) e^{Ls}
\end{align*}\]
とし、$ s = \omega i $ を代入して実部 $ \delta_r $と虚部$ \delta_i $で関数を表現すると
\[\begin{align*}
\delta^* (\omega) &= \delta_r(s) + i \delta_i(s) \\ &=(K_p \omega i + K_i) + (-\tau \omega^2 + \omega i) e^{L\omega i}
\end{align*}\]
\[
\delta_r(\omega) = K_i - \tau \omega^2 cos(L \omega) - \omega sin (L \omega)
\]
\[
\delta_i(\omega) =  \omega (K_p - \tau \omega sin(L \omega) + cos(L \omega))
\]

ここでエルミート・ビーラー(Hermite-Biehler)の定理から、

1. $ \delta_r(s) $の解と$\delta_i(s)$の解が、すべて実数で、重解はなく、互いに隔離している.

【注:解の隔離】
方程式 f(z) = 0 の解と方程式 g(z) = 0 の解が互いに解を隔離するとは、どちらの方程式も重解を持たず、大小の順において f(z) = 0 の隣接する 2 つの解の間に g(z) = 0 の1 つの解があり、また g(z) = 0 の隣接する 2 つの解の間に f(z) = 0 の 1 つの解があることを言う。

[3] あってよかった複素数
http://izumi-math.jp/F_Yasuda/complex_number/good.pdf

さらに
2. $ z=L \omega$として、
\[
\frac{d\delta_i}{dz}\delta_r - \delta_i \frac{d\delta_r}{dz}>0
\]

以上の1, 2を満たす範囲でG(s)の解は複素平面の左半面にあると言え、G(s)は安定.

これらから $ K_p, K_i $の関係を一般的に解くの私にはできませんでした.
(なので数値解析でお試してみます)


2. 定理を試してみる

実際に$ \tau = 1 $ $ L = 1 $の条件で(つまり$G(s) = \frac{1}{s+1}e^{-s}$)、
安定になる$ K_p, K_i $の範囲(境界線)を数値解析で求めたのが、こちらの範囲.

PI制御のチューニング結果を比較するため、ジーグラ・ニコルスの限界感度法(ZN限界感度法)やジーグラ・ニコルスのステップ応答法(ZNステップ応答法)で求まる結果もプロットしました.

両結果とも安定領域にきっちり入っていますね.
計算した条件ではステップ応答法の方がおとなしいチューニング結果になっていそうです.


もう少し掘り下げて勉強してみます.

参考文献

[1] ヘヴィサイドの展開定理
ヘヴィサイドの展開定理は有理型関数にも拡張できるとのことなので、大丈夫だと思います.

[2] Generalizations of the Hermite–Biehler theorem http://www.sciencedirect.com/science/article/pii/S0024379599000695

[3] PID Controllers for Systems with Time-Delay
http://msc.berkeley.edu/PID/modernPID3-delay.pdf

[4] あってよかった複素数
http://izumi-math.jp/F_Yasuda/complex_number/good.pdf

2015/07/12

PI制御チューニング : Ziegler Nicholsの限界感度法 - むだ時間ありの1次遅れ系への適用

PI制御チューニング : Ziegler Nicholsの限界感度法の続きです.

0. むだ時間 - 導入

システムの入力に対する出力の時間遅れの1種で、時間領域で以下のように表現でき、
\[
y(t) = x(t - T_{d})
\]
伝達関数では、
\[
G(s) = e^{-sT_{d}}
\]
と表される.

このような性質が付与されたシステムを制御する時、
制御コントローラにしてみれば、制御してから結果に反映されるまで時間がかかる状態になるので、
ゲインを高く設定しすぎるとハンチングを起こすなどそわそわしている状態になります。

こういった性質は制御分野以外でも最近とりあげられていて、
「世界はシステムで動く」というドネラさんの本でも、
例えば在庫管理的な場面でむだ時間が生じると現象が振動的になるという話が出ていました。
(制御工学で社会分析している印象を受ける本です. 応用としては面白いかも.)



むだ時間は実世界では必ず存在し、その影響が制御設計やチューニングを難しくします


1. 1次遅れ系とむだ時間

むだ時間が付与された1次遅れ系のベクトル軌跡を見てみる.
1次遅れ系の時定数を時定数1(sec), むだ時間を1(sec)として伝達関数は、
\[G_2(s)=\frac{1}{s+1}e^{-s}\]

これを複素平面で表現するため極形式にすると、
\[G_2(\omega i) = \frac{1}{\omega i + 1}e^{-\omega i }\]
から
\[
z=|G_2(\omega i) | \cdot exp( i \cdot arg(G_2(\omega i))) = \frac{1}{\sqrt{\omega^2+1}} exp(i(-arctan(\omega)+(-\omega)))
\]

これを図で表すと分かる通り、時間遅れのない状態(水色)に時間遅れをつけるとオレンジの線のようにぐるぐるができ、不安定点(Re=-1)に近づき、より不安定になっていることがわかります.
水色($ \frac{1}{s+1}$)、オレンジ($\frac{1}{s+1}e^{-s} $)
Note: お絵かきにはエクセルを使い、ω = [0.01, 100] (rad/s)の範囲で計算しました.

$ \omega $r$ \theta $RealIm
0.01 ~ 100
= 1/sqrt(A2^2+1)=- atan2(1, A2) - A2=B2 * cos(C2)=B2 * sin(C2)


2. むだ時間1次遅れ系のPI制御へジーグラ・ニコルス限界感度法を適用

先の$ P(s) = \frac{1}{\tau s+1} $ として、プラント $ P_2(s) =  P(s) e^{-Ls} $ にPI制御 $ C(s) = K_p + \frac{K_i}{s} $ を施すと、
\[
 G(s) = \frac{P \cdot C \cdot e^{-Ls}}{1+P \cdot C \cdot e^{-Ls}} = \frac{K_p e^{-Ls} s + K_i e^{-Ls}}{\tau s^2 + (1+K_p e^{-Ls})s + K_i e^{-Ls}}
\]

$ K_u $を求めるために、$ K_i = 0 $とすると、
\[
 G(s) = \frac{K_p e^{-Ls}}{\tau s + 1+K_p e^{-Ls}}
\]

これに対し、当初、$ K_u $を求めるために$ e^{-Ls}$のパデ近似を用いて代数的に範囲を求めようと考えたのですが、文献[2]で検証している通り、特に1次,2次などの低次のパデ近似では良い解にたどりつけないようです.

$ K_u $を求めるにあたってナイキスト線図を使って幾何的に考えることにします.

まず、$ K_p = 1 $として
\[
 G(s) = \frac{e^{-Ls}}{\tau s + 1+K_p e^{-Ls}}
\]
のナイキスト線図を描きます.

この時
\[ z=|G(\omega i) | \cdot exp( i \cdot arg(G(\omega i))) = \frac{1}{\sqrt{(\tau \omega)^2+1}} exp(i(-arctan(\tau \omega)+(-L \omega))) \]
を描くことになるが、以下図中の緑の点の座標を求めます.


$ \omega $が大きくなるほどゼロ点に近づくので、
zが複素平面左半面の実軸上と交わる正で最小の$ \omega_0 $を求めればよい.

従って 、$ -arctan(\tau \omega_0) + (-L \omega_0) = -\pi $で、最小のものが求める周波数です.
※ この時、$ z(\omega_0) = - \frac{1}{\sqrt{(\tau \omega_0)^2 +1}} $

ここで $ K_p = \sqrt{(\tau \omega_0)^2 +1} $とするとベクトル線図は以下の図の水色の線のようになり、これが $ K_u $となります.


拡大した図
まとめ:
以上の検討から$ K_p $の限界である$ K_u $は、
\[
K_u = \sqrt{(\tau \omega_0)^2 +1}
\]
ただし、
\[
arctan(\tau \omega_0)  = \pi - (L \omega_0)
\]
($ \tau \cdot \omega_0 $は、(0,$\pi/2$)の範囲 )

また、この時、振動周期(s)は$T_u = \frac{2\pi}{\omega_0} $である.

$ \omega_0 $は数値解析で求めるほかありませんが、二分法など使うと良いでしょう.

二分法 - 高精度計算サイト http://keisan.casio.jp/exec/system/1180917668

2.1. お試し計算

例えば、\[G(s)=\frac{1}{s+1}e^{-s}\]として$ K_u $を求めると、$ K_u \approx  2.26 (\omega_0 \approx 2.029; T_u \approx 3.097sec) $と求まります.

シミュレーションで、この値を使ってステップ応答を求めると以下の結果になります.


ジーグラ・ニコルスの限界感度法に従って、
制御の種類$ K_p $$ K_d $$ K_i $
P制御$ 0.5 K_u $00
PI制御$ 0.45 K_u $0$ \frac{0.45 K_u}{0.83 T_u} $
PID制御$ 0.6 K_u $$ 0.075 T_u K_u $$ \frac{0.6 K_u}{0.5 T_u} $


から$ K_p = 1.018 $, $ K_i = 0.396 $としてステップ応答を求めたものがこちら.


やっとZNの限界感度法がうまく行きました.

参考

[1] The Time Delay
http://lpsa.swarthmore.edu/BackGround/TimeDelay/TimeDelay.html

[2] PID Controllers for Systems with Time-Delay
http://msc.berkeley.edu/PID/modernPID3-delay.pdf


2015/07/04

PI制御チューニング : Ziegler Nicholsの限界感度法 - 2次遅れ系への適用 (むだ時間なし)

PI制御チューニング : Ziegler Nicholsの限界感度法の続きです.

2. ジーグラ・ニコルスの限界感度法と2次遅れ系

同様に、むだ時間のない2次遅れ系にPI制御を適用することを考える.
まず、
\[ P_1(s) = \frac{\omega_n^2}{s^2+2 \zeta \omega_n s+\omega_n^2} \]
\[ C_1(s) = K_p + \frac{K_i}{s} \]
従って、フィードバック系の伝達関数は、
\[ G_1(s) = \frac{P_1 C_1}{1+P_1 C_1} = \frac{\omega_n^2 K_p s + \omega_n^2 K_i}{s^3+2 \zeta \omega_n s^2+\omega_n^2(1 +K_p)s+\omega_n^2 K_i} \]

0.-1の段取りに従い $ K_i = 0 $として式を整理すると2次遅れ系の形式になる.
\[ G_1(s) = \frac{\omega_n^2 K_p}{s^2+2 \zeta \omega_n s+\omega_n^2 (1+K_p)} \]

ここで新たにP制御を加味した減衰比が定義でき、
\[ \zeta_{1} = \frac{\zeta}{\sqrt{1 + K_p}} \]

従って$ \zeta_{1}=0 $であれば振動的になると言えるのですが([5]最下部)、これが起こるのは $ K_p $ → ∞ の時だけ.
便宜上、 $ \zeta_{1} = 0.02 $ (ほどほどに振動的な減衰振動)を選んで $ K_u $ $ T_u $を決めてやると、
\[ K_u = \frac{\zeta^2}{\zeta_1^2} -1= (2500 \zeta^2 -1) \]
\[ \frac{1}{ T_u} = \frac{1}{2 \pi} (\omega_n \sqrt{1 + K_u}) \sqrt{1-\zeta_1^2} = \frac{1}{2 \pi} \omega_n  \sqrt{1 + K_u} \sqrt{1-(\frac{\zeta}{\sqrt{1 + K_p}})^2} = \frac{1}{2 \pi} \omega_n \sqrt{1-\zeta^2+K_u} \]
となり、
\[ K_p = 0.45 (2500 \zeta^2 - 1) \approx 1125 \zeta^2 \]
\[ K_i  = \frac{1}{2 \pi 0.83} \omega_n K_p \sqrt{1-\zeta^2+K_u}  \approx 10786 \omega_n \zeta^3 \]
と求まります.
(結構大きいゲインですね…)

2.1 お試しのシミュレーション

チューニングの成果を確認がてらシミュレーションしてみることにした.
シミュレーションにはScilab 5.4.1を使用. コードはページ末に付けました.

制御対象 : $ \omega_n = 1 $ [rad/s],  $ \zeta = 1 $ の2次遅れ系.
\[ P(s) = \frac{1}{s^2+2s+1} \]

まず、$ K_u $ を前節に従って求めると $ K_u = 2500 $.

ステップ応答の減衰がいやに早く感じますが結果は以下(コードは[6]).
理論とほぼ等しく約8Hzの振動.

これから、$ Kp = 0.45 Ku = 1125 $, $ Ki = \frac{Kp}{0.83Tu} =  10786 $と求まり、これを適用した結果がこちら.
発振していますね…
(Scilabのコードは[7])

原因を検討すると、特性方程式$ 1+PC=0 $の解(=閉ループ伝達関数の極)が(3.5±34.3i, -9.0 +5.8i)で、
複素平面右側に極がいることが分かります. 発振してしまうわけです.

参考

[5] The 2nd-order lag element
http://www.atp.ruhr-uni-bochum.de/rt1/syscontrol/node31.html

[6] Kuのシミュレーション用コード(Scilab)
s=%s;

w = 1
zeta = 1

Ku = 2500 * zeta^2
Kp = 0.45*Ku
Ki =  10786 * zeta^2 * w

t=linspace(0,3,1000);

G = (w^2*Ku)/(s^2+2*zeta*w*s+w^2*(Ku+1))

sys=syslin('c',G)    
y=csim('step',t,sys);
plot(t,y)
xtitle('Step Response','Time(sec)','')

[7] ジーグラ・ニコルス法適用のコード(Scilab)
s=%s;

w = 1
zeta = 1

Ku = 2500 * zeta^2
Kp = 0.45*Ku
Ki =  10786 * zeta^2 * w

t=linspace(0,3,1000);

G = (w^2*Kp*s+w^2*Ki)/(s^3+2*zeta*w*s^2+w^2*(Kp+1)*s+w^2*Ki)

sys=syslin('c',G)    
y=csim('step',t,sys);
plot(t,y)
xtitle('Step Response','Time(sec)','')

2015/02/01

PI制御チューニング : Ziegler Nicholsの限界感度法

PID制御のチューニングのテクニックとして、
Ziegler Nicholsの限界感度法(ジーグラ・ニコルス法)が有名なので、
チューニング結果がどのようなものか確認してみることにしました.

ジーグラ・ニコルス法が適用できる前提条件は、
[1]の文献では特に前提がないもの (参考例ではむだ時間のない3次遅れ系)
[2]によるとむだ時間 + 1次遅れ系 ($ \frac{K}{1+Ts}e^{-Ls} $)のときのみ
[3]ではむだ時間+1次遅れ系($ \frac{K}{1+Ts}e^{-Ls} $)ないしはむだ時間+積分系($ \frac{K}{s}e^{-Ls} $)
としている.

むだ時間をパデ近似で有利関数表現するとどこまでも次数を上げることができるので、
[1], [2], [3]を読み合わせるて解釈すると、
むだ時間がある1次遅れ系相当か、
ある程度の次数の有理関数表現ができるプラントになら適用できるということだろう.
(後者を言ってしまうとなんでも…か…?)

追記 : ジーグラ・ニコルス法は、(むだ時間のない) 積分系・1次遅れ系・2次遅れ系に適用できないとのこと.
以降の検討で「1. ジーグラ・ニコルスの限界感度法と1次遅れ系」「2. ジーグラ・ニコルスの限界感度法と2次遅れ系」は、それを裏付ける結果になっている.
むだ時間を含む1次遅れ系では有効ですが結果については「3. ジーグラ・ニコルスの限界感度法とむだ時間のある1次遅れ系」を参照ください.
Note that the Ziegler-Nichols’ closed loop method can be applied only to processes having a time delay or having dynamics of order higher than 3. Here are a few examples of process transfer function models for which the method can not be used:
  • $ H(s) = \frac{K}{s} $ (integrator) (10.5)
  • $ H(s) = \frac{K}{T s + 1} $ (first order system) (10.6)
  • $ H(s) = \frac{K}{(\frac{s}{ω_0})^2 + 2ζ\frac{s}{ω_0} + 1}$ (second order system) (10.7)
http://home.hit.no/~hansha/documents/control/theory/tuning_pid_controller.pdf


ジーグラ・ニコルスの限界感度法は、ステップ応答法やCHR法と行ったチューニング手法と違い、むだ時間がある系のチューニングを明示的には前提としていないものの実態は、むだ時間のある系を対象としていることがわかります。

むだ時間を考えないでチューニングする場合は

などを適用すると良いと思います。
-追記ここまで-


0. ジーグラ・ニコルスの限界感度法の段取り

  1. $ K_i = 0, K_d =0 $として$ K_p $ を0からだんだん大きくしていき、
    系が振動的になったところで止める
    このときの $ K_p $ を $ K_u $とし、
    振動周期を $ T_u $とする
  2. $ K_u, T_u $を使って、$ K_p, K_d, K_i $を定める
    決め方は表の通り
制御の種類$ K_p $$ K_d $$ K_i $
P制御$ 0.5 K_u $00
PI制御$ 0.45 K_u $0$ \frac{0.45 K_u}{0.83 T_u} $
PID制御$ 0.6 K_u $$ 0.075 T_u K_u $$ \frac{0.6 K_u}{0.5 T_u} $



1. ジーグラ・ニコルスの限界感度法と1次遅れ系

まず、むだ時間のない1次遅れ系にPI制御を適用することを考えてみる.

この時、
\[ P_0(s) = \frac{1}{1+Ts} \]
\[ C_0(s) = K_p + \frac{K_i}{s} \]
従って、フィードバック系の伝達関数は、
\[ G_0(s) = \frac{P_0 C_0}{1+P_0 C_0} = \frac{K_p s + K_i}{Ts^2+(1+K_p)s+K_i} \]

すると、段取りに従い、安定限界を求めるべく$ K_i = 0 $とすると、
\[ G_0(s) = \frac{K_p}{Ts + (K_p +1)} \]
とただの1次遅れ系なってしまい安定限界を求めるどころではない.

つまり、この時は、$ K_u $ を求めようにも、
 $ K_p $ を上げたいだけ上げても安定な状態になってしまうので、
ジーグラ・ニコルスの限界感度法を適用するのは少し難しい.

だが実際に1次遅れ系にPI制御を適用するとき($ K_i > 0 $)は、
\[ G_0(s) = \frac{K_p s + K_i}{Ts^2+(1+K_p)s+K_i} = \frac{(2\zeta \omega_n - 1) s + \omega_n^2}{s^2+2 \zeta \omega_n s+\omega_n^2} \]
と、閉ループ伝達関数は2次遅れ系になる.
つまり、$ \zeta = \frac{K_p + 1}{2T \sqrt{K_i}}$ を小さく取り過ぎると振動が持続するので注意.

   以上のことをひとことで言うと、1次遅れ系で、むだ時間なしの場合はジーグラ・ニコルスの限界感度法でKp, Kiは求まらない
   (但し書きとしてアクチュエータの出力制限がないとか、離散化が無視できるとか云々)

2. ジーグラ・ニコルスの限界感度法と2次遅れ系

エントリーが大きくなってきたので、こちらへ分配 ⇒

>> PI制御チューニング : Ziegler Nicholsの限界感度法 - 2次遅れ系への適用 (むだ時間なし)


3. ジーグラ・ニコルスの限界感度法とむだ時間のある1次遅れ系

ジーグラ・ニコルスがそれなりに上手く行く条件です.



参考

[1] 古典制御 - 吉川 恒夫
      

[2] Ziegler-Nichols Tuning Rules for PID
 http://www.mstarlabs.com/control/znrule.html

[3]制御工学 (JSMEテキストシリーズ)
       

[4] The Time Delay
http://lpsa.swarthmore.edu/BackGround/TimeDelay/TimeDelay.html


2状態1入力システム:最適制御入力 (2次遅れ系の状態フィードバック制御)

バネマスダンパ系に対応した状態制御を考えるために、
2状態1入力システムに対する状態フィードバック+積分制御する系を考える
(1状態1入力システムと同様)

バネマスダンパ系プラントを扱ったエントリーと同じM, c, kを使い
$ A\equiv \begin{bmatrix}  0 & 1 \\ -\frac{c}{M}&-\frac{k}{M} \end{bmatrix}\equiv \begin{bmatrix}  0 & 1 \\ -a_1&-a_2 \end{bmatrix}$ ,$ B\equiv  \begin{bmatrix} 0 \\ \frac{b}{M} \end{bmatrix}\equiv  \begin{bmatrix} 0 \\ K \end{bmatrix}$
そして、$ X \equiv \begin{bmatrix} x \\ \dot{x} \end{bmatrix} $とすると、

$ X $を行列のまま扱うと少しわかりにくいので
同伴型[1]に展開してブロック線図を描くと

この時、
\[
\frac{d}{dt} \begin{bmatrix} z \\ x \\ \dot{x} \end{bmatrix}
= \begin{bmatrix}
0 & 1 & 0\\
0 & 0 & 1\\
0 & -a_1 & -a_2
\end{bmatrix}
\begin{bmatrix} z \\ x \\ \dot{x} \end{bmatrix}
+\begin{bmatrix} 0 \\ 0 \\ K \end{bmatrix}u \]

この系に対して最適制御入力を求める。
\[J = \int_{0}^{\infty}\left \{ \begin{bmatrix} z & x & \dot{x} \end{bmatrix}\begin{bmatrix}
1 & 0 & 0 \\
0 & q_1 & 0 \\
0 & 0 & q_2 \\
\end{bmatrix} \begin{bmatrix} z \\ x \\ \dot{x} \end{bmatrix}+ru^2\right \}dt\]

この方程式を満たす対称正定行列 Pはリカッチ方程式
\[A^T P+PA+Q-PBR^{-1}B^TP=0\]
を満たす。

ただし、
\[P\equiv \begin{bmatrix} X_1 & X_2 & X_3 \\ X_2 & X_4 & X_5 \\ X_3 & X_5 & X_6 \\ \end{bmatrix}, X_1,X_2,X_3,X_4,X_5,X_6>0\]

これを解いていくと
\[\begin{bmatrix}
-\frac{K^2 X_3^2}{r}+1 & -\frac{K^2 X_3 X_5}{r}-a_1 X_3 + X_1 & -\frac{K^2 X_3 X_6}{r}-a_2 X_3 + X_2 \\
-\frac{K^2 X_3 X_5}{r}-a_1 X_3 + X_1 & -\frac{K^2 X_5^2}{r}-2 a_1 X_5 + 2 X_2 + q_1 & -\frac{K^2 X_5 X_6}{r}-a_1 X_6 - a_2 X_5 + X_3 + X_4 \\
-\frac{K^2 X_3 X_6}{r}-a_2 X_3 + X_2 &  -\frac{K^2 X_5 X_6}{r}-a_1 X_6 - a_2 X_5 + X_3 + X_4 &  -\frac{K^2 X_6^2}{r}-2 a_2 X_6 + 2 X_5 +q_2  \\
\end{bmatrix}=0 \]

これの一般解が複雑なので、式で表すと、
\[ X_3 = \frac{\sqrt{r}}{K} \]
\[ X_5 = \frac{K^2 X_6^2}{2r}+a_2X_6-\frac{q_2}{2}\]
\[ \frac{K^6}{4r^3}X_6^4 + \frac{K^4a_2}{r^2}X_6^3
+(\frac{K^2 a_2^2}{r}+\frac{a_1 K^2}{r}+\frac{q_2 K^4}{2r^2})X_6^2\\
+(\frac{K^2a_2q_2}{r}+2a_1 a_2+\frac{2}{\sqrt{r}})X_6
+(-\frac{2\sqrt{r}a_2}{K}-q_1+a_1 q_2)=0\]

$ X_6 $は4次方程式. この解を求めるには…>> 四次方程式の解 - 高精度計算サイト 

そして、$ X_6 $は以下を満たすように選ぶ.
\[ (X_2=) \frac{X_6}{\sqrt{r}}+\frac{a_2 \sqrt{r}}{K} > 0 \]
\[ (X_1=) \frac{K X_4}{\sqrt{r}}+\frac{a_1 \sqrt{r}}{K}>0 \]
\[ (X_5=) \frac{K^2 X_4 X_6}{r}+a_1 X_6 + a_2 X_4 - \frac{\sqrt{r}}{K}>0 \]

最適制御フィードバックゲインは
\[FB=-R^{-1}B^TP=-\frac{K}{r}\begin{bmatrix} X_3 & X_5 & X_6 \end{bmatrix}  \]
この行行列の1項目からPID制御でいうところの$ K_i, K_p, K_d $に相当する.
(この場合、最適制御はPID制御チューニングの方法の一種ともいえます)

[1] : 制御工学 (JSMEテキストシリーズ) など
         

2015/01/11

1状態1入力システム:最適制御入力 (1次遅れ系の状態フィードバック制御)

1状態1入力システムに対して定常偏差を補償できるよう
状態フィードバック + 積分制御する系を考えます[1]

\[\frac{d}{dt}\binom{z}{x}=\begin{bmatrix}
0 & 1\\
0 & -a
\end{bmatrix}\binom{z}{x}+\binom{0}{K}u\]

この系に対して最適制御入力を求める
(※ この時プラントは $ \tau = 1/a $, ゲイン $ 1/a $ の1次遅れ系)

\[J = \int_{0}^{\infty}\left \{ (z, x)\begin{bmatrix}
1 & 0\\
0 & q
\end{bmatrix} \binom{z}{x}+ru^2\right \}dt\]

この方程式を満たす対称正定行列 Pはリカッチ方程式
\[A^T P+PA+Q-PBR^{-1}B^TP=0\]
を満たす。[1] ($ A\equiv \begin{bmatrix} 0 & 1 \\ 0 & -a \end{bmatrix}$ ,$ B\equiv  \binom{0}{K}$)
ただし、
\[P\equiv \begin{bmatrix}
X_1 & X_2\\
X_2 & X_3
\end{bmatrix}, X_1,X_2,X_3>0\]

これを解いていくと、
\[\begin{bmatrix}
-\frac{K^2X_2^2}{r}+1 & \frac{K^2X_2X_3}{r}-aX_2+X_1\\
\frac{K^2X_2X_3}{r}-aX_2+X_1 & -\frac{K^2X_3^2}{r}+2X_2-2aX_3+q
\end{bmatrix}=0\]

\[P = \frac{1}{K}\begin{bmatrix}
\sqrt{a^2r+2r^{\frac{1}{2}}+q} & \sqrt{r} \\
\sqrt{r} & -ar+ \sqrt{a^2r^2+2r^{\frac{3}{2}}+rq}
\end{bmatrix}\]
が求まる。

したがって、最適制御フィードバックゲインは
\[FB=-R^{-1}B^TP\\
=\begin{bmatrix}
-\frac{1}{\sqrt{r}} & a-\sqrt{a^2+\frac{2}{\sqrt{r}}+\frac{q}{r}}
\end{bmatrix}\]

このときPI制御のゲイン $ K_p, K_i $と比較すると以下のように一致する.
\[K_p = - a+\sqrt{a^2+\frac{2}{\sqrt{r}}+\frac{q}{r}}\]
\[K_i = \frac{1}{\sqrt{r}}\]

サーボ系で考える場合、目標の与え方によってI+P制御となったり、PI制御となったり違いが出ます.

1. 通常の状態フィードバック(I+P制御相当):
\[u=\begin{bmatrix}
K_i & K_p
\end{bmatrix}\left(\binom{r}{0}-\binom{z}{x}\right)\]



2. PI制御相当(のはず):
\[u=\begin{bmatrix}
K_i & K_p
\end{bmatrix}\left(\binom{r}{r}-\binom{z}{x}\right)\]

参考

[1]: Integral action in state feedback control - Prof. Alberto Bemporad
        http://cse.lab.imtlucca.it/~bemporad/teaching/ac/pdf/08-integral-action.pdf

[2] : 制御工学 (JSMEテキストシリーズ) など
         

2015/01/04

連続系PID制御 : DCモータの速度制御

プラントをDCモータに限定して、電圧(PWMなど)で速度制御する場合を考えます.
今回は、PID制御とI+PD制御について伝達関数を求めます.

1. PID制御の場合

このシステムのブロック線図は、

記号LRKeJ
定義・単位電機子インダクタンス
[H]
電機子抵抗
[Ohm]
誘起電圧定数
[V/(rad/s)]

トルク定数
[Nm/A]

※ 両者等しい 
モータ軸イナーシャ
[kg m^2]

これを等価変換すると、バネマスダンパ系にPID制御を適用した場合と同等のシステムと考えられます.

この時の負荷トルクTd[Nm]および速度指令値ω*[rad/s]に対する応答は、

\[ \omega = \frac{K_e K_d s^2+K_e K_p s+K_e K_i}{LJs^3+(RJ+K_e K_d)s^2+(K_e^2+K_e K_p)s+K_e K_i}\omega^{*} - \frac{s}{LJs^3+(RJ+K_e K_d)s^2+(K_e^2+K_e K_p)s+K_e K_i} (Ls+R)Td \]


2. I+PD制御の場合

ブロック線図は、

これを同様に等価変換すると

この時の負荷トルクTd[Nm]および速度指令値ω*[rad/s]に対する応答は、

\[ \omega = \frac{K_e K_i}{LJs^3+(RJ+K_e K_d)s^2+(K_e^2+K_e K_p)s+K_e K_i}\omega^{*} - \frac{s}{LJs^3+(RJ+K_e K_d)s^2+(K_e^2+K_e K_p)s+K_e K_i}(Ls+R)Td \]

となります.

同じ$ K_p, K_d, K_i $のパラメータを使った場合で比較すると、
PID制御と比較してI+PD制御は
  ・ 目標追従の面で少しおとなしい応答
  ・ 外乱に対しては同じ応答性を有する
といえそう.

その分、逆に、その分、$ K_p, K_d, K_i $をもう少し過激にチューニングして、
外乱に対する応答を良くする考え方もできるのではないだろうか?

I+PD制御という構成は状態フィードバック制御で良く使われる考えだと思います.