Pythonで制御 その１

いままで、制御問題を解いたり、シミュレーションするときにはMatlabをもっぱら使ってきたが、最近Pythonでも同じようなことができることに気づいたので、プロモートしてみたい。 まだまだ、機能のバラエティはMatlabよりも限られている感じ受けるものの、フリーなのがなによりも良い。

python-control.readthedocs.io

こちらのコマンドでインストールしてみる。

pip install slycot   # optional
pip install control

スモークテストとして、バネマスダンパ系のシミュレーションをやってみたい。下の図のマスｍに右方向の力を加えると、ビヨーンと振動を繰り返しながらも、ある一定の横位置に落ち着くはずだ。

リマインダとして、

---

状態空間で表現したときの、このシステムのダイナミクス（入力である力を与えたときに、変位と速度がどのように変わるか表すモデル）は、

状態Xのうち、位置のほうを観測することにすると、出力方程式は、

なお、状態Xは変位と速度から成り立っている。

---

上のPython Controlのホームページからデモ用のコードを借りてみる。

import os
import matplotlib.pyplot as plt   # MATLAB plotting functions
from control.matlab import *  # MATLAB-like functions

# Parameters defining the system
m = 250.0           # system mass
k = 40.0            # spring constant
b = 60.0            # damping constant

# System matrices
A = [[0, 1.], [-k/m, -b/m]]
B = [[0], [1/m]]
C = [[1., 0]]
sys = ss(A, B, C, 0)

と、ここまで書いてみてMatlabとシンタックスが全く同じである。。Pythonの場合は、ラインの末尾に；が不要なくらいが違いだろうか。 力を突然右方向に加えたと想定して、位置が時間ともにどのように変化するか調べてみる。必要なコードはMatlabでステップ応答を書くときとほぼ同じだ。

plt.figure(1)
yout, T = step(sys)
plt.plot(T.T, yout.T)
plt.show(block=False)

ステップ応答

一度大きくオーバーシュートしたのち、ある値に収束していくのがわかる。収束した0.025[m]というのが、このバネダンパ系に1[N]の力を加えたときに、定常的にバネが伸びる量ということになる。次に、

plt.figure(2)
mag, phase, om = bode(sys, logspace(-2, 2), Plot=True)
plt.show(block=False)

とすると、ボード線図がかける。ここもMatlabとシンタックスがほぼ同じだ。

ボード線図

ピークは0.06[Hz]くらいのところにある。図１の時刻歴の応答が大体１７秒周期、つまりで1/17(=0.059)[Hz]で波打っているので、さきほどの時刻歴と矛盾していない。

ところで、 昨今、Matlabを使わないと簡単にはできなかった様々なことが、Pythonによって可能になっていていると感じることが増えている。学生のときからMatlabをヘビーに使ってきた身ではあるが、最近Matlabが色々な分野から駆逐されていくのではないかと危惧している。実際には、時代とともにニーズを反映したより良いツールが現れて、これまでのものを淘汰していくのは望ましいことだろう。ただし、エンジニアとしては従来のスキルセットが売れなくなってくるので、新しいツールにどんどん適用していかなければならないのが辛いところである。