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」が連続系と異なるところで、

センサの値をマイコンなんかに取り込むために必須の処理です．

[Wikipedia : A/D変換回路]

このデジタル変換で、「量子化」(値の離散化)と「標本化」(時間の離散化)をいっぺんにやってしまうのですが [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

デジタル制御ではセンサのAD変換値を平滑化するために移動平均を使います．この移動平均を連続系で考えてみるために考察します．

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

カットオフ周波数$ f_c $ は、おおよそ下記で表されます．

\[ f_c = \frac{0.443}{\sqrt{N^2 - 1}} \cdot f_s \]

ステップ応答の比較

まずステップ応答を比べてみます．

計算条件：カットオフ周波数44.3Hzで合わせたフィルタ

① サンプリング周波数 1kHz, 10個の移動平均
② サンプリング周波数 2kHz, 20個の移動平均
③ カットオフ周波数 44.3Hzの１次ローパスフィルタ $ \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]を参考にしました)

ボード線図の見どころは

1. 先にあげた公式はゲインが-3.0dBになる周波数を求めている
   (4つの曲線はすべて44.3Hz, -3.0dBを通る)
2. 移動平均のフィルタとしての特性は、移動平均数 ÷ サンプリング周波数で決まる．
3. 0～100Hzの範囲での移動平均は1次のローパスフィルタより強く、2次のローパスフィルタに似ている
   (先のステップ応答で見られた通り）
4. 100～500Hzまで見てみると、ところどころ減衰が悪い周波数があり、その周波数では1次LPF程度のレベルで考えたほうが良い

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

100 Hz前後までのボード線図

もう少しひいて1kHz前後までのボード線図

係数の理屈

パターン①

参考文献[1]や[2]を参照して、
サンプリング周波数 $ f_s $で無次元化した周波数 $ \omega = 2 \pi \frac{f}{f_s} $ に対して、 移動平均の伝達関数は

\[ H(\omega) = \frac{1}{N} \frac{e^{-i \omega N/2}}{e^{-i \omega/2}} \frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)} \]

これからゲインを計算すると

\[ |H(\omega)| = \frac{1}{N} \left|\frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}\right| \]

です。これを$ \omega = 0 $ 周りでテイラー展開して(1/N*Sin(x∗N/2)/Sin(x/2) expands at x=0 - Wolfram|Alpha)

\[ |H(\omega)| \approx 1 + \frac{1}{24} (1 - N^2) \omega^2 + O(\omega^4) \]

するとカットオフ周波数 $ f_c $は$ \omega_c = 2 \pi \frac{f_c}{f_s} $ とおいて

\[ \begin{array}{rrl} & \frac{1}{\sqrt{2}} & = 1 + \frac{1}{24} (1 - N^2) \omega^2 \\ \Leftrightarrow & \frac{1}{\sqrt{2}} &= 1 +\frac{1}{24} (1 - N^2) \left(2 \pi \frac{f_c}{f_s}\right)^2 \\ \Leftrightarrow & f_c^2 &= \frac{24 (\sqrt{2} - 1)}{4 \sqrt{2} \pi^2 (N^2 - 1)} f_s^2 \\ \Leftrightarrow & f_c & \approx \frac{0.4219 \dots}{ \sqrt{N^2 - 1}} f_s \end{array} \]

ここで得た0.4219..に対して、テーラー展開の展開点 $ \omega = 0 $ から離れてしまうことによる誤差を勘案して調節した秘伝の係数が 0.44294.. ということらしい。 よく見つけるわ…。

パターン②

ある定義域 $ 0 < \omega \leq \pi/2N $ で

\[ \begin{array}{rl} |H(\omega)|&= \frac{1}{N} \left|\frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}\right| \\ &\approx \frac{1}{N \omega}{\sqrt{2(1-\cos N \omega)}} \end{array} \]

近似状況の例: plot 1/N*Sin(N*x/2)/Sin(x/2) and 1/(x*N)*sqrt(2-2*cos(x*N)) at N = 4 between x = 0 and 6 - Wolfram|Alpha

したがって、

\[ \begin{array}{rrl} & \frac{1}{\sqrt{2}} &= \frac{1}{N \omega}{\sqrt{2(1-\cos N \omega)}} \\ \Rightarrow & N \omega &\approx 2.78311 \dots \\ \Leftrightarrow & f_c &\approx \frac{2.78311 \dots / 2 \pi}{N} f_s \\ &&= \frac{0.42946 \dots}{N} f_s \end{array} \]

参考文献

[1] What is the cut-off frequency of a moving average filter? - Signal Processing Stack Exchange

[2] デジタルフィルタを計算する

[3] filter design - 3dB-Cut off frequency of moving average - Signal Processing Stack Exchange

より良い近似を求めて論争が繰り広げられていますが、0.4429..意外に良いんですよね…。

プラントにむだ時間を持つ系の制御が難しいことを、
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{aligned} 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{aligned}\] もし、この伝達関数を部分分数分解できて、ヘヴィサイドの展開定理のように、
\[ 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{aligned} \delta^*(s) &= e^{Ls} \delta(s) \\ &= (K_p \cdot s + K_i) + (\tau s^2 + s) e^{Ls} \end{aligned}\] とし、$ s = \omega i $ を代入して実部 $ \delta_r $と虚部$ \delta_i $で関数を表現すると
\[\begin{aligned} \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{aligned}\] \[ \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