無料で実践ロバスト制御② ーまずは古典制御ー
平田氏の著書「実践ロバスト制御」↓を無料のPython-Controlで勉強する試み②。
以下これを「教科書」と呼ぶ。Python-Controlのインストール
Pythonを使うならAnacondaをインストールするのが簡単そう。
「Pythonによる制御工学」↓
http://y373.sakura.ne.jp/minami/wp-content/uploads/2019/06/pyctrl_anaconda20190603.pdf
私はWindows10のPCを使っているが,問題なくインストールできた。
Jupyter Labの動作確認
上記でJupyter Labがおすすめのようだったので,Jupyter Labを使うことにする。
私の環境ではAnaconda Navigatorから起動しようとしたらエラーが出たが,Anaconda Navigatorを「管理者として実行」したら起動できた。
とりあえずPython-Control本家のサイトのサンプル↓をJupyter Labにコピーして実行してみた。問題なし。
Examples — Python Control Systems Library dev documentation
アップデート
Anaconda NavigatorやJupyter Labのアップデートができない(途中で何も応答なく止まってしまう)問題があったのだが,↓の「インストール済みパッケージの一括更新など」の通り下記を実行したら更新できた。
https://www.kkaneko.jp/tools/win/anaconda3.html
conda config --remove channels conda-forge conda upgrade -y --all conda clean -y --packages
実践ロバスト制御開始(まだロバストではない)
準備ができたので,教科書の内容に移る。
まずは1.3節の例で,剛体(図1.4(a))のPD制御の話。
制御対象は剛体
とし,が-0.5のときと-2.5のときのステップ応答と相補感度関数を求めている(図1.5(a)と(b))。
教科書の「実行1.1」を移植
上記を求めるためのMATLABのサンプルコードが掲載されている(実行1.1)。これをPython-Controlでやってみた。
変更箇所
・最初に,Pythonで使うライブラリをimport
import numpy as np import matplotlib.pyplot as plt from control.matlab import *
・^ → **
Pythonではべき乗は **
・figure( ) → plt.figure( )
図の描画は,matplotlib.pyplotを使う。
・step(Tn1,'--',Tn2,10) →
(yout, T ) = step( Tn1, np.linspace(0, 10, 100) )
(yout2, T2) = step( Tn2, np.linspace(0, 10, 100) )
plt.plot( T, yout, '--', T2, yout2 )
ステップ応答を求める関数はcontrol.matlabのstep。
・上記のようにplotで2つ同時に書いたらシンプルだが,凡例をつけるには,下記のように2つに分けてlabelを設定しておく必要があるようだ。
→plt.plot( T, yout, '--', label='$\u03b1 = -0.5$')
plt.plot( T2, yout2, label='$\u03b1 = -2.5$' )
plt.legend()
ギリシャ文字のコードは↓
【保存版】ギリシャ文字 Latex, Pythonコマンドリスト | OPTY LIFE(オプティライフ)
・ylim([0 1.4]) → plt.ylim( 0, 1.4 )
・ボード線図は,Python-Controlにはゲインだけプロットする関数がないので,最初,
bodemag(Tn1,'--',Tn2) → bode(Tn1,Tn2)
としてみたが,これだと位相も表示されてしまうし,細かい設定ができなさそう。bodeではplot=Falseとして,グラフは書かず,いったんmagなどの変数に保存し,semilogxで描けば,ゲインだけプロットしたり破線などの設定ができた。
→mag,phase,om = bode(Tn1, logspace(-1,4), plot=False);
mag2,phase2,om2 = bode(Tn2, logspace(-1,4), plot=False);
plt.semilogx(om, 20*np.log10(mag), '--');
plt.semilogx(om2, 20*np.log10(mag2));
移植したコード「実行1.1」
import numpy as np import matplotlib.pyplot as plt from control.matlab import * s=tf('s') M = 1 Pn = 1/(M*s**2) alpha = -0.5; K1 = alpha**2*M - 2*alpha*M*s/(0.01*s+1); alpha = -2.5; K2 = alpha**2*M - 2*alpha*M*s/(0.01*s+1); Tn1 = feedback(Pn*K1,1); Tn2 = feedback(Pn*K2,1); plt.figure(1) (yout, T ) = step( Tn1, np.linspace(0, 10, 100) ) (yout2, T2) = step( Tn2, np.linspace(0, 10, 100) ) plt.plot( T, yout, '--', label='$\u03b1 = -0.5$') plt.plot( T2, yout2, label='$\u03b1 = -2.5$' ) plt.ylim( 0, 1.4 ) plt.legend() plt.figure(2) mag, phase, om = bode(Tn1, logspace(-1,4), plot=False); mag2,phase2,om2 = bode(Tn2, logspace(-1,4), plot=False); plt.semilogx( om, 20*np.log10( mag), '--'); plt.semilogx(om2, 20*np.log10(mag2));
実行結果
教科書の図1.5(a)(b)と(だいたい)同じ図を表示できた。
(ちなみに,JupyterLabのグラフをコピーするには,シフトを押しながら右クリック)
古典制御についてちょっと掘下げる
本筋の話とはそれるが,実応用を考える場合,性能をもっと良くするにはどうすればよいか,どれだけ良くできるか,という観点が重要になる。上図の応答は悪くはないが,良くもない。この結果を見て,古典制御だとこれくらいしかできない,と思わないように補足しておく。
まず,以下のイメージを伝えたい。
なぜ明確に書かれていないかというと,理論的に明確には言えないからである。しかし,実験してみると明らかに上記傾向が経験できる。発振限界付近の挙動については,理論と実際が一致しない場合が多い。実際,今回の例のように単純化した剛体の理論では,ゲインを大きくしても不安定にはならず,ほぼパーフェクトな制御が可能である。ちょっといじってパーフェクトな制御にしてみよう。
上記で極の絶対値を大きくし,-100と-10000にしてみるとステップ応答は以下のようになる。
が -100のとき(オレンジ実線)はかなりパーフェクトに近く,-10000(青破線)だと発振してイマイチになっている。周波数特性↓を見ると,
カットオフ周波数あたりでピークができており,ゲインを大きくするほどピークが大きくなっている。「フィードバックゲインが大きいほど性能が上がる」と書いたが,ここでは悪化している。これはなぜかというと,微分というのが厳密には実現できない(「プロパーではない」という)ので,近似微分にしているためだ。0.01のところを0.0001にすると,
この様になり,だいぶパーフェクトに近づいた。実際には微分は基本的に前回値と今回値の差分から求める(オイラー微分)ため,この0.01や0.0001という数字は,制御周期によって決まる。制御周期を短くするほど性能が上がりそうであることがわかる。
さて、制御周期が十分に速く,微分が近似ではなくとみなせるとしても,零点があるので周波数特性は0dBより少し持ち上がる。そこで,もっとわかりやすくするために,上記PD制御に以下のように速度FFを加えてみる。
こうすると指令から応答までの伝達関数はと,分子が定数項だけになり,分母を重根などとおけば単純な2次のローパスフィルタの特性になる。重根の場合,極の絶対値がカットオフ周波数[rad/s]になり,この極(すなわちフィードバックゲイン)は理論上いくらでも大きくできるのでパーフェクトな制御ができてしまう。実際には,上記のように微分が厳密には実現できないことや,無駄時間,無視した高周波の共振の影響(スピルオーバ)によって,ゲインを上げていくとどこかで発散する。どのくらいのゲインで発散するのかを,理論的に予測することはなかなか難しく,理論と実際が合わないことが多い。この辺がロバスト制御が生まれた背景とも言えるだろう。