Butterfly_Effect( )

地面と別れる方法

線形システムの安定性を調べる

Dynamics

octaveでisstable( )する話．
追記：octave-forgeのcontrolコンポーネントの中にisstable()関数がありました．そちらを使うことを推奨します．

線形システムの内部安定性

フライトシミュレーションの記事でも扱ったけれど，線形なシステムは一般に

\displaystyle \dot{\mathbf{x}} = \mathbf{Ax+Bu}

という形で表記できる．特に \mathbf{u}=0のときには

\displaystyle \dot{\mathbf{x}} = \mathbf{Ax}

となるわけだけど，この時にシステムが安定かどうか（つまり内部安定性）は行列 \mathbf{A}によって決定付けられる．手計算でやるならラウスやフルビッツの方法を使うけれど，計算機があるならそんな手間は不要．そのまま特性方程式の根を求めてやって，根の実部が全部負だったら安定ということになる．

そういった処理を行うとき，Matlabだとisstable( )関数で一発なんだけど，octaveには残念ながらisstable( )関数が無い．poly( )関数で特性多項式を出して，roots( )関数で根を求めて，根の一覧を眺めてナルホドしてもいいんだけど，ちょっと面倒くさい．

octaveでisstable( )する

というわけでisstable( )関数を書いた．行列 \mathbf{A}の特性方程式の根の実部がすべて負であれば安定なので，それを書いてやるだけ．安定ならば1, 不安定ならば0を返す関数になっている．関数スクリプトなので，ファイル名は"isstable.m"にすること．

関数スクリプト

% isstable.m
function T_F = isstable(A)
	%特性多項式を計算
	pl = poly(A);

	%根を計算
	r = roots(pl);

	%安定判別. 根の実部に正のものがあれば不安定．
	T_F = (length(find(real(r)>0))==0);

	if (T_F==true)
		printf('the given system is stable.\n')
	else
		printf('the given system is unstable.\n')
	end
end

使い方

行列をisstable( )関数に突っ込む．

octave:1> A = [-1,2;-4,-5];
octave:2> isstable(A);
the given system is stable.

安定らしい．確認してみる．

octave:3> roots(poly(A))
ans =

  -3.0000 + 2.0000i
  -3.0000 - 2.0000i

実部は確かに負だ．

以上．

Octaveの精義―フリーの高機能数値計算ツールを使いこなす

Octaveの精義―フリーの高機能数値計算ツールを使いこなす

フィードバック制御入門 (システム制御工学シリーズ)

フィードバック制御入門 (システム制御工学シリーズ)