Python에서 Bode 다이어그램과 시스템의 시간 응답을 찾습니다.

제어 공학의 문제로 다음과 같이 출제되었기 때문에, sympy에 의해 수치 계산시켰기 때문에 비망록으로서 기록합니다. (과제 내용과 수치를 변경하고 있습니다)

질문 1)
다음 전달 함수의 Bode 다이어그램을 곱합니다.
$$G(s) =\frac{10}{s(s+4)(3s+1)} $$

from sympy import *​

def main():​
    x = Symbol("x")         # symbolとして使う変数の宣言
    y = Symbol("y")
    f = x* (x+4) *(3*x+1)         # 関数f(x)の定義
    f1 = expand(f)          # 関数f(x)を展開
    f2 = factor(f1)         # 関数f(x)を因数分解
    print("f = "+str(f))    # 計算結果の表示
    print("f1 = "+str(f1))
    print("f2 = "+str(f2))
​
if __name__ == '__main__':
    main()

f = x*(x + 4)(3*x + 1)
f1 = 3*x3+13*x2+4*x
f2 = x(x + 4)*(3*x + 1)
가 출력됩니다. (단지의 전개이므로 머리로 계산해도 문제 없습니다.) 이것보다 이하의 파라미터에 계수를 넣어 갑니다.

from control.matlab import *
from matplotlib import pyplot as plt

# 伝達関数のパラメータ
num = [ 10]     # 分子の係数
den = [3, 13, 4,0]     # 分母の係数
sys = tf(num, den)  # 伝達関数モデルの作成
bode(sys)           # ボード線図のプロット
plt.show()

보드 다이어그램이 출력되었습니다.

질문 2)
다음 시스템을 제어 할 때 시스템의 시간 응답을 찾습니다.

\frac{dx(t)}{dt} = \begin{pmatrix} 
0 & 1 \\
-4 & -5\\
\end{pmatrix}x(t) + \begin{pmatrix} 0 \\ 1 \\ 
\end{pmatrix}u(t) \\
x(0) = \begin{pmatrix} 1 \\ 0 \\ \end{pmatrix}, u(t) = Heaviside(t)

from control.matlab import *
from matplotlib import pyplot as plt
import sympy as sym
import numpy as np

def main():
    var("a:z")
    var("A:Z")
 var("X0")
    A = [[0,1], [-4,-5]]
    B = [[0], [1]]
    X0 = [[1],[0]]
    A = sym.Matrix(A)  
    X = (sym.eye(2)*s - A).inv()    
    Y = inverse_laplace_transform(X,s,t) #exp(At)  
    X0 = sym.Matrix(X0)#これはx0
   #print(Y *X0 )#第一項

    #積分の項を計算していく    
    var("tau")#積分のためのtau
    Z = inverse_laplace_transform(X,s,t-tau)#再度計算
    B = sym.Matrix(B)#b
    W = integrate(Z * B ,(tau,0,t))#積分　eAt * b * u(t)
    #print(W)
    Answer = Y * X0 + W

    print(sym.simplify(Answer[0,0]))
    print(sym.simplify(Answer[1,0]))

if __name__ == "__main__":
    main()

Heaviside(t)/4 + exp(-t)*Heaviside(t) - exp(-4*t)*Heaviside(t)/4
-exp(-t)*Heaviside(t) + exp(-4*t)*Heaviside(t)
라고 출력되었습니다. 이것으로부터, $$ x(t)= [x_1(t), x_2(t)]$$로 한 경우,
$$x_1(t) =\frac{1}{4} + exp(-t) -\frac{exp(-4t)}{4}$$

$$x_2(t) = -exp(-t) + exp(-4t)$$
이것이 해답이지만,이 함수를 염두에두고 플롯합니다.

import sympy as sym
var("t")
sym.plotting.plot((1/4 + exp(-t) - exp(-4*t)/4,(t,0,10)))
sym.plotting.plot(-exp(-t) + exp(-4*t),(t,0,10))

도움이 된 기사
htps : // / l hr thm. 조호.んふぉ / 세이 죠 리론 / py 텐 콘 t 로 루 시무 치온 /
htps : // 코 m / 응 _ 그 싸 / ms / 아다 16 에 29 에 f8282498 b7
+sympy에서 수치 계산에 관한 다양한 기사