はじめに
現代制御では、微分方程式を状態方程式に変換してその特性である安定性などを考察する。このとき、位置と速度を合わせたベクトルを扱うので、状態方程式は最低でも2×2の行列の演算になる。今回は、2×2の行列の演算の場合について、代数的に求めた解と数値解を比較することで状態方程式の理解を深めることを目的とする。
状態方程式の解
状態方程式を以下のように定義する。
両辺に、を掛けると以下のようになる。
指数行列の微分の逆より、
ここでt→τとしてτ=0からτ=tまでを積分すると以下のようになる。
ゆえに、これを展開すると以下のようになる。
最終的には、以下のような式になる。
ちなみに、ラプラス変換と畳込み積分を用いれば、上式を導くことができる。特殊な場合については下記の記事に記してある。
例題設定
例題として、以下のような場合で考察してみる。
としたとき、
と状態方程式は表せるものとする。
解の表示
この解について、Pythonでグラフ上にプロットするプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.linalg import expm
from scipy import integrate
from scipy.integrate import quad
#2次元のグラフ（散布図）を描くための設定
fig = plt.figure()
ax = fig.add_subplot()
A = np.array([[0,1 ], [-5,-8]])
b=np.array([1,1]).T
c=[1,0]
u=1
x_0=np.array([0,0]).T
def f(t1):
return expm(A*(t-t1))@b*u
def f1(t1):
return f(t1)[0]
def f2(t1):
return f(t1)[1]
for i in range(100):
t=i/10
seki1 = quad(f1, 0, t)[0]
seki2 = quad(f2,0,t)[0]
x_ans=x_0@expm(A*t)+np.array([seki1,seki2]).T
y=c@x_ans
ax.scatter(t, y, color='blue')
plt.xlabel("時刻t")
plt.ylabel("応答値y")
plt.show()
このプログラムにより表示されるグラフは以下のようになる。
数値解
一方で、以下のようなプログラムを書くことによって近似的にだが数値解を求めた。
細かいアルゴリズムは、以下の記事を参照してほしい。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
A = np.array([[0,1 ], [-5,-8]])
b=np.array([1,1]).T
c=[1,0]
u=1
t=0
x=0
delta_t=0.1
x=np.array([0,0]).T
#2次元のグラフ（散布図）を描くための設定
fig = plt.figure()
ax = fig.add_subplot()
while(t<10):
dxdt=A@x+b*u
x=x+dxdt*delta_t
t=t+delta_t
y=c@x
#print(y)
ax.scatter(t, y, color='blue')
plt.xlabel("時刻t")
plt.ylabel("応答値y")
plt.show()
このプログラムを実行すると以下のようなグラフが出力される。
このように、数値解と解析解はほぼ同じような応答を示すということが分かった。
まとめ
今回は、現代制御にいて、状態方程式の解の公式を導出した。また、それが正しいかどうか一次遅れ系の状態方程式を用いて、数値解と比較することによって確認した。電験では、指数行列の性質を用いて攻めるというよりは、状態方程式をラプラス変換して攻めるといった問題が多いようである。なぜなら、それによって古典制御と関連付けられるからであろう。今回のプログラムにおいて、行列を積分するプログラムを少し工夫すればよりスマートなコードになるのではないかと感じた。
コメント