一般化ヘッセンベルグ形式を用いた
ディスクリプタシステムの伝達関数の数値計算法
○江角貴宏
古賀雅伸 （九州工業大学）
A Method for Numerical Computation of Transfer Function of Descriptor System
by Using Generalized Hessenberg Form
∗T. Esumi and M. Koga (Kyushu Institute of Technology)
Abstract– The purpose of this paper is to propose a new numerical method which computes the function
of descriptor system without solving the generalized eigenvalue problem and gives very accurate solutions.
Due to computation of the generalized eigenvalue and eigenvector, the existing methods for computing
the transfer function of descriptor system suffer from numerical errors. The proposed method transforms
coefficient matrices E and A into block diagonal matrices which are similar to coefficient matrices of descriptor
system in weierstrass standard form. The transfer function is directly computed from the obtained coefficient
matrices. Some illustrative examples show the validity of the proposed method.
Key Words: Descriptor System, Transfer Function, Hessenberg
1
はじめに
ディスクリプタシステム表現は、微分要素を扱うこ
とができるので非プロパーな線形システムを表現する
ことができ、さらに代数拘束条件を同時に含めること
ができるので状態空間表現より幅広いシステムを扱う
ことができる 1) ．また、非プロパーなシステムや代数
ループを含む結合システムをディスクリプタシステム
として表現する手法が提案されている 2) ．
状態空間表現から伝達関数を求めるにはファディー
ブのアルゴリズムが用いられている．一方、ディスク
リプタシステム表現からの伝達関数を求める手法は一
般化固有値や一般化固有ベクトルを求める必要がある
ため、数値計算誤差が発生しやすい
本論文では、ディスクリプタシステムの一般化ヘッセ
ンベルグ形式を用いてディスクリプタシステム表現か
ら伝達関数を求めるための数値計算法を提案する．ディ
スクリプタシステムの一般化ヘッセンベルグ形式とは
係数行列 A と E がそれぞれ上ヘッセンベルグ行列、上
三角行列となるように正則変換した形式である．提案
手法では一般化固有値を必要としないので、数値計算
誤差が発生しにくく、高速である．例題を用いて提案
手法の有効性を示す．
2
問題設定
本論文ではディスクリプタシステム
E x˙
= Ax + Bu
y
= Cx + Du
(1)
の伝達関数を求める数値計算について考える．ただし、
x は n 次のディスクリプタベクトル、u は m 次の入力
ベクトル、y は p 次の出力ベクトル、E(n × n), A(n ×
n), B(n × m), C(p × n), D(p × m) は実行列である．E
は正則ではなく、rankE = r < n と仮定する．また、
det(sE − A) ̸≡ 0 とする．このとき伝達関数は
G(s) = D + C(sE − A)−1 B (2)
となる．E が正則でなければ (sE − A)−1 は一般にプ
ロパーでないので、伝達関数は非プロパーになる可能
性がある．
3
伝達関数の計算手法
ディスクリプタシステムの伝達関数は式 (2) で表現
できるが、ペンシル (sE − A) の逆行列の計算方法が異
なるいくつかの数値計算法がある．以下では、それら
の手法を誤差と計算量の観点で考察する．
3.1 ワイヤストラス標準形
本節ではディスクリプタシステムの伝達関数の係数
と関係が深いワイヤストラス標準形について説明する.
一般化固有値問題から得られる有限周波数に対する
一般化固有ベクトルからなる行列 Vs ∈ Rn×r と無限周
波数に対する一般化固有ベクトルからなる行列 Vf ∈
Rn×(n−r) を用いて
[
]−1
[
]
P T = EVs AVs
,
Q = Vs Vf
(3)
を定義する．これらを用いるとペンシル (sE − A) は
[
]
sI − As
0
T
P [sE − A]Q =
0
sAf − I
のように有限固有値に対するモードと無限固有値に対
応するモードに分解できる 3) . ただし As ∈ Rr×r は有
限周波数に対する一般化固有値からなるジョルダン形
式の行列で、Af ∈ R(n−r)×(n−r) は無限周波数に対す
る一般化固有値に対するジョルダン形式のベキ零行列
である.
式 (3) で定義した P と Q による変換と状態変数変換
[
]
[
] ηs
x(t) = Qη(t) = Vs Vf
ηf
を行うとワイヤストラス標準形
[
][
] [
][
] [
]
I 0
η˙s
As 0
ηs
Bs
=
+
u
0 Af
η˙f
0 I
ηf
Bf
[
]
[
] ηs
+ Du
y = Cs Cf
ηf
が得られる．
このとき G(s) は
ならば
G(s) = Cs (sI − As )−1 Bs + D − Cf (I − sAf )−1 Bf
を満たす ν∞ だけの無限固有ベクトルを選択する. ν∞
は無限固有値に対する幾何学的重複度であり、このベ
k
クトル v∞j
に対して漸化式
1
Ev∞j
= 0,
のように計算できる. ここで Af がベキ零行列であるこ
とから、
(I − sAf )
−1
2
= I + sAf + s
A2f
+ ··· + s
n−r−1
ここで Af はベキ零行列である．
3.2 一般化固有値問題に基づく手法
3.1 節で示した行列 P と Q を求めるには、一般化固
有値、すなわち
rank(λE − A) < n
を満たす λ ∈ C を求める必要がある．一般化固有値を
λ1 , λ2 , . . . , λl とするとき、ペンシルの階数が
rank(λi E − A) = n − νi ,
νi > 0
であるとする. ただし、νi は一般化固有値 λi に対する
幾何学的重複度である．νi だけ一般化固有ベクトルを
1
(λi E − A)vij
= 0,
j = 1, 2, . . . , νi
を満たすように選ぶ. ただし λi の代数的重複度を σi
とすると、σi > νi の場合には一意に定まらないので漸
化式
k−1
k
(λi E − A)vij
= Evij
,
k−1
k
Ev∞j
= −Av∞j
,
An−r−1
f
である．
また完全なワイヤストラス標準形でなくても以下に
定義する擬似ワイヤストラス形式に変形すれば上の計
算で伝達関数を求めることができる．
[
][
] [
][
] [
]
I 0
η˙s
ηs
Bs
As 0
=
+
u
η˙f
ηf
0 Af
0 I
Bf
[
]
] ηs
[
y = Cs Cf
+ Du
ηf
k = 2, 3, . . . , µij
j = 1, 2, . . . , ν∞
k = 2, 3, . . . , µ∞j
によって無限拡張ベクトルを求める必要がある. この
k
ようにして求まったベクトル群 {v∞j
}j,k は無限固有値
に対応する線形部分空間の基底を構成しており、その
次元は
ν∞
∑
n−r =
µ∞j
j=1
となる. ここで
µν
1
2
Vf = [v∞1
v∞2
· · · v∞νl l ]
である.
この手法において固有値の代数的重複度を求める際、
固有値が重複しているかいないかの判定を計算機で行う
のは丸め誤差等が含まれるので非常に困難である．また
拡張固有ベクトルを求める際、計算過程で行列の rank
と Kernel を計算する．一般に、行列の rank と Kernel
を計算するには特異値分解を行い、特異値が 0 である
かどうかの判定を行う必要があり、この判定によって
rank と Kernel の値が変わるよってこの手法は計算機
で行うには適していない.
3.3 一般化ファディーブのアルゴリズムを用いる手法
ディスクリプタシステム表現から伝達関数を計算す
る scilab の関数 des2tf では一般化ファディーブのアル
ゴリズムが用いられている 4) ．状態空間表現からプロ
パーな伝達関数を求める際に (sI − A) の逆行列の係数
を求めるときに用いられるファディーブのアルゴリズ
ムを拡張した (sE − A) の逆行列の係数を求めることが
できるアルゴリズムを次に示す.
アルゴリズム
det(sE − A) ̸= 0 のとき
を繰り返し解いて拡張固有ベクトルを求める必要があ
k
}i,j,k
る. このようにして求まったベクトルの集合 {vij
は有限固有値 {λi }i に対応する線形部分空間の基底を
構成しており、その次元は
ただし
νi
l ∑
∑
Bs
= B0 + sB1 + · · · + sn−1 Bn−1
Bf
= sB1 + s2 B2 + · · · + sm−1 Bm−1 + sm Bm
r=
µij
i=1 j=1
(sE − A)−1 =
a
となる. ここで
Vs = [v1 v2 · · · vνl ]
Vi = [vi1 vi2 · · · viνi ]
′
′
Bs
− Bf
a
′
= α0 + sα1 + · · · + sn−1 αn−1 + sn
このとき
(sE − A)−1 = · · · +
Si
+ Pi + sDi + · · ·
s
µν
1
2
Vij = [vij
vij
· · · vij l ]
である.
このとき
ν∞ = n − rankE > 0
とすると
Bn−1
= Si ,
αn−1
= −tr(ESi ASi EBn−1 )
′
Bn−2
=
αn−2
=
Si ASi EBn−1 + αn−1 ,
1
− tr(ESi ASi EBn−2 )
2
分解した後、Householder 変換または Givens Rotation
を用いれば簡単に行うことができる 7) .
..
.
Bn−k
αn−k
=
=
Si ASi Bn−k+1 + αn−k+1 ,
1
− tr(ESi ASi EBn−k )
k
..
.
B1
=
α1
=
B0
=
α0
=
′
B1
′
B2
′
B3
′
Bm
=
Si ASi EB2 + α2 ,
1
−
tr(ESi ASi EB2 )
n−1
Si ASi EB1 + α1 ,
1
− tr(ESi ASi EB1 )
n
Pi ADi
=
B1 ADi
=
..
.
B2 ADi
=
Bm−1 ADi
A1 = P1 AQ1 =
′
···
..
.
..
.
0
..
.
..
.
0
···
..
···
.
..
..
.
···
···
× ···
.
× ..
..
.
0
.. . .
.
.
0 ···
×
..
.
..
.
..
.
.
0
···
..
..
×
···
.
.
0
..
.
×
×
..
.
..
.
..
.
×
′
ここで計算量は O(n3 ) である．
′
この手法では Si , Pi , Di を求めるために一般化固有空
間の基底を求める必要がある. その計算を計算機で行
う場合誤差が発生する場合がある.
3.4 零点と極を用いた手法
ディスクリプタシステム表現から伝達関数表現に変
換する Matlab の関数 tf では零点と極を計算する手法
が用いられている 5) ．ディスクリプタ表現の係数行列
A, B, C, D, E から各入出力における零点と極を求める
ことができれば伝達関数表現に変換することができる．
零点と極を求めるということは一般化固有値を求め
るということと等しい．よって数値計算上安定しない
場合がある．
3.5 実シュア分解を用いた手法
実シュア分解を用いてディスクリプタシステム表現
から伝達関数を求める手法が提案されている 6) ．本論
文で提案する手法はこの手法を改善したものである．本
手法との違いは、ペンシルに対して実シュア分解を行
うという点である．
実シュア分解を行う際にペンシルを一般化ヘッセン
ベルグ形式に変形する. そこから A を擬似上三角行列
にする際に解が収束しないことがある．また収束した
としてもその計算を行うために本論文で提案する手法
より計算量が増え、誤差が発生する.
4
E1 = P1 EQ1 =
×
提案手法
4.2
無限固有値成分の入れ替え
すべての無限固有値成分を右下に集めることを考え
る．すなわち、係数行列 E の対角成分の 0 を右下に押し
下げ、それに対応する係数行列 A の成分を上三角化し、
E2 = P2 E1 Q2 =
A2 = P2 A1 Q2 =
×
0
..
.
..
.
..
.
..
.
0
×
..
.
..
.
...
···
..
.
..
.
..
.
...
× ··· ···
.
× ..
.. ..
.
.
0
.. . .
. ×
.
..
..
.
.
..
.
0 ... ...
···
···
···
×
..
.
..
.
..
.
..
.
..
.
..
.
×
..
..
. 0
.
..
..
.
. ×
... ... 0 0
··· ··· ··· ×
..
.
..
.
..
×
.
..
..
.
0
.
..
.. ..
..
.
.
. .
... 0
0 ×
本章では本論文で提案する手法について述べる．ディ
スクリプタシステムを完全なワイヤストラス標準形に
変形しなくても擬似ワイヤストラス形式に変形すれば
伝達関数を計算することができる．その形に変形するた
めにペンシルを一般化ヘッセンベルグ形式に変形する．
のように変換できること以下に示す．
4.1 一般化ヘッセンベルグ形式への変形
ディスクリプタシステム表現の係数行列 E と A を一般
化ヘッセンベルグ形式に変換する．この変換は E を QR
Givens Rotation を用いて E と A の構造を崩さない
ようにして入れ替えを行う手順を以下に示す．ここで、
すでに集められた無限固有値成分が一つで E(k, k) = 0
のように A(k + 2, k) に 0 でない成分が現れる．A を
再びヘッセンベルグ形式にするために j = k として
A(k + 2, k) = 0 とすると
の場合を考える．
E=
A=
k
× ··· ··· ··· ··· ···
..
.
0
.. . .
. ×
.
..
..
. 0
.
..
..
. ×
.
..
.. ..
.
.
.
..
..
.
.
0 ··· ··· ··· ··· ···
×
×
0
..
.
..
.
..
.
..
.
0
···
..
.
..
.
..
.
···
..
..
..
···
···
.
.
..
.
.
..
.
..
.
.
..
.
..
..
···
···
···
.
···
···
×
..
.
..
.
..
.
..
.
..
.
× ×
0 0
··· · ×
..
.
..
.
..
.
..
.
..
..
.
.
× × ×
0 0 ×
E
A = AZk,k+1 = Hessen1
k
ここで対角成分が X である上三角行列を Trig(X)、上
ヘッセンベルグ行列で p 行から n 行が上三角化されて
いる行列を Hessenp と定義すると E と A は
E
=
A =
= EZk,k+1 = Trig(×, · · · , ×, ×, 0, 0, ×, · · · , ×, 0)
Trig(×, · · · , ×, 0, ×, · · · , ×, 0)
Hessen1
となる．
2 × 2 の Givens Rotation 行列を Qi,i+1 とする．E の
i と i + 1 行に対して左からかけて E(i + 1, i + 1) = 0
とする．また Givens Rotation 行列 Zj,j+1 は A の j と
j + 1 列に対して右からかけて A(i + 1, j) = 0 とする．
まず i = k として E(k + 1, k + 1) = 0 とすることを
考える．このとき A に対しても同様の変換行列をかけ
ると A(k + 1, k − 1) に 0 でない成分が現れる．p 行 q
列の成分が 0 でなく, それ以外の成分が 0 である行列を
Yp,q とすると
E
= Qk,k+1 E = Trig(×, · · · , ×, 0, 0, ×, · · · , ×, 0)
A
= Qk,k+1 A = Hessen1 + Yk+1,k−1
となる．A を再びヘッセンベルグ形式にするために j =
k − 1 として A(k + 1, k − 1) = 0 とすると
となる．E の対角成分に 0 が現れるまで上と同様の変
換を i を増加させて行うと対角成分に連続する 2 つの
0 がスライドする．最終的に
E
= Trig(×, · · · , ×, 0, 0, 0)
A
= Hessen1
のように E と A を変換できる．このとき E の n 行目
の成分がすべて 0 であるので右からどんな変換行列を
かけても E の n 行目には影響がないので、j = n − 1
として A(n, n − 1) = 0 とすると
E
= EZn−1,n = Trig(×, · · · , ×, 0, 0)
A
= AZn−1,n = Hessen2
となる．E の右下に 0 を押し下げるのと同時に対応す
る A の成分を対角化すること、つまり無限固有値成分
を右下に集めることができた．
ここでたかだか n2 /4 回入れ替えが必要である．一回
の入れ替えで 8n 回の計算が必要である．よって計算量
は O(n3 ) である．
4.3 ペンシルのブロック対角化
A2 と E2 を
[
]
E11 E12
= E2
0
E22
[
]
A11 A12
= A2
0
A22
のように分割する．このとき A11 は上ヘッセンベルグ
行列、E22 はベキ零行列となっている．A22 , E11 は上
三角行列で対角成分がすべて 0 でないので正則である．
ここで
[
]
[
]
Ir
L
Ir
R
P3 =
, Q3 =
0 In−r
0 In−r
とし、
[
E11
P3
0
[
A11
P3
0
E12
E22
A12
A22
]
[
Q3 =
]
[
Q3 =
E11
0
0
E22
A11
0
0
A22
]
= E3
]
= A3
を満たす L と R を求めることを考える．
E
= EZk−1,k = Trig(×, · · · , ×, 0, 0, ×, · · · , ×, 0)
E11 R + LE22 = −E12
(4)
A
= AZk−1,k = Hessen1
A11 R + LA22 = −A12
(5)
となる．A と E の構造を崩すことなく E(k +1, k +1) =
0 とすることができた．
次に i = k + 1 として E(k + 2, k + 2) = 0 とするこ
とを考える．このとき先ほどと同様に
E
=
A =
Qk+1,k+2 E = Trig(×, · · · , ×, 0, 0, 0, ×, · · · , ×, 0)
Qk+1,k+2 A = Hessen1 + Yk+2,k
式 (4),(5) はクロネッカ積を用いて以下のように表すこ
とができる.
[
][
] [
]
T
Ir ⊗ E11 E22
⊗ In−r
vec(R)
−vec(E12 )
=
Ir ⊗ A11 AT22 ⊗ In−r
vec(L)
−vec(A12 )
この線形方程式を解くと L と R を求めることができ
る．このとき計算量は O((r(n − r))3 ) である．
4.4 擬似ワイヤストラス形式への変形
A22 と E11 を単位行列にすることを考える．
[ −1
]
E11
0
P4 =
0
A−1
22
を左からかけると
[
P4 E3 =
[
P4 A3 =
I
0
0
A−1
22 E22
−1
A22
E11
0
0
I
式 (6) から正しく伝達関数を計算できていることが確認
できる参考のため実行環境が CPU AMD Phenom II X4
945 3.00GHz, メモリ 3.87GB, OS Windows7(64bit)、
Java JDK1.6 23 で実行時間を計測したところ 75ms で
あった．
また Fig.1 に示したディスクリプタシステム表現か
ら scilab version 5.3.0 の des2tf 関数によって伝達関数
を求めた結果を以下に示す．
✓
✏
]
]
警告:
行列は単数に近いです、または、最悪の比率です。
状態 = 0.0000D+00
ans =
column 1
2
3
3.632D-08-1.0D-07s+0.9999998s-3.130D-09s
----------------------------------------2
-3.632D-08+0.9999998s+s
2
10000000s+s
------------1
column 2
1.000D-07s
----------------------2
-3.632D-08+0.9999998s+s
1
1
となる. ここで A22 は上三角行列であるので A−1
22 は上三
角行列である．E22 はベキ零行列であるので A−1
22 E22 は
ベキ零行列である. よってペンシルを目的の形に変形で
きたことが確認できる．P = P4 P3 P2 P1 , Q = Q1 Q2 Q3
として与えられたシステムのディスクリプタ表現に対
して左から P 、右から Q をかけると伝達関数が計算で
きる形に変形できる.
例題による評価
5
5.1 例題 1: 他のツールとの比較
Fig.1 のシステムについて考える．このシステムの
ディスクリプタシステム表現を制御系のモデリングプ
ラットフォーム Jamox8) で求めると Fig.1 の結果が得
られる．
✒
✑
警告が出て正しく計算できていないことが確認できる．
これは G3 のブロックの極が 1 と 10000000 で大きく異
なるため一般化固有ベクトルの基底を正しく求めるこ
とができなかったためであると考えられる．また上と
で実行したところ、実行時間は 237ms であった．
また Fig.1 に示したディスクリプタシステム表現か
ら Matlab version R2007b の tf 関数によって伝達関数
を求め、極零相殺を行った結果を以下に示す．
✓
✏
零点/極/ゲイン from input 1 to output...
s
#1: ----(s+1)
#2:
Fig. 1: Improper System(Jamox)
このシステムの伝達関数は
[
G(s) =
s
s+1
s2 +10000000s
1
1
s2 +10000001s+10000000
1
1
s (s+1e007)
零点/極/ゲイン from input 2 to output...
(s+3.638e-005)
#1: ----------------s (s+1e007) (s+1)
]
(6)
#2:
1
✒
✑
in2 から out1 の伝達零点が正しく求めれてないことが
確認できる．また上と同じ環境で実行したところ、実
行時間は 31ms であった．
である．
次に提案手法によって伝達関数を求めた結果を以下
に示す．
✓
✏5.2 例題 2: 実システムへの適用
Fig.2 に 3 リンク平面マニュピレーターのモデル図を
===
(
2 x
2) RaMatrix ===
示す．このモデルはビルの窓ガラスの清掃ロボットに
[
( 1)
] [
( 2)
]
s
1
用いられ、接面が A から B の間をスライドするように
( 1) --------------- ----------------------各トルクに入力電圧を加える．運動方程式は以下のよ
s+1
s^2+10000001*s+10000000
うになる．
s^2+10000000*s
( 2) --------------1
✒
1
----------------------1
˙ + Gz (θ) = uz + F T µ
Mz (θ)¨
z + Cz (θ, θ)
z
✑
ψz (θ) = 0
✓
✏
===
ans
(
3
x
3)
RaMatrix ===
( 1)
]
0.0144102
1) -------------------------------------s^2+0.103232*s+0.0529961
[
(
(
0.10211*s^2-5.14222E-08*s+1.00923
2) -------------------------------------s^2+0.103232*s+0.0529961
(
0.0720509*s^2+1.00871E-06
3) -------------------------------------s^2+0.103232*s+0.0529961
[
Fig. 2: Three-Link Planar Manipulator
ただし
z = [x y ϕ]
T
ϕ = θ1 + θ2 + θ3
ここで uz ∈ R3 はトルクによる制御入力、Fz =
δψz (z)/δz 、µ ∈ R2 はラグランジュの未定乗数 、
˙ ∈ R3 は
Mz (θ) ∈ R3×3 は質量マトリクス、Cz (θ, θ)
遠心力とコリオリ力によるベクトル、Gz (θ) は重力に
よるベクトル、ψz (θ) は拘束条件を表している．
これを
zω = [l l +
M0 δ¨
z + D0 δ z˙ + K0 δz = S0 δu + F0T δµ
F0 δz = 0
ただし、δz, δ z,
˙ δ¨
z は平衡点からの差である．このモデル
は以下のようなディスクリプタシステム表現で表せる．
E x˙ = Ax + Bu
I
E= 0
0
0
M0
0
y = Cx + Du
0
0
0 , A = −K0
0
F0
0
0 1
B = S0 , C = 0 0
0
0 0
0
0
0
0
0
0
0
F0T ,
0
I
−D0
0
0
0
0
0
0
0
0
1
0
0
0 ,
1
D=0
これらの係数行列から本手法を用いて伝達関数を求め
た結果を以下に示す．ただし、上式の導出方法と各パ
ラメータは参考文献を参照 1) ．
本手法は実システムに対しても有効であり、伝達関
数表現でモデルを表すことによって入出力関係が明白
になっている.
(
1)
(
2)
0.549093*s^2+0.0761144*s-1.82121
-----------------------------------s^2+0.103232*s+0.0529961
(
3)
- 0.13281*s^2-1.85934E-06
----------------------------------s^2+0.103232*s+0.0529961
[
(
∆l T
∆l T
0] , zω = [l l +
0]
2
2
で線形近似する．ただし ∆l は AB 間の距離である．こ
のとき以下の線形方程式を得る．
( 2)
]
-0.0265619
-----------------------------------s^2+0.103232*s+0.0529961
(
-0.651204*s^2-0.0761145*s+0.811986
2) --------------------------------------s^2+0.103232*s+0.0529961
(
-0.939242*s^2-0.103233*s-0.0529953
3) --------------------------------------s^2+0.103232*s+0.0529961
✒
6
( 3)
]
0.0121518
1) --------------------------------------s^2+0.103232*s+0.0529961
✑
おわりに
本論文では、ディスクリプタシステムの伝達関数を
一般化ヘッセンベルグ形式を用いて求めるための数値
計算法を提案した．本手法は一般化固有値や一般化固
有ベクトルを求めないため誤差が発生しにくく高速で
ある．提案した手法と以前から存在する手法の比較を
行った．
参考文献
1) Guang-Ren Duan. Analysis and Design of Descriptor
Linear Systems. Springer, 2010.
2) 阿南悟, 古賀雅伸. 非プロパーなシステムを含む結合線形
システムのディスクプリタ形式を用いた数式状態空間実
現. 第 38 回制御理論シンポジウム, pp. 167–170, 2009.
3) 汐月哲夫. 表現方法としてのディスクリプタシステム―
その可能性を探る. システム/制御/情報, Vol. 39, No. 5,
pp. 219–224, 1995.
4) INRIA. Scilab,
http://www.scilab.org/.
5) Inc. The Math Works. Matlab,
http://www.mathworks.com/products/matlab/.
6) Andra Varga. Numerical algorithm and software tools
for analysis and modeling of descriptor systems. prepr.
of 2nd IFAC Workshop on System Structure and Control, pp. 392–395, 1992.
7) Gene H. Golub and Charles F. Van Loan. Matrix
Computations Third Edition. Johns Hopikins Univ Pr,
1996.
8) 石倉雄飛. Eclipse rcp を用いた制御系開発支援統合環
境 jamox の開発. 第 54 回システム制御情報学会, pp.
585–586, 2010.
© Copyright 2024 ExpyDoc
