Pythonで学ぶ制御工学 第11弾：安定性

#Pythonで学ぶ制御工学< 安定性 >

##はじめに
基本的な制御工学をPythonで実装し，復習も兼ねて制御工学への理解をより深めることが目的である．
その第11弾として「安定性」を扱う．

##安定性
まず，伝達関数についての安定性について，図を使っての説明を以下に示す．

BIBO安定の定義を学んだところで，次に安定性を判定するために極との関係性を含めた説明を以下の図に示す．

これで，BIBO安定の定義とその必要十分条件を学んだ．しかしながら，これは伝達関数から得られる安定性の判断である．
続いて示すのは，状態空間モデルにおける安定性の判断である．先ほど同様に，以下の図にて説明を示す．

これで，漸近安定を学んだ．また，漸近安定とBIBO安定との関係性についても知ることができた．

##実装
先ほどの説明に従い，いくつか簡単な例を参考にPythonで実装して，安定性についての理解を深める．プログラムはBIBO安定と漸近安定の2つある．以下にそれぞれのソースコードとそのときの出力を示す．

#####ソースコード①：BIBO安定

BIBO_stability.py

"""
2021/03/04
@Yuya Shimizu

入出力安定（BIBO安定）
"""
from control import tf, pole

##伝達関数の作成
#1次遅れ系
P1 = tf([0, 1], [1, 1])     #1/(s+1)
P2 = tf([0, 1], [-1, 1])     #1/(-s+1)

#2次遅れ系
P3 = tf([0, 1], [1, 0.05, 1])     #1/(s^2+0.05s+1)
P4 = tf([0, 1], [1, -0.05, 1])     #1/(s^2-0.05s+1)


##極の計算
#1次遅れ系
pole1 = pole(P1)
pole2 = pole(P2)

#2次遅れ系
pole3 = pole(P3)
pole4 = pole(P4)

print(f"<極によるBIBO安定判定>\n-----1次遅れ系-----\n(例1){P1}\n極:s={pole1[0]}\t\t⇒\tBIBO安定である")
print(f"\n(例2){P2}\n極:s={pole2[0]}\t\t⇒\tBIBO安定でない")
print(f"\n-----2次遅れ系-----\n(例3){P3}\n極:s={pole3}\t\t⇒\tBIBO安定である")
print(f"\n(例4){P4}\n極:s={pole4}\t\t⇒\tBIBO安定でない")

##別記法
"""
#極が分母多項式の根ということより，次のようにして記述することもできる

from control import tfdata
import numpy as np

#分子多項式と分母多項式に分ける．
Dp = []
P = [P1, P2, P3, P4]
for p in P:
    [[N]], [[D]] = tfdata(p)
    Dp.append(D)

Pole = []
for dp in Dp:
    root = np.roots(dp)
    Pole.append(root)

for i in range(len(Pole)):
    print(f"\n(例{i+1}){P[i]}\n極:s={Pole[i]}\n")
"""

#####出力

<極によるBIBO安定判定>
-----1次遅れ系-----
(例1)
  1
-----
s + 1

極:s=-1.0		⇒	BIBO安定である

(例2)
  1
------
-s + 1

極:s=1.0		⇒	BIBO安定でない

-----2次遅れ系-----
(例3)
        1
----------------
s^2 + 0.05 s + 1

極:s=[-0.025+0.99968745j -0.025-0.99968745j]		⇒	BIBO安定である

(例4)
        1
----------------
s^2 - 0.05 s + 1

極:s=[0.025+0.99968745j 0.025-0.99968745j]		⇒	BIBO安定でない

※ソースコードにおいて，コメントアウトしている箇所は出力結果に表示していない．

#####ソースコード②：漸近安定

Asymptote_stability.py

"""
2021/03/04
@Yuya Shimizu

漸近安定
"""
from for_plot import plot_set   #別ファイルに定義
from control import tf, pole
import matplotlib.pyplot as plt
import numpy as np

"""
##固有値の計算
A = np.array([[0, 1],
                     [-4, -5]])
eigenvalue = np.linalg.eigvals(A)
"""
print(f"<固有値の計算>\nA = {A}\t\tAの固有値  →  {eigenvalue}")

##位相図面

w = 1.5     #刻み幅の設定
Y, X = np.mgrid[-w:w:100j, -w:w:100j]

A = np.array([[0, 1],
                     [-4, -5]])
eig, vec = np.linalg.eig(A)     #行列Aの固有ベクトルvecと固有値eig

#dx/dt = Axを用いて，dx/dtの成分を計算
U = A[0, 0]*X + A[0, 1]*Y
V = A[1, 0]*X + A[1, 1]*Y

t = np.arange(-1.5, 1.5, 0.01)

##描画
fig, ax = plt.subplots()

#行列Aが実固有値を持つ場合のみ不変集合をプロット
if eig.imag[0] == 0 and eig.imag[1] == 0:
    ax.plot(t, (vec[1, 0]/vec[0,0])*t, ls='--')
    ax.plot(t, (vec[1, 1]/vec[0,1])*t, ls='--')

#xの位相図面をプロット
ax.streamplot(X, Y, U, V, density=0.7, color='black')
plot_set(ax, '$x_1$', '$x_2$')
ax.set_ylim(-1.5, 1.5)
ax.set_title("Phase Drawing for state 'x'")

plt.show()

#####出力

※この図ついては，後述する．

##位相図面
先ほどの実装結果において，図が現れた．これは，位相図面と呼び，状態空間モデルの状態 $x$ の軌跡をプロットした図である．以下に改めて，安定時の図と不安定時の図を示す．なお，行列*$A$*は次に示すとおりである．
　　　　　　　　　　　　　　　　安定

A=\begin{bmatrix}
0 & 1 \\
-4 & -5
\end{bmatrix}

　　　　　　　　　　　　　　　　不安定

A=\begin{bmatrix}
0 & 1 \\
-4 & 5
\end{bmatrix}

#####安定

#####不安定

直感的に分かりやすいのは，黒の矢印が示すの状態の遷移である．安定ならば，$[0, 0]$へ収束し，不安定ならば，外へ発散していくような形となっている．また，橙色と青色の破線で示されているのは，固有ベクトルであり，次の式が成り立つ．
　　　　　　　　　　　　　　　　　　　$\dot{x}=Ax=\lambda x$
この直線の上に一度乗ってしまえば，あとはその直線に沿って状態が遷移する．この直線のことを不変空間という．

##感想
安定について，BIBO安定の復習ができた．漸近安定という名前は知らなかったが，これについても復習することができた．さらに，位相図面というものは初めて知った．これは，状態の遷移も分かり，直感的に理解しやすい図であると思い，役立ちそうだと感じた．これで，システムの安定性について理解を深められたことと思う．システム設計するときには，極配置など安定性を意識することが求められることが予想されるため，このあたりは大切になってくる．

##参考文献
Pyhtonによる制御工学入門　　南 祐樹　著　　オーム社