力学系へのPI制御の利用 - how to code something

こんな状態だとする

x' 方向の運動方程式
$\ddot{x'}(t)=-kx'+mg sin\theta(t) = -kx'(t) +mg \phi(t)$

Laplace変換(初期値 x'(t)=0 phi(t)=0)
$s^2 X(s) = - \frac{k}{m}X(s) +g \Phi (s)$

伝達関数
$G(s)=\frac {X(s)}{\Phi (s)}= \frac{g}{s^2 + k/m}$

ステップ(0.5)入力（Φ=0.5→θを60度固定）

t=0:0.01:10

u=0 
for i=1:length(t) 
step(i)=0.5 
end 

clf
G=9.8/(%s^2+1)
plot2d(t,csim(step',t,syslin("c",G)),1)

周期1のサイン波を10秒間加える

t=0:0.01:10

for i=1:10*100+1 #iは１から始める（Scilabの配列のインデックスが１からなので）
u(i) = abs(sin(2*%pi*(i-1)/100))
end

plot2d(t,u)

こんな波形

このときの入力時間応答

plot2d(t,csim(u',t,syslin("c",G)))

振動周波数は変わらないが、加振によって振幅が増えている

比例制御を入れた時の安定性をみるために、根軌跡(root locus)を調べる
ゲインKpを上げた時の極の振る舞いが分かる（右半面に極が移動すると不安定）
evansコマンドが使える

evans(G)

不安定化はしないがY軸から離れるので振動成分が増すことがわかる
これに対するPI制御系を設定
$C(s)=Kp(1+\frac{1}{T_i {\rm s}})$

C=1+1/s

制御器を入れた時のナイキスト線図（一巡伝達関数のベクトル軌跡）を見る

nyquist(syslin("c",C*G))

よくわかんないので(-1+j0)付近を拡大

一応-1+j0を左に見てるから安定といえるのか
ナイキスト線図自体があまり見る形でないのでなんとも言えない　
制御系を入れた後の根軌跡

ゲインをあげると不安定になってしまうことが分かる
ステップ応答を見る

plot2d(t,csim(step',t,syslin("c",C*G)),1)

発散しておりフィードバック制御系が成り立っていない
サイン波で加振した時の時間応答

plot2d(t,csim(u',t,syslin("c",G)))

同じく発散

次に変位xではなくv（速度）を出力変数にとって考えてみる
伝達関数
$G(s)=\frac {V(s)}{\Phi (s)}= \frac{g s^2}{s^2 + \frac{k}{m}}$

G=s^2/(s^2+1)

ステップ応答とサイン波で加振した時の時間応答

こっちの方がわかりやすくはある
PI制御則を入れた時の根軌跡

evans(C*G)

比例ゲインKpをあげると収束そうな気配がある
一巡伝達関数の極

->roots(denom(C*G)+numer(C*G)) 
ans  =
 
  - 0.25 + 0.6614378i  
  - 0.25 - 0.6614378i

右半面には極は存在しない=不安定極は無い
一巡伝達関数のベクトル軌跡＝ナイキスト線図

nyquist(syslin("c",C*G))

ナイキスト線図Γは-1+j0を左に見ているので安定と言える
(これも軌跡が円弧状じゃないので良くわからない)
速度を0に収束すべくゲインを上げてみる

clf  
plot2d(t,csim(step',t,syslin("c",C*G)),1)
plot2d(t,csim(step',t,syslin("c",10*C*G)),2)
plot2d(t,csim(step',t,syslin("c",100*C*G)),3)
legend(["Kp=1" "Kp=10" "Kp=100"])

根軌跡からはゲイン→大で振動→収束って思ったけどそんなに甘くはなかった。
次に積分時間を長くして振動成分をカットすることを考える

clf 
plot2d(t,csim(step',t,syslin("c",G*(1+1/s))),1)
plot2d(t,csim(step',t,syslin("c",G*(1+1/s*0.1))),2)
plot2d(t,csim(step',t,syslin("c",G*(1+1/s*0.01))),3)
legend(["Ti=1" "Ti=10" "Ti=100"])

TI>10ではさほど違いは見られない
サイン波入力による加振に適用してみる

clf 
plot2d(t,csim(u',t,syslin("c",G*(1+1/s))),1)
plot2d(t,csim(u',t,syslin("c",G*(1+1/s*0.1))),2)
plot2d(t,csim(u',t,syslin("c",G*(1+1/s*0.01))),3) 
legend(["Ti=1" "Ti=10" "Ti=100"])

これもTi>10での改善は乏しい

最後に閉ループ系の伝達関数のボード線図

bode(syslin("c",C*G/(1+C*G)),0,10000,0.01)

ハイパスフィルタ的な特性を持っている