Python-Control 制御工学 基本機能

本記事では、制御工学の計算に使える Python-Control の基本的な機能について記載します。Python-Control は Python 環境で Matlab のように演算することができます。そのため、Matlab を使ったことがある方であれば、似ている部分があると思います。私自身は Matlab にはほとんど触れたことがないので、本記事は Matlab の経験がない方でも分かるように書いていきます。

環境

・windows 10（64bit）
・python 3.8.2（32bit版）
・エディタ Visual Studio Code 1.45.1

Python-Control 環境構築

Python-Control 動作環境を構築する際に、slycot のインストールで苦労したため、以下の記事に困った箇所と解決方法を記載しました。

ライブラリ情報

Python-Control ライブラリ内の関数については下記サイトに情報がまとめられています。また、Python-Control は Matlab 風に作られているので、参考として Matlab の資料をみても分かりやすいと思います（ただし全く同じではないので注意してください）。

Python-Control Functions — Python Control Systems Library 0.6d documentation

MATLAB Compatibility Module — Python Control Systems Library 0.6d documentation

伝達関数モデル作成

まず手始めに、伝達関数モデルを作成してみます。ライブラリ「control」から「matlab」をインポートします。関数 tf([num], [den]) を使うことで伝達関数（Transfer Function）モデルを作成できます。numには伝達関数の分子、denには分母を入れます。

# Python上のコード
from control import matlab

Ga = matlab.tf([1], [1, 2])
Gb = matlab.tf([1], [1, 2, 3])
Gc = matlab.tf([1,2], [1, -2, 3, -4])

print(Ga)
print(Gb)
print(Gc)

結果

  1
-----
s + 2


     1
-------------
s^2 + 2 s + 3


       s + 2
---------------------
s^3 - 2 s^2 + 3 s - 4

また、次のようにパラメータを定義すれば2次系の伝達関数モデルとして扱いやすくなると思います。

# Python上のコード
zeta    = 2
wn      = 3

G = matlab.tf([wn*wn], [1, 2*zeta*wn, wn*wn])
print(G)

結果

      9
--------------
s^2 + 12 s + 9

ボード線図作成

ボード線図を作図するには関数 bode( ) を使用します。括弧の中には伝達関数を入れます。また、図を出力するためにライブラリ「matplotlib」から「pyplot」をインポートします。ライブラリがインストールされていなければコマンドプロンプトからpipを使ってインストールします。

# Python上のコード
from control import matlab
from matplotlib import pyplot as plt

zeta    = 2
wn      = 3

G = matlab.tf([wn*wn], [1, 2*zeta*wn, wn*wn])

matlab.bode(G)           # ボード線図のプロット
plt.show()

結果
パラメータは適当ですが、2次系の伝達関数でゲインが -40dB となっており、位相は 0° から -180° まで変化しているので正しく出力できているのではないかと思います。

ゲイン余裕・位相余裕

ゲイン余裕と位相余裕を出力してみます。
ゲイン余裕は、ボード線図において、位相が -180° のときにゲインが0dBに対して何dBあるのかを示したものです。位相余裕は、ボード線図において、ゲインが0dBのとき、位相が -180° に対して何° あるのかを示したものです。
ゲイン余裕と位相余裕を出力するには関数 margin(sys) または margin(mag, phase, w) を使います。今回は margin(sys) で使ってみます。

margin(sys) で使う場合の引数
sysには単入力単出力（SISO）システムの伝達関数を入れます。

margin(mag, phase, w) で使う場合の引数
ゲイン（magnitude）、位相（phase 単位はdeg）、周波数（frequencies 単位は rad/sec ）をボード線図の周波数応答データから入力します。

返り値 Return
演算結果として次のパラメータが返ってきます。
gm：ゲイン余裕 
pm：位相余裕 [deg] 
wg：ゲイン余裕の周波数（位相 -180° の周波数）
wp：位相余裕の周波数（ゲイン 0dB の周波数）

# Python上のコード
from control import matlab
from matplotlib import pyplot as plt

sys = matlab.tf([1], [1, 2, 1, 0])
print(sys)

gm, pm, wg, wp = matlab.margin(sys)
print(gm, pm, wg, wp)

matlab.bode(sys)           # ボード線図のプロット
plt.show()

結果
ゲイン余裕は 2.0dBです。位相余裕は 21.39° です。ボード線図をみてみると合っていそうな感じです。

       1
---------------
s^3 + 2 s^2 + s

2.0 21.386389751875043 1.0 0.6823278038280193

注意事項
先ほど、ボード線図作成で使った2次系の伝達関数モデルでは、ゲイン余裕と位相余裕が計算されず、「inf inf nan nan」と出力されました。推測ですが、ゲインが 0dBを跨ぎ、位相が -180° を跨いでいないと上手く計算してくれないのではないかと思います。

ナイキスト線図作成

ナイキスト線図を作図してみます。ナイキスト線図は、周波数応答 G( jω) の実部を横軸に、虚部を縦軸にとる極座標系において、角周波数 ω を 0 から ∞ まで変化させた軌跡です。関数 nyquist(sys) を使います。引数sysには伝達関数を入れます。

# Python上のコード
from control import matlab
from matplotlib import pyplot as plt

