Butterfly_Effect( )

僕はとにかく地面と別れる方法を探しているんだ

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

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の精義―フリーの高機能数値計算ツールを使いこなす

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

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