伝達関数と状態方程式 - how to code something

１入力１出力線形時不変系の状態方程式

$\dot{x} = \bold{A}x + bu$
$y = cx$

の伝達関数Gを求める。
x(t)の初期条件をx(0)=0として上式をLaplace変換する。

$s \bold{x}(s) = \bold{A}\bold{X}(s) + b U(s)$
$Y(s) = c \bold{X}(s)$

これより、
$\bold{X}(s) = (s \bold{I} - \bold{A} )^{-1} b U(s)$
したがって
$Y(s) = c(s I - A)^{-1} b U(s)$

$G(s) = Y(s) / U(s)$ より、
$G(s) = (s I -A)^{-1} b = \frac{c \bullet adj(s I -A) b}{det(s I - A )}$

これが、伝達関数と状態方程式を結びつける関数（制御のしょっぱなに習った）
これは、「状態方程式から伝達関数を導く」問題。

逆に、伝達関数が次のように与えられているとき
$G(s) = \frac{c_{n-1}s^{n-1} + c_{n-2}s^{n-2} + \ldots + c_0}{s^n + a_{n-1}s^{n-1} + a_{n-2}s^{n-2} + \ldots +a_0}$

$(A,b,c)$ の伝達関数が上記のGに等しくなるように $(A,b,c)$ を定めることを実現という。

以下、Scilabによる計算。
状態方程式から伝達関数を得る。

A=[0 1;-2 -3];b = [0;1];c=[1 2];d=0;
//syslinは線形システム
ss_sys = syslin('c',A,b,c,d);
//ss2tfで状態方程式を伝達関数に変換する
tf_sys = ss2tf(ss_sys)

逆に、伝達関数から状態方程式を得る。

s=%s
//numは分子多項式、denは分母多項式
num = 2*s +1;den=s^2+3*s+2;
//syslinは線形システムの定義。cは連続時関係を表す（離散ならd）
tf_sys = syslin('c',num,den)
//tf2ssで状態方程式に変換
ss_sys = tf2ss(tf_sys)