sys = matlab.tf([2, 5, 1], [1, 2, 3])
print(sys)

matlab.nyquist(sys)
plt.show()

結果

2 s^2 + 5 s + 1
---------------
s^2 + 2 s + 3

ステップ応答波形

ステップ応答波形を出力してみます。伝達関数モデルは2次系とします。ステップ応答の確認には関数 step(sys, t) を使います。引数sysには伝達関数を入れます。引数tには時間を入れます。今回、時間は数値演算ライブラリ numpyを使っています。

# Python上のコード

# 数値演算ライブラリ
import numpy as np
# 制御工学ライブラリ
from control import matlab
# グラフ描画
from matplotlib import pyplot as plt

# 伝達関数モデル作成
zeta    = 0.6
wn      = 3.0

G = matlab.tf([wn*wn], [1.0, 2.0*zeta*wn, wn*wn])
print(G)

# ステップ応答
t = np.linspace(0, 5, 1000)
(yout, T) = matlab.step(G, t)
plt.plot(T, yout)
plt.grid()
plt.axhline(1, color="b", linestyle="--")
plt.xlim(0, 5)
plt.show()

結果
パラメータは適当ですが、2次系の伝達関数で ζ＝0.6 としたのでちょとオーバーシュートした波形になります（知識は浅いですが、一応そのくらいは知っています笑）。実際の出力結果もそのような波形になっています。

       9
---------------
s^2 + 4.8 s + 9

ブロック線図（直列）

ブロック線図の直列を演算してみます。直列に接続されるブロック線図の伝達関数を導出するには関数 series( ) を使用します。

# Pytho上のコード
# 制御工学ライブラリ
from control import matlab

# ブロック線図 直列
Ga = matlab.tf([1], [1, 1])
Gb = matlab.tf([2, 2], [2])

print(Ga)
print(Gb)
print(Ga * Gb)

Gab = matlab.series(Ga, Gb)

print(Gab)

結果
ポイントとして、関数 series( ) を使わなくても「Ga*Gb」でも計算できていることが分かります。

  1
-----
s + 1


2 s + 2
-------
  2


2 s + 2
-------
2 s + 2


2 s + 2
-------
2 s + 2

ブロック線図（並列）

ブロック線図の並列を演算してみます。並列に接続されるブロック線図の伝達関数を導出するには関数 parallel( ) を使用します。

# Pytho上のコード
# 制御工学ライブラリ
from control import matlab

# ブロック線図 並列
Ga = matlab.tf([1], [1, 1])
Gb = matlab.tf([2, 2], [2])

print(Ga)
print(Gb)
print(Ga + Gb)

Gab = matlab.parallel(Ga, Gb)

print(Gab)

結果
ポイントとして、関数 parallel( ) を使わなくても「Ga+Gb」でも計算できていることが分かります。

  1
-----
s + 1


2 s + 2
-------
  2


2 s^2 + 4 s + 4
---------------
   2 s + 2


2 s^2 + 4 s + 4
---------------
   2 s + 2

ブロック線図（フィードバック）

ブロック線図のフィードバックを演算してみます。フィードバック接続されるブロック線図の閉ループ伝達関数を導出するには関数 feedback( ) を使用します。

# Pytho上のコード
# 制御工学ライブラリ
from control import matlab

# ブロック線図 フィードバック
Ga = matlab.tf([1], [1, 1])
Gb = matlab.tf([2, 2], [2])
print(Ga)
print(Gb)

Gab = matlab.feedback(Ga, Gb)
print(Gab)

結果

  1
-----
s + 1


2 s + 2
-------
  2


  2
-------
4 s + 4

注意事項
注意事項として、Python-Controlのfeedback関数では、ネガティブフィードバックになっています。

まとめ

ここまでに記載した機能を使えば、伝達関数モデルを複数作成し、システムの閉ループ伝達関数を導出して、ステップ応答を導出することができると思います。また、ボード線図を作図して、ゲイン余裕や位相余裕を見ればシステムの安定性も検討できると思います。さらに、PI補償器やPID補償器を入れればシステムの制御設計も活用できるのではないかと思います。

Python Control 関連記事

本記事では Python Control の基本的な扱い方について触れました。ここから発展させた内容の記事も書いていますので、ぜひ参考にしてみてください。

下記の記事では、電気回路を例にシステムのモデリング、制御系の設計について書いています。さらに、Python Control に加えて、回路シミュレータのPSIM を用いたシミュレーションまでを一連の流れで書いているので、ぜひ読んでみてください。

あとがき
（Visual Studio Codeの良い所）

私は初めてPythonを触ったときは、Spiderというエディタを使っていました。たしか、Anacondaの環境をインストールするとセットで付いてきたと思います。Spiderで使いにくかったのはコード補完などがない点です。
そんな時、別のプログラミング学習で Visual Studio Code を使う機会がありました。コード補完されるし、様々な言語をひとつのエディタで扱えるのでとても使いやすいと感じました。そこでPythonでも使ってみることにしたら、とても便利でした！特にPython-Controlの関数を扱う際、関数にカーソルを当てれば引数や返り値の情報を表示してくれるのがかなり良いです（下図）。初心者には使いやすいと感じました。ちなみに Visual Studio Code は無料で使えます。

以上