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