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


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)','')