Pythonプログラミング（ステップ７・配列・線形代数・制御理論）

連続時間の線形なシステム

世の中の動きを伴う現象の多くは、線形の微分方程式を使って（近似的に）表現することができる。 時刻$t$でのシステムの状態（あるいは出力）を変数 $$ \boldsymbol{x}(t) = \begin{pmatrix} x_{n-1}(t) \\ x_{n-2}(t) \\ \vdots \\ x_{0} \\ \end{pmatrix} $$ 表し、時間に依存した入力や外乱を $$ \boldsymbol{u}(t) = \begin{pmatrix} u_{n-1}(t) \\ u_{n-2}(t) \\ \vdots \\ u_0 \\ \end{pmatrix} $$ で表すことにする。 そのとき、行列$A$を使って、線形な微分方程式 $$ \frac{d}{dt} \boldsymbol{x}(t) = A \boldsymbol{x}(t) + \boldsymbol{u}(t) \tag{1} $$ によって表現できるようなシステムをここでは考える。

微分方程式の理論（あるいは制御理論）が教えるところによれば、 入力の無い場合（$\boldsymbol{u}(t)=0$)、 ある初期値から出発して、長時間的に系が安定である（全ての$x_i(t)$が有限に留まる）ためには、 行列$A$の全ての固有値の実部が負でなければならない。

高階の微分方程式 $$ \frac{d^n y(t)}{dt^n} + a_{n-1} \frac{d^{n-1} y(t)}{dt^{n-1}} + \cdots + a_1 \frac{d y(t)}{dt} = - a_0 y(t) + u(t) \tag{2} $$ で表されるようなシステムの場合は $$ x_{n-1} = \frac{d^{n-1}}{dt^{n-1}}y, \; x_{n-2} = \frac{d^{n-2}}{dt^{n-2}}y, \; \cdots, \; x_{1} = \frac{d}{dt}y, \; x_{0} = y \tag{3} $$ と置き直すと、連立の一階微分方程式 $$ \begin{eqnarray} \frac{d}{dt} x_{n-1} & = & - a_{n-1} x_{n-1} - a_{n-2} x_{n-2} - \cdots - a_0 x_{0} + u \\ \frac{d}{dt} x_{n-2} & = & x_{n-1} \\ \frac{d}{dt} x_{n-3} & = & x_{n-2} \\ \vdots & & \\ \frac{d}{dt} x_0 & = & x_1 \\ \end{eqnarray} \tag{4} $$ に書き直せるので、行列 $$ A = \begin{pmatrix} -a_{n-1} & -a_{n-2} & \cdots & -a_1 & -a_0 \\ 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & & & & \\ 0 & 0 & \cdots & 1 & 0 \\ \end{pmatrix} \tag{5} $$ とベクトル $$ \boldsymbol{u}(t)= \begin{pmatrix} u(t) \\ 0 \\ 0 \\ \vdots \\ 0 \\ \end{pmatrix} \tag{6} $$ を使って、やはり $$ \frac{d}{dt} \boldsymbol{x}(t) = A \boldsymbol{x}(t) + \boldsymbol{u}(t) $$ と表現できる。

この場合、$A$の固有値は、固有多項式 $$ \begin{vmatrix} -a_{n-1} - \lambda & -a_{n-2} & \cdots & -a_1 & -a_0 \\ 1 & - \lambda & \cdots & 0 & 0\\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & & & & \\ 0 & 0 & \cdots & 1 & - \lambda \\ \end{vmatrix} = \lambda^2 + a_{n-1} \lambda^{n-1} + \cdots + a_1 \lambda + a_0 = 0 \tag{7} $$ の解として得ることができる。

　練習： 行列式の計算

式(7)の左辺の行列式が右辺の多項式で表されることを確認してみなさい。

　ヒント

　練習： 固有値計算

NumPyの固有値計算機能を使って、式(5)の行列の固有値を計算して表示するコードを書いてみなさい。

　解説: 古典制御理論との対応

古典制御理論では、系の挙動をラプラス変換を使って議論するのが一般的である。 実際に微分方程式 (2)をラプラス変換すると $$ \left[ s^n + a_{n-1} s^{n-1} + \cdots + a_1 s \right] Y(s) = -a_0 Y(s) + U(s) + Y_0(s) \tag{9} $$ となる（$y(t)$のラプラス変換を$Y(s)$, $u(t)$のそれを$U(s)$、初期条件を$Y_0(s)$と、それぞれ置いた）。 ここで、 $$ \begin{eqnarray} G(s) & = & \frac{1}{s^n + a_{n-1} s^{n-1} + \cdots + a_1 s} \\ H(s) & = & a_0 \end{eqnarray} \tag{10} $$ と置き直すと、形式的に、入力と出力の関係を $$ Y(s) = \frac{G(s)}{1+G(s)H(s)} \left( U(s) + Y_0(s) \right) \tag{11} $$ とまとめることができる。 ここで、伝達関数（$\frac{G}{1+GH}$）の極の実部が系の安定性を決めることになるが、それには$1+GH$のゼロ点を調べればよい。 そして、$1+GH$のゼロ点を与えるのが $$ s^n + a_{n-1} s^{n-1} + \cdots + a_1 s + a_0= 0 \tag{12} $$ となる。すなわち、この解$s$は、(5)式の行列の固有値に他ならない。

システムの安定性

システムの安定性を知るには、固有方程式の解を求めれば良いが、次数が大きい場合、実際の計算は大変な手間となる。 そこで、比較的少ない工数で安定性を判別するための手続きがいくつか考案されてきた。

そのひとつに、フルビッツの安定判定法がある。

系の特性方程式が $$ s^n + a_{n-1} s^{n-1} + \cdots + a_1 s + a_0 = 0 $$ のとき、最高次の項から数えて、偶数番目の項と奇数番目の項の係数をずらしながら並べた行列 $$ H = \begin{pmatrix} a_{n-1} & a_{n-3} & a_{n-5} & a_{n-7} & \cdots & \cdots & 0 & 0 \\ 1 & a_{n-2} & a_{n-4} & a_{n-6} & \cdots & \cdots & 0 & 0 \\ 0 & a_{n-1} & a_{n-3} & a_{n-5} & a_{n-7} & \cdots & 0 & 0 \\ 0 & 1 & a_{n-2} & a_{n-4} & a_{n-6} & \cdots & 0 & 0 \\ 0 & 0 & a_{n-1} & a_{n-3} & a_{n-5} & a_{n-7} & \cdots & 0 \\ 0 & 0 & 1 & a_{n-2} & a_{n-4} & a_{n-6} & \cdots & 0 \\ \vdots & & & & & & & \\ 0 & 0 & 0 & 0 & \cdots & a_4 & a_2 & a_0 \\ \end{pmatrix} $$ を考えると、その部分行列の行列式の符号によって安定性が判定できる。

詳細は制御工学等の教科書にあたるとして、その手順をPythonでなぞってみたのが以下のコードである：

# coding: utf-8

import numpy as np

coef = np.array([1,2,2,3,1], dtype='float64')

N = len(coef)-1
H = np.zeros((N,N))

evens = coef[1::2]
odds = coef[0::2]

evens = np.pad(evens,[0,N-len(evens)],'constant',constant_values=0)
odds = np.pad(odds,[0,N-len(odds)],'constant',constant_values=0)

for i in range(0,N):
    if i%2==0:
        H[i,i//2:] = evens[:N-i//2]
    else:
        H[i,i//2:] = odds[:N-i//2]

print('H=\n',H)

for n in range(1,N+1):
    D =H[:n,:n]
    det = np.linalg.det(D)
    print('D=\n',D)
    print('det=',det)