MATLAB,Python,Scilab,Julia比較 その59【状態空間モデル⑰】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia-backnumber/

はじめに
登場人物
【再掲】微分解決済みの状態空間モデル
Pythonコード
シミュレーション結果
まとめ

はじめに

前回は状態空間モデルの演算をベクトル行列で行ったものをMATLABで実施。

今回は、これのPython版

登場人物

博識フクロウのフクさん

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1

エンジニア歴8年の太郎くん

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1

【再掲】微分解決済みの状態空間モデル

フクさん

まずは、微分解決済みの状態空間モデルの再掲。

状態方程式
\(\boldsymbol{x}(t+\Delta t)=\boldsymbol{x}(t)+\{A\boldsymbol{x}(t)+B\boldsymbol{u}(t)\}\Delta t \)

出力方程式
\(\boldsymbol{y}(t+\Delta T)=C\boldsymbol{x}(t+\Delta t)+D\boldsymbol{u}(t)\)

太郎くん

ベクトル、行列の演算だけでいけるのは分かったから、Python(Numpy)でも似た感じでいけそうだね。

Pythonコード

フクさん

そしてコード。

import numpy as np
import matplotlib.pyplot as plt

def statespacemodel(A,B,C,D,u,dt,x):
  # 状態方程式
  x = x + (A@x + B@u) * dt
  
  # 出力方程式
  y = C@x + D@u
  return x,y


m=1

A=np.array([[0,0],[1,0]])
B=np.array([[1/m],[0]])
C=np.array([[1,0],[0,1]])
D=np.array([[0],[0]])
dt=0.001

t = np.linspace(0, 10, 10000) # 時間(横)軸
u = np.zeros((1,10000));	# 入力信号生成
u[0][5000:10000]=1	# 5秒後に0から1へ
y=np.zeros((2,len(t)))
x=np.zeros((2,1))

for i in range(0, len(t)):
  x,y[:,[i]] = statespacemodel(A,B,C,D,u[:,[i]],dt,x)

plt.plot(t.T,y.T)
plt.plot(t.T,u.T, "--b")
plt.grid();

plt.show()

太郎くん

これも演算としてはMATLABと一緒だね。

フクさん

以前もいったが、Numpyの場合、内積は「@」になるんで、そこは注意が必要だな。