Pythonで学ぶ制御工学 第1弾：Pythonモジュールのまとめ

#Pythonで学ぶ制御工学< Pythonモジュールのまとめ >

はじめに

基本的な制御工学をPythonで実装し，復習も兼ねて制御工学への理解をより深めることが目的である．
その第1弾としてPythonモジュールをまとめる．

Pythonモジュール

制御工学を学ぶにあたって，ここでは5つのモジュールをそれぞれソースコードとそのときの出力を示すことで簡単にまとめておく

Numpy

数値計算の基本パッケージ
効率よく，高速にさまざま数値計算や統計処理，信号処理を行うことができる．

ソースコード

learn_numpy.py

"""
2021/02/10
@Yuya Shimizu

Numpyについての簡単なまとめ
"""

import numpy as np

##基本的な数値計算

#平方根
Sqrt = np.sqrt(4)
print(f"<平方根>\n{Sqrt}\n")

#絶対値
Abs = np.abs(-5)
print(f"<絶対値>\n{Abs}\n")

#三角関数
deg = 30
rad = np.deg2rad(deg)
Sin = np.sin(rad)
Cos = np.cos(rad)
Tan = np.tan(rad)
Arcsin = np.arcsin(Sin)
Arccos = np.arccos(Cos)
print(f"<三角関数>\ndeg = {deg} → rad = {rad}\nsin = {Sin}\ncos = {Cos}\ntan = {Tan}\narcsin = {Arcsin}\narccos = {Arccos}\n")

#指数
Exp = np.exp(1)
print(f"<指数>\n{Exp}\n")

#対数
Nlog = np.log(Exp)
Log = np.log10(10)
print(f"<対数>\n{Nlog}(自然対数)\n{Log}(常用対数)\n")

#四捨五入
Round = np.round(5.55)
print(f"<四捨五入>\n{Round}\n")

#円周率
Pi = np.pi
print(f"<円周率>\n{Pi}\n")

#複素数
Con = -5+2j
Conj = np.conj(Con)
Real = np.real(Con)
Imaginary = np.imag(Con)
print(f"<複素数>\n{Con}\n共役な複素数: {Conj}\n実部: {Real}\n虚部: {Imaginary}\n")


##ベクトルや行列の演算

#定義
A = np.array([[1, 2],
                     [-3, 4]
                    ])
print(f"<行列(ベクトル)の定義>\n{A}\n")

#転置
Trans = A.T
print(f"<転置>\n{Trans}\n")

#逆行列
Inverse = np.linalg.inv(A)
print(f"<逆行列>\n{Inverse}\n")

#行列式
Determinant = np.linalg.det(A)
print(f"<行列式>\n{Determinant}\n")

#ランク(階数)
Rank = np.linalg.matrix_rank(A)
print(f"<ランク(階数)>\n{Rank}\n")

#固有値と固有ベクトル
w, v = np.linalg.eig(A)
print(f"<固有値と固有ベクトル>\neigenvalue = {w}\neigenvector = \n{v}\n")

#ノルム
x = np.array([1, 2])
Norm = np.linalg.norm(x)
print(f"<ノルム>\n{Norm}\n")

#数列の作成
Numerical_sequence = np.arange(0, 10, 1)
print(f"<数列の作成>\n{Numerical_sequence}\n")

出力

<平方根>
2.0

<絶対値>
5

<三角関数>
deg = 30 → rad = 0.5235987755982988
sin = 0.49999999999999994
cos = 0.8660254037844387
tan = 0.5773502691896257
arcsin = 0.5235987755982988
arccos = 0.5235987755982987

<指数>
2.718281828459045

<対数>
1.0(自然対数)
1.0(常用対数)

<四捨五入>
6.0

<円周率>
3.141592653589793

<複素数>
(-5+2j)
共役な複素数: (-5-2j)
実部: -5.0
虚部: 2.0

<行列(ベクトル)の定義>
[[ 1  2]
 [-3  4]]

<転置>
[[ 1 -3]
 [ 2  4]]

<逆行列>
[[ 0.4 -0.2]
 [ 0.3  0.1]]

<行列式>
10.000000000000002

<ランク(階数)>
2

<固有値と固有ベクトル>
eigenvalue = [2.5+1.93649167j 2.5-1.93649167j]
eigenvector = 
[[0.38729833-0.5j 0.38729833+0.5j]
 [0.77459667+0.j  0.77459667-0.j ]]

<ノルム>
2.23606797749979

<数列の作成>
[0 1 2 3 4 5 6 7 8 9]

Matplotlib

グラフ描画のためのモジュール
さまざまなタイプのグラフやアニメーションを作成することができる．
グラフ作成例（公式）: https://matplotlib.org/2.0.2/gallery.html

ソースコード

learn_matplotlib.py

"""
2021/02/11
@Yuya Shimizu

Matplotlibについての簡単なまとめ
"""
from matplotlib import pyplot as plt

import numpy as np

##グラフ作成
#データの準備
x = np.arange(0, 4*np.pi, 0.1)
y = np.sin(x)

#描画
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()


##オブジェクト指向のグラフ作成
"""先ほどの作成では，グラフを細部までカスタマイズして，よりきれいなグラフを作る時には不十分"""
#FigureとAxesオブジェクトの作成
fig, ax = plt.subplots()

#Axesオブジェクトの中にグラフを作成
ax.plot(x, y)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid()
fig.show()


##複数の図を比較描画
#2行1列という形で図の配置レイアウトを指定して，FigureとAxesオブジェクトを作成
fig, ax = plt.subplots(2, 1)

x = np.arange(0, 4*np.pi, 0.1)
y = np.sin(x)
z = np.cos(x)
w = y + z

#1つ目のグラフを作成
ax[0].plot(x, y, linestyle = '-', label = 'sin', color = 'k')
ax[0].plot(x, z, linestyle = '-.', label = 'cos', color = 'k')
ax[0].set_xlabel('x')
ax[0].set_ylabel('y')
ax[0].set_xlim(0, 4*np.pi)
ax[0].grid()
ax[0].legend()

#2つ目のグラフを作成
ax[1].plot(x, w, color = 'k', marker = '.')
ax[1].set_xlabel('x')
ax[1].set_ylabel('w')
ax[1].set_xlim(0, 4*np.pi)
ax[1].grid(linestyle = ':')

#グラフが重なりすぎないようにする
fig.tight_layout()
fig.show()


##図を保存
fig.savefig('fig1.pdf')

出力

はじめの2つの描画同じものを描画している．

比較描画したものが次である．

Scipy

数値計算アルゴリズムのパッケージ
信号処理や最適化，統計のための関数が用意されている．また，制御系解析・設計のための関数も含まれる．
以下に示すのは，微分方程式を解くプログラムである．扱う微分方程式は次式である．

\dot y(t) = -\frac{1}{5}y(t) + \frac{1}{5}u(t)\\

いま入力は次のとおりである．

u(t)=
\left\{
\begin{array}{ll}
\ 0& (t < 10)\\
\ 1& (t \ge 10)
\end{array}
\right.

ソースコード

learn_scipy.py

"""
2021/02/11
@Yuya Shimizu

Scipyについての簡単なまとめ
"""
#微分方程式を解くためのodeintをインポート
from scipy.integrate import odeint

import numpy as np
from matplotlib import pyplot as plt


"""dy = -1/5 y + 1/5 u
    この微分方程式の数値積分を実行する
    いま，入力を次のようにする
    u = 0 (t<10) or u = 1 (t>=10)
"""

#微分方程式の定義
def system(y, t):
    if t < 10.0:
        u = 0.0
    else:
        u = 1.0

    dydt = (-y + u)/5.0

    return dydt


#初期値と時間を設定して，微分方程式を解く
y0 = 0.5
t = np.arange(0, 40, 0.04)
y = odeint(system, y0, t)


#グラフを描画
fig, ax = plt.subplots()

ax.plot(t, y, label = 'y')
ax.plot(t, 1 * (t>=10), linestyle = '--', color = 'k', label = 'u')
ax.set_title('Differential equation: dydt = (-y + u)/5.0')
ax.set_xlabel('t')
ax.set_ylabel('y, u')
ax.legend(loc = 'best')
ax.grid(linestyle = ':')

fig.show()

fig.savefig('scipy_微分方程式.jpg')

出力

Sympy

数式処理モジュール
変数を文字のまま扱い，さまざまな計算を行うことができる．
(例)
式の展開や因数分解，微分，積分，ラプラス変換など

ソースコード

learn_sympy.py

"""
2021/02/12
@Yuya Shimizu

Sympyについての簡単なまとめ
数式処理モジュール
変数を文字のまま扱い，さまざまな計算を行うことができる
(例)
式の展開や因数分解，微分，積分，ラプラス変換まで実行可能
"""

import sympy as sp

sp.init_printing()  #これを実行しておけば，LaTeX形式の出力が表示される

##数式処理

#方程式を解く
x = sp.Symbol('x')
root = sp.solve(2*x**2 + 5*x + 3, x)
print(f"<方程式を解く>\n2*x**2 + 5*x + 3 = 0\tx = {root}\n")

#式の展開
f = sp.expand((x+1)*(x+2)**2, x)
print(f"<式の展開>\n(x+1)*(x+2)**2 = {f}\n")

#因数分解
g = sp.factor(f, x)
print(f"<式の展開>\n{f} = {g}\n")

#テイラー展開
sin_t = sp.series(sp.sin(x), x)
cos_t = sp.series(sp.cos(x), x)
print(f"<テイラー展開>\nsin(x) = {sin_t}\n")
print(f"cos(x) = {cos_t}\n")

#部分分数分解
a = sp.apart(1/((5*x)*(x+1)))
print(f"<部分分数分解>\n1/((5*x)*(x+1)) = {a}\n")

#ラプラス変換
s, t = sp.symbols('s t')
w = sp.Symbol('w', real=True)   #real=Trueは必要でないと，逆ラプラス変換のときにおかしな表示となる
laplace = sp.laplace_transform(sp.sin(w * t), t, s)
print(f"<ラプラス変換>\nsin(s) = {laplace[0]}\n")

#逆ラプラス変換
inverse_laplace = sp.inverse_laplace_transform(laplace[0], s, t)
print(f"<逆ラプラス変換>\n{laplace[0]} = {inverse_laplace}\n")

出力

<方程式を解く>
2*x**2 + 5*x + 3 = 0    x = [-3/2, -1]

<式の展開>
(x+1)*(x+2)**2 = x**3 + 5*x**2 + 8*x + 4

<式の展開>
x**3 + 5*x**2 + 8*x + 4 = (x + 1)*(x + 2)**2

<テイラー展開>
sin(x) = x - x**3/6 + x**5/120 + O(x**6)

cos(x) = 1 - x**2/2 + x**4/24 + O(x**6)

<部分分数分解>
1/((5*x)*(x+1)) = -1/(5*(x + 1)) + 1/(5*x)

<ラプラス変換>
sin(s) = w/(s**2 + w**2)

<逆ラプラス変換>
w/(s**2 + w**2) = sin(t*w)*Heaviside(t)

Control

制御工学における重要な関数などを扱える
(例)
伝達関数モデルの定義，状態空間モデルの定義　など

関数はたくさんあるため，詳細は次のサイトを参考にすればよい．Python Control System Library (https://python-control.readthedocs.io/en/0.8.4/)
※このURLは更新されて表示されない場合がある．そのときは，Python Control System Libraryで検索すれば見つかる．

ソースコード

learn_control.py

"""
2021/02/12
@Yuya Shimizu

controlについての簡単なまとめ
制御工学における重要な関数などを扱える
(例)
伝達関数モデルの定義，状態空間モデルの定義　など
"""
import control

#MATLABのような関数を使用するときには...
from control.matlab import *

特に何も示していない．先ほど記述したとおり，詳細はURL参照．なお，ソースコードにもあるように，MATLABのような関数の使用も可能のようだ．

感想

今のところ何も分からず，とりあえずモジュールを使う準備だけを行った．これから制御工学の学習で使っていく中でその仕様にも慣れていくと期待している．

参考文献

Pyhtonによる制御工学入門　　南 祐樹　著　　オーム社