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制御という構成は状態フィードバック制御で良く使われる考えだと思います.