"""
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安定でない

"""
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()

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

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