manvaのエンジニアリング魂

エンジニアリング・ものづくり・DIYをもっと身近にするためのブログ。インスピレーションを刺激します。

無料で実践ロバスト制御② ーまずは古典制御ー

平田氏の著書「実践ロバスト制御」↓を無料のPython-Controlで勉強する試み②。

以下これを「教科書」と呼ぶ。

Python-Controlのインストール

Pythonを使うならAnacondaをインストールするのが簡単そう。
Pythonによる制御工学」↓

の著者の南氏のサイト↓にAnacondaやPython-Controlのインストールの仕方が書かれており,これ以上書くこともない。
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制御の話。
制御対象は剛体

\displaystyle{
y=\frac{1}{Ms^2}u
}
として,PD制御
\displaystyle{
u=(k_p+sk_d)(r-y)
}
をしたときの指令rから応答yまでの伝達関数を求めると
\displaystyle{
\frac{y}{r}=\frac{\frac{k_d}{M}s+\frac{k_p}{M}}{s^2+\frac{k_d}{M}s+\frac{k_p}{M}}
}
となる。この分母が重根\alphaとなるようにゲインk_p, k_dを決めている。
M=1とし,\alphaが-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のグラフをコピーするには,シフトを押しながら右クリック)

古典制御についてちょっと掘下げる

本筋の話とはそれるが,実応用を考える場合,性能をもっと良くするにはどうすればよいか,どれだけ良くできるか,という観点が重要になる。上図の応答は悪くはないが,良くもない。この結果を見て,古典制御だとこれくらいしかできない,と思わないように補足しておく。
まず,以下のイメージを伝えたい。

フィードバックゲインは大きいほど性能が上がるが,大きすぎると不安定になり,発散する
知っている人にとっては当たり前と思うだろうが,教科書(平田氏の著書のことではなく一般の)で理論だけ学んでいると,意外に明確には書かれていなかったりするので,そのせいで何だかピンと来ないまま勉強している初学者も多いのではないか(私がそうだった)。
なぜ明確に書かれていないかというと,理論的に明確には言えないからである。しかし,実験してみると明らかに上記傾向が経験できる。発振限界付近の挙動については,理論と実際が一致しない場合が多い。実際,今回の例のように単純化した剛体の理論では,ゲインを大きくしても不安定にはならず,ほぼパーフェクトな制御が可能である。ちょっといじってパーフェクトな制御にしてみよう。
上記で極\alphaの絶対値を大きくし,-100と-10000にしてみるとステップ応答は以下のようになる。

\alphaが -100のとき(オレンジ実線)はかなりパーフェクトに近く,-10000(青破線)だと発振してイマイチになっている。周波数特性↓を見ると,

カットオフ周波数あたりでピークができており,ゲインを大きくするほどピークが大きくなっている。「フィードバックゲインが大きいほど性能が上がる」と書いたが,ここでは悪化している。これはなぜかというと,微分sというのが厳密には実現できない(「プロパーではない」という)ので,近似微分
\displaystyle{
\frac{s}{0.01s+1}
}
にしているためだ。0.01のところを0.0001にすると,


この様になり,だいぶパーフェクトに近づいた。実際には微分は基本的に前回値と今回値の差分から求める(オイラー微分)ため,この0.01や0.0001という数字は,制御周期によって決まる。制御周期を短くするほど性能が上がりそうであることがわかる。

さて、制御周期が十分に速く,微分が近似ではなくsとみなせるとしても,零点があるので周波数特性は0dBより少し持ち上がる。そこで,もっとわかりやすくするために,上記PD制御に以下のように速度FFを加えてみる。

\displaystyle{
u=(k_p+sk_d)(r-y)-s k_dr
}
こうすると指令から応答までの伝達関数
\displaystyle{
\frac{y}{r}=\frac{\frac{k_p}{M}}{s^2+\frac{k_d}{M}s+\frac{k_p}{M}}
}
と,分子が定数項だけになり,分母を重根などとおけば単純な2次のローパスフィルタの特性になる。重根の場合,極\alphaの絶対値がカットオフ周波数[rad/s]になり,この極(すなわちフィードバックゲインk_p, k_d)は理論上いくらでも大きくできるのでパーフェクトな制御ができてしまう。実際には,上記のように微分が厳密には実現できないことや,無駄時間,無視した高周波の共振の影響(スピルオーバ)によって,ゲインを上げていくとどこかで発散する。どのくらいのゲインで発散するのかを,理論的に予測することはなかなか難しく,理論と実際が合わないことが多い。この辺がロバスト制御が生まれた背景とも言えるだろう。

MATLABPython-Controlの比較。

上記のステップ応答をMATLABでプロットしてみると↓のようになる。

横軸は拡大している。Python-Controlではパーフェクトに見えたが,MATLABではオーバーシュートが出ている。MATLABは計算誤差が大きくなりそうなところの刻み幅を自動的に細かくしてくれたりするので,たぶんMATLABの方が正確なんだろう。やはりMATLABは高くても選ばれるだけの価値があるということだ。