WO2005015737A1

WO2005015737A1 - システム推定方法及びプログラム及び記録媒体、システム推定装置

Info

Publication number: WO2005015737A1
Application number: PCT/JP2004/011568
Authority: WO
Inventors: Kiyoshi Nishiyama
Original assignee: Japan Science And Technology Agency
Priority date: 2003-08-11
Filing date: 2004-08-05
Publication date: 2005-02-17
Also published as: JPWO2005015737A1; EP2560281A2; EP1657818A1; EP2560281A3; EP1657818B1; US20070185693A1; JP4444919B2; CN1836372A; CN1836372B; US7480595B2; EP2560281B1; EP1657818A4

Abstract

忘却係数を理論的に最適に決定できる推定方法を確立すると共に、その数値的に安定な推定アルゴリズムと高速アルゴリズムを開発する。まず、処理部は、記憶部又は入力部ら上限値γfを読み出し又は入力する(S101)。処理部は、式(15)によって忘却係数ρを決定する(S103)。その後、処理部は、忘却係数ρに基づき、式(10)～式(13)のハイパーH∞フィルタを実行する(S105)。処理部101は、式(17)(あるいは、後述の式(18))の存在条件を計算し(S107)、その存在条件がすべての時刻で満たされれば(S109)、γfをΔγだけ小さくして同じ処理を繰り返す(S111)。一方、あるγfで存在条件を満たさなくなったとき(S109)、そのγfにΔγを加えたものをγfの最適値γf opとして出力部に出力及び/又は記憶部に記憶する(S113)。

Description

明細書システム推定方法及びプログラム及び記録媒体、システム推定装置技術分野

本発明は、システム推定方法及びプログラム及び記録媒体、システム推定装置に係り、特に、 H_∞評価基準に基づいて開発されたハイパー H_∞フィルタの高速 H_∞フィルタリングアルゴリズムを用いて、状態推定のロバスト化と忘却係数の最適化を同時に実現するシステム推定方法及びプログラム及び記録媒体、システム推定装置に関する。

背景技術

一般に、システム推定とは、入出力データに基づいてシステムの入出力関係の数理モデル (伝達関数、あるいはインパルス応答など)のパラメータを推定することである。代表的な応用例として、国際通信におけるエコーキャンセラ、データ通信における自動等化器、音響システムにおけるエコーキャンセラや音場再生および自動車などにおけるアクティブ騒音制御などがある。詳細は、非特許文献 1 . 1 993年電子情報通信学会「ディジタル信号処理ハンドブック」等参照。 (基本原理）

図 8に、システム推定のための構成図の例を示す（未知システムは IIR (Infinite Impulse Response)フィレタで表現してもよしゝ） ₀

このシステムは、未知システム 1、適応フィルタ 2を備える。また、適応フィルタ 2は、 FIRディジタルフィルタ 3、適応アルゴリズム 4を有する。

以下に、未知システム 1を同定する出力誤差方式の一例を説明する。ここで、 ujま未知システム 1の入力、 d_kは所望信号であるシステムの出力、 cT_kはフィルタの出力である。（なお、」は、推定値の意味であり、文字の真上に付されるものであるが、入力の都合上文字の右上に記載する。以下同様。 )

未知システムのパラメータとしては、一般にインパルス応答が用いられるので、適応フィルタは図の評価誤差 e_k = d_k— cT_kを最小にするように適応アルゴリズムによって FI Rディジタルフィルタ 3の係数を調節する。

また、従来、システムのパラメータ (状態）の推定には、誤差共分散行列の更新式 (リカツチ方程式）に基づくカルマンフィルタが広く用いられて来た。詳細は、非特許文献 2. S. Haykin: Adaptive filter theory, Prentice— Hall(1996)など Ίこ示されている。

以下にカルマンフィルタの基本原理について説明する。

次式のように、状態空間モデルで表された線形システム

x_k+1 = p~^1/2x_k, y_k=H_kx_k+v_k (1)

| _kは、状態の誤差共分散行列∑ | ，を用いて次のように得られる。

ただし、 ά₀ト i

= ε₀Ι, ε₀ > 0 (5) x_k:状態ベクトルまたは単に状態；未知であり、これが推定の対象となる。

y_k:観測信号；フィルタの入力となり、既知である。

H_k:観測行列；既知である。

v_k:観測雑音；未知である。 P：忘却係数；一般に試行錯誤で決定される。

K_k :フィルタゲイン；行列∑ I _k から得られる。

∑

| _kの誤差の共分散行列に対応；リカツチ方程式により得られる。

∑ \_{+1 ] k} : x i | _kの誤差の共分散行列に対応；リカツチ方程式により得られる。 ∑ _{| 0}：初期状態の共分散行列に対応；本来未知であるが、便宜上 ε ₀Ιが用いられる。

また、本発明者は、既に高速 Η_∞フィルタによるシステム同定アルゴリズムを提案した（特許文献 1参照）。これは、システム同定のために新たに Η_∞評価基準を定め、これに基づくハイパー Η_∞フィルタの高速アルゴリズムを開発すると共に、この高速 1~1_∞フィルタリングアルゴリズムに基づく高速時変システム同定方法である。高速 Η_∞フィルタリングアルゴリズムは、単位時間ステップ当たり計算量 Ο (Ν)で急激に変化する時変システムの追跡が可能である。上限値の極限で高速カルマンフィルタリングアルゴリズムと完全に一致する。このようなシステム同定により時不変および時変システムの高速実時間同定および推定を実現することができる。

なお、システム推定の分野で通常知られる方法として、例えば、非特許文献 2、 3参照のこと。

(エコーキャンセラへの適用例） '

国際電話など長距離電話回線では，信号増幅などの理由から 4線式回線が用いられている。一方、加入者回線は比較的短距離なので、 2線式回線が使用されている。図 9に通信系とエコーについての説明図を示す。 2線式回線と 4線式回線の接続部には図示のようにハイブリッドトランスが導入され、インピーダンス整合が行われている。このインピーダンス整合が完全であれば、話者 Βからの信号 (音声)は話者 Αのみに到達する。し力、し、一般に整合を完全とするのはむずかしく、受信信号の一部は 4 線式回線に漏れ、増幅された後、再び受信者 (話者 A)に戻ると云った現象が起こる。これがエコー（echo)である。エコーは、伝送距離が長くなるにつれて (遅延時間が長くなるにつれて)影響が大きくなリ、著しく通話の品質を劣化させる (パルス伝送においては近距離であってもエコーによる通話品質の劣化は大きく影響する)。

図 1 0に、エコーキャンセラの原理図を示す。

そこで、図示のようにエコーキャンセラ (echo canceller) を導入し、直接観測可能な受信信号とエコーを用いてエコーパスのインパルス応答を逐次推定し、それを利用して得た疑似エコーを実際のエコーから差し引くことによってエコーを打ち消し、その除去を図っている。

エコーパスのインパルス応答の推定は、残留エコー e_kの平均 2乗誤差が最小になるように行われる。このとき、エコーパスの推定を妨害する要素は、回線雑音と話者 Aからの信号 (音声)である。一般に、話者 2人が同時に話し始めた (ダブルトーク)ときはィンパルス応答の推定を中断する。また、ハイブリッドトランスのインパルス応答長は 50 [ms]程度なので、サンプリング周期を 1 25 [〃s]とするとエコーパスのインパルス応答の次数は実際は 400程度となる。

非特許文献 1

1 993年電子情報通信学会「ディジタル信号処理ハンドブック」

非特許文献 2

S. Haykin : Adaptive filter theory, Prentice— Hall O 996)

非特許文献 3

B. Hassibi, A. H. Sayed, and T. Kailath: Indefinite-Quadratic Estimation and Control", SIAM (1 996)

特許文献 1

特開 2002— 1 351 7 1号公報

発明の開示

し力、しなカら、式（1 )〜（5)のような従来の忘却係数 pを入れたカルマンフィルタでは、忘却係数 Pの値は試行錯誤で決定しなければならず非常に手間が掛かった。さらに、決定された忘却 P係数の値が果たして最適な値であるかどうか判定する手段も無かった。

また、カルマンフィルタで用いる誤差共分散行列は、本来、零でない任意のベクトノレとの 2次形式が常に正（以下、「正定」という。）であるがコンピュータにより単精度で計算した場合にはその 2次形式が負（以下、「負定」という。）となり、数値的に不安定になることが知られている。また、計算量が 0 (N²) (あるいは 0 (N³) )であるため、状態ベクトル x_kの次元 Nが大きい場合、 1ステップ当たりの演算回数が急激に増大し、実時間処理には適さなかった。

本発明は、以上の点に鑑み、忘却係数を理論的に最適に決定できる推定方法を確立すると共に、その数値的に安定な推定アルゴリズムと高速アルゴリズムを開発することを目的とする。また、本発明は、通信システムや音響システムにおけるエコーキヤンセラ、音場再生又は騒音制御などに適用することができるシステム推定方法を提供することを目的とする。

本発明は、上記課題を解決するために、新たに考案した H_∞最適化手法を用いて忘却係数が最適決定可能な状態推定アルゴリズムを導出する。さらに、常に正定であるべき誤差共分散行列の代わりに、その因数行列を更新することによって数値的に安定な推定アルゴリズムと高速アルゴリズムを開発する。

本発明の第 1の解決手段によると、

次式で表される状態空間モデルに対して、

x_{k+ 1} = F_kx_k + G_kw_k

y_k= H_kx_k+v_k

z_k= H_kx_k

ここで、

x_k :状態ベクトルまたは単に状態

w_k :システム雑音

v_{k :}観測雑音 y_k :観測信号

z_k :出力 1目"^

F_k：システムのダイナミックス

G_k :駆動行列

評価基準として忘却係数 Pで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値 r_fに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数 pの最適化を同時に行うためのシステム推定方法及びプログラム及び該プログラムを記録したコンビユータ読み取り可能な記録媒体であって、処理部は、上限値 r_f、フィルタの入力である観測信号 y_k、観測行列 H_kを含む値を記憶部又は入力部から入力するステップと、処理部は、前記上限値？^に従い、状態空間モデルに関連する忘却係数 pを決定するステップと、

処理部は、記憶部から初期値又はある時刻の観測行列 H_kを含む値を読み取り、前記忘却係数; 0を用いて次式で表されるハイパー H_∞フィルタを実行するステップと、

X k I k一 f k- 1 ^X k - 1 I k - 1 + ks, yk一 HkF_k一 i X _k- l | k- 1

ここで、

観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

F_k：システムのダイナミックス

K_S, _K :フィ (レタゲイン

処理部は、ハイパー H_∞フィルタに関する求められた値を記憶部に記憶するステップと、処理部は、観測行列とフィルタゲイン K_s,；により、前記上限値 r_f及び前記忘却係数) 0に基づく存在条件を計算するステップと、処理部は、上限値 r_fを小さくしていき前記ハイパー H_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、

を含む前記システム推定方法、各ステツラをコンピュータに実行させるためのシステム推定プログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。

また、本発明の第 2の解決手段によると、

次式で表される状態空間モデルに対して、

Xk+ i = F_kx_k + G_kw_k

y_k= H_kx_k + v_k

z_k = H_kx_k

ここで、

x_k :状態ベクトルまたは単に状態

w_k :システム雑音

v_{k :}観測雑音

y_k :観測信号

z_{k :}出力信号

F_k :システムのダイナミックス

G_{k :}駆動行列

評価基準として忘却係数 Pで重み付けされた外乱からフィルタ誤差への最大エネルギ一ゲインを予め与えられた上限値 r _fに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数 Pの最適化を同時に行うためのシステム推定装置であって、

推定アルゴリズムを実行する処理部と、

前記処理部により読み取り及び又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、

を備え、前記処理部は、上限値 _f、フィルタの入力である観測信号 y_k、観測行列 H_kを含む値を記憶部又は入力部から入力すること、前記処理部は、前記上限値 r_fに従い、状態空間モデルに関連する忘却係数 Pを決定すること、

前記処理部は、記憶部から初期値又はある時刻の観測行列 H_kを含む値を読み取リ、前記忘却係数; 0を用いて次式で表されるハイパー H_∞フィルタを実行すること、

X k I k^— F k- 1 ^X k - 1 I k - 1 +^3、 1一 hkF_k— ·] Χ k - l | k- 1

ここで、

|_k :観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

F_k :システムのダイナミックス

K_{s k} :フィルタゲイン

前記処理部は、ハイパー H_∞フィルタに関する求められた値を記憶部に記憶するこ

前記処理部は、観測行列とフィルタゲイン K_s,；により、前記上限値 r _f及び前記忘却係数 )に基づく存在条件を計算すること、前記処理部は、上限値 r_fを小さくしていき前記ハイパー H_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶すること

を備えた前記システム推定装置が提供される。

本発明の推定方法は忘却係数を最適に決定することが可能であり、かつアルゴリズムは単精度でも安定に動作可能であるため、低コストで高い性能が実現できる。一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。図面の簡単な説明

図 1は、本実施の形態に関するハードウヱァの構成図である。

図 2は、 H_∞フィルタの口バス卜化と忘却係数 pの最適化についてのフローチャートで c6る。

図 3は、図 2中の H_∞フィルタ（S 1 05)のアルゴリズムについてのフローチャートである。

図 4は、定理 2の平方根アレアルゴリズムの説明図である。

図 5は、定理 3の数値的に安定な高速アルゴリズムのフローチャートである。図 6は、インパルス応答山₌。²³の値を示す図である。

図 7は、定理 3の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果である。

図 8は、システム推定のための構成図である。

図 9は、通信系とエコーについての説明図である。

図 1 0は、エコーキャンセラの原理図である。

発明を実施するための最良の形態

以下に、本発明の実施の形態について説明する。

1 .記号の説明

まず、本発明の実施の形態で用いる主な記号及びその既知又は未知について説明する。

x_k :状態ベクトルまたは単に状態；未知であり、これが推定の対象となる。

X。：初期状態；未知である。

w_k :システム雑音；未知である。

v_k :観測雑音；未知である。 y_k:観測信号;フィルタの入力となり、既知である。

z_k:出力信号；未知である。

F_k:システムのダイナミックス；既知である。

G_k:駆動行列；実行時に既知となる。

H_k:観測行列；既知である。

観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値；フィルタ方程式によつて与えられる。

x _{+1 |k}：観測記号 y。〜y_kまで用いた時刻 k+1の状態 x_k+1の推定値；フィルタ方程式によって与えられる。

χ^Λ ₀|₀:状態の初期推定値；本来未知であるが、便宜上 0が用いられる。

誤差の共分散行列に対応；リカツチ方程式によって与えられる。

I _kの誤差の共分散行列に対応；リカツチ方程式によって与えられる。

∑ ！ 0:初期状態の共分散行列に対応；本来未知であるが、便宜上 ε。1が用いられる。

K_s, _k：フィルタゲイン；行列∑ I ，から得られる。

P：忘却係数;定理"!〜 3の場合、 "T_fが決まれば) 0=1—ズ（： T_f)より自動的に決定される。 e_f._i:フィルタ誤差 R_e,_k:補助変数

なお、記号の上に付される" ^Λ"、 V'は、推定値の意味である。また、 "〜"、 "一"、 " U"等は、便宜上付加した記号である。これらの記号は、入力の都合上、文字の右上に記載する力数式で示すように、文字の真上に記載されたものと同一である。また、 x、w、H、G, K、R、 ∑、等は行列であり、数式で示すように太文字で記されるものであるが、入力の都合上、普通の文字で記載する。 2.システム推定のハードウェア及びプログラム

本システム推定方法又はシステム推定装置'システムは、その各手順をコンビユータに実行させるためのシステム推定プログラム、システム推定プログラムを記録したコンピュータ読み取り可能な記録媒体、システム推定プログラムを含みコンピュータの内部メモリにロード可能なプログラム製品、そのプログラムを含むサーバ等のコンビュータ、等により提供されることができる。

図 1は、本実施の形態に関するハードウェアの構成図である。

このハードウェアは、中央処理装置（CPU)である処理部 101、入力部 102、出力咅 15103、表示部 104及び記憶部 105を有する。また、処理部 101、入力部 102、出力部 103、表示部 104及び記憶部 105は、スター又はバス等の適宜の接続手段で接続されている。記憶部 105は、システム推定される「1.記号の説明」で示した既知のデータが必要に応じて記憶される。また、未知 "既知のデータゃ訐算されたハイパ —H_∞フィルタに関するデータ'その他のデータが処理部 101により、必要に応じて書込み及び Z又は読み出しされる。

3.忘却係数が最適決定可能なハイパー H_∞フィルタ

(定理 1)

次式のような状態空間モデルを考える。 xk+i = F_kx_k + G_kw_k, w_k,x_k e TZ^N (6) Vk = H_kx_k + v_k, y_k,v_k e Tl (7) z_k = _kx_k, z_ke1l, H_keTZ^]xN, k = 0,l,...,L (8)

このような状態空間モデルに対して、次式のような H_∞評価基準を提案する。

sup く 7/

x^{o wi vi}}

(9) この H_∞評価基準を満たす状態推定値 x ！ _k (あるいは出力推定値 z^v _k , _K)は、次のレベル r _fのハイパー H_∞フィルタによって与えられる。

(10)

(11)

(12) ∑k\k

T

k^k\k^ k IP (13) ここで、

ef,i =

=∑o

P 0

(14) o -pyj J '

< p = 1 - x(j_f) < 1, 7/ > 1 (15) である。なお、式（ 1 1 )はフィルタ方程式、式（ 1 2)はフィルタゲイン、式（ 1 3)はリカツチ方程式をそれぞれ示す。

また、駆動行列 G_Kは次のように生成される。

Xilf)

また、上述の高速 H_∞フィルタの追従能力を向上するためには、上限値 r_fは次の存在条件を満たすように出来るだけ小さく設定する。

ただし、（r_f)は、； c ( i ) = i、（∞) =〇を満たす r_fの単調減衰関数である。定理 1の特徴は状態推定のロバスト化と忘却係数 pの最適化を同時に行っている点にある。

図 2に、 H_∞フィルタのロバスト化と忘却係数 yOの最適化についてのフローチャートを示す。ここで、ブロック「EXC>0」： H_∞フィルタの存在条件、

△ 7·:正の実数、である。まず、処理部 101は、記憶部 105又は入力部 102から上限値 r_fを読み出し又は入力する（S101)。この例では、 r_f>>1 が与えられる。処理部 101は、式（15)によって忘却係数 pを決定する（S103)。その後、処理部 101は、忘却係数 yOに基づき、式（10)〜式（13)のハイパー H_∞フィルタを実行する（S105)。処理部 101は、式 (17) (あるいは、後述の式（18))の右辺（これを EXCとする）を計算し（S107)、その存在条件がすべての時刻で満たされれば（sio9)、 r_fを△ rだけ小さくして同じ処理を繰り返す（si 11)。一方、ある r_fで存在条件を満たさなくなったとき（si 09)、その r_fに△ rを加えたものを r_fの最適値 r_f°^pとして出力部 103に出力及び z又は記憶部 105に記憶する（S113)。なお、この例では、△ rを加えているが、それ以外の予め設定された値を加えてもよし、。この最適化のプロセスを r一イタレーシヨンと呼ぶ c なお、処理部 101は、 H_∞フィルタ計算ステップ S105及び存在条件の計算ステップ S 107等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部 1 05に記憶し、また、記憶部 105から読み出すようにしてもよい。

ハイパー H_∞フィルタが存在条件を満たすとき、式（9)の不等式は常に満たされる。よって、式（9)の分母の外乱エネルギーが限定される場合、式（9)の分子の 2乗推定誤差の総和は有界となり、ある時刻以降の推定誤差が 0になる。これは、 "T_fをより小さく出来れば、状態 x_kの変化に推定

_kが速やかに追従できることを意味する。ここで、定理 1のハイパー H_∞フィルタのアルゴリズムは通常の H_∞フィルタのものとは異なることに注意されたし、。また、 "T_f→∞のとき、 p=1、G_k=0となり、定理 1の 1~1_∞フィルタのアルゴリズムはカルマンフィルタのアルゴリズムと一致する。

図 3に、図 2中の（ハイパー） H_∞フィルタ（S105)のアルゴリズムについてのフローチヤー卜を示す。ハイパー H_∞フィルタリングアルゴリズムは以下のように要約することができる。

[ステップ S201 ] 処理部 101は、記憶部 105から再帰式の初期条件を読み出し、又は、初期条件を入力部 102から入力し、図示のように定める。なお、 Lはあらかじめ定められた最大データ数を示す。

[ステップ S203] 処理部 101は、時刻 kと最大データ数 Lとを比較する。処理部 101 は、時刻 kが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。または、必要に応じて再スタ —卜を行ってもよい。）

[ステップ S205] 処理部 101は、フィルタゲイン K_s,_kを式（12)を用いて計算する。

[ステップ S207] 処理部 101は、式（11 )のハイパー H_∞フィルタのフィルタ方程式を更新する。

[ステップ S209] 処理部 101は、誤差の共分散行列に対応する項∑ _|k、 ∑ ₊₁₁ _kを式（13)のリカツチ方程式を用いて計算する。

[ステップ S211] 時刻 kを進ませて（k = k+1)、ステップ S203に戻り、データがある限り続ける。

なお、処理部 101は、 H_∞フィルタ計算ステップ S205〜S209等の各ステップで求めた適宜の中間値及び最終値、存在条件の値等を必要に応じて適宜記憶部" !05に記憶し、また、記憶部 105から読み出すようにしてもよい。

(スカラー存在条件）

ところで、式（17)の存在条件の判定には 0(N²)の計算量が必要であった。しかし、次の条件を用いれば計算量 O(N)で定理 1の1~)_∞フィルタの存在性、すなわち式（9)を検証することができる。

系 1 :スカラー存在条件

次の存在条件を用いれば計算量 O(N)でハイパー H_∞フィルタの存在性が判定できる。

+ f f > 0, i = 0，… , (18) で、

ただし、 K_{s i}は式（12)で求めたフィルタゲインである < (証明）

以下に系 1の証明を説明する。

2 2の行列 _kの特性方程式

= 0

を解けば、 R_{e k}の固有値 λ；が次のように得られる。

ただ

である ( もし、

であれば、行列 R_e, _kの 2つの固有値の 1つは正となり、もう 1つは負となり、行列 R_kと

R_e, _kは同じイナ一シャをもつ。これより、 τ Hkh._k HuK 一

を用いれば、式（18)の存在条件が得られる。ここで、 H_kK_{s k}の計算量は O(N)である。

4.数値的に安定な状態推定アルゴリズム

上述のハイパー H_∞フィルタは、

単位時間ステップ当たりの計算量は 0 ( N²)となる、すなわち、 N²に比例する算術演算が必要となる。ここで、 Nは状態ベクトル x_kの次元である。よって、 x_kの次元が增加するにつれて本フィルタの実行に要する計算時間は急速に増大する。また、誤差共分散行列は、その性質から常に正定でなければならないが、数値的には負定になる場合がある。特に、単精度で計算した場合はこの傾向は顕著となる。このとき、フィルタは不安定となることが知られている。よって、アルゴリズムの実用化および低コスト化のためには、単精度 (例： 32bit)でも動作可能な状態推定アルゴリズムの開発が望まれる。

そこで、次に、

^Ke, k^— e. k^Jl ^K k、 k I k-1 厶 k I k- 1 ^e- k | k- 1

に着目して、数値的に安定化した定理 1の1~1_∞フィルタ（平方根アレイアルゴリズム)を定理 2に示す。ただし、ここでは簡単のため F_k = Iとしたが、 F_k≠Iの場合も同様に求めることができる。以下に、数値的に安定な状態推定アルゴリズムを実現するための、ハイパー H_∞フィルタを示す。

(定理 2)

ただし、

R_k

Re,k (23)

であり、 Θ ( ）は J-ュニタリ行列、すなわち θ(&μθ (ん) ^Γ = «7を満たし、 - ① 、 jは単位行列である。また、は行列 ii の 1列目の列ベクトルを表す。なお、式（21)、（22)において、 —¹およびは削除可能である。

図 4に、定理 2の平方根アレイアルゴリズムの説明図を示す。この計算アルゴリズムは、図 2に示した定理 1のフローチャート中の H_∞フィルタの計算（S105)で用いることができる。

本推定アルゴリズムは、∑ | ，をリカツチ型の更新式で求める代わりに、その因数行列∑ ^^^(Σ ^の平方根行列）を jーュニタリ変換に基づく更新式で求めている。このとき生じる 1ー1ブロック行列と 2— 1ブロック行列からフィルタゲイン _kを図示のように求めている。このため、、1/2 ^1/2 K_s k k- k k-1 k k- >〇となり、の正定性は保証され、数値的に安定化できる。なお、定理 2の Η_∞フィルタの単位ステップ当たりの計算量は 0(Ν²)のままである。

なお、図 4において、」 ¹は削除可能である。

まず、処理部 101は、式（22)の左辺の行列式の各要素に含まれる項を記憶部 10 5から読み出し又は内部メモリ等から得て、 Jーュニタリ変換を実行する（S301)。処理部 101は、求めた式（22)の右辺の行列式の要素からシステムゲイン K_k、 K_s, を式（21 )に基づき計算する（S303、 S305)。処理部 101は、式（20)に基づき状態推定値 x | _kを計算する（S307)。

5.状態推定のための数値的に安定な高速アルゴリズム

上述のように、定理 2の1~1_∞フィルタの単位ステップ当たりの計算量は 0(N²)のままである。そこで、計算量の対策として、 ϋ=^₊₁Ψ, H_k=[u(k), ■■■， u(0), 0， ■ ■, 0]のとき、 =[ ， Ο^τ]^τの 1ステップ予測誤差の共分散行列 I_{k+1 |k}が +i|ん-― ^∑kk-x^^T = ~L_kR-lL L_k = (24)

0 を満たすことを利用して、 _{k+1 |k}の代わりに次元の低い] ^ (すなわち L~_k)を更新することを考える。ここで、 R_r,_k =

と表されることに注意すれば次の定理 3が得られる。

(定理 3) xk\k =

K_s,k = ¾ (：, !)/¾(!, 1) , K_k= p (K_kRj)Rl_k (62) β,

0

ここで、 e(k)は任意の j-ュニタリ行列であり、 = ₊ が成り立つ。

ただし、

9 ¹5 0 - ¹ 「1 0

R^RiJ^, Rl =

0 J Jl =

0 = -^fe|fc-l^fc|fc-l

一 1

Re,k (23)

なお、定理 3の証明は、後述する。

上式は、 K— _k( = P一 ^1/2K_k)の代わりに Κ,-についても整理することができる _c さらに、次の J—ュニタリ行列

を用いれば定理 4の高速化した状態推定アルゴリズムが得られる。ただし、 Ψはシフ卜行列を表す。

(定理 4)

if_s,_fc (_:,l)/i?_e,_fc(l,l) (26)

れ e,fc+l R C_k+1L_kR_rtkL_k C_k+1 (29)

Rr,k一 L_k C_k+1K_{e k}Ck+iL, (30) ただし、

k+l

， H_k+l = [u_k+1 u(k + l- N)j = [u(k + 1) u_k], Hi = [u(l),0,...,0] も +1 J

-^e.l = -H P o

I + ό^'ψέ^， Hj = _1|0 = diag{p², ³,..., +²}, p = l - _X(_V)

0 -frrj

1 0

Ln = 0 0 一 1 0

, ^0 = 0, x_olo = ic₀ , K_k = p' K_k (31) 0 p—^N

0 1 であり、 diagi'}は対角行列、 R_e,_k+1(1, 1)は行列 R_e,_k+1の 1—1成分をそれぞれ表す。また、上式は K—_kの代わりに K_kに関しても整理できる。

本高速アルゴリズムは、次の因数分解 - ∑^_ΧΨ^Τ = -L_kR~lL (32) における L~_keR^(N+1)x2の更新によってフィルタゲイン K_s, _kを求めているので、単位ステツプ当たりの計算量は 0 (N + 1 )で済む。ここで、次式に注意されたい。

図 5に、定理 3の数値的に安定な高速アルゴリズムのフローチャートを示す。この高速アルゴリズムは図 2の!"!_∞フィルタの計算ステップ（S1 05)に組み込まれ、 T一イタレーシヨンによって最適化される。よって、存在条件が満たされる間は r_fは除々に減少されるが、満たされなくなった時点で、図示のように： _fは増加される。

H_∞フィルタリングアルゴリズムは以下のように要約することができる。

[ステップ S401 ] 処理部 1 01は、再帰式の初期条件を図示のように定める。なお、 Lは最大データ数を示す。

[ステップ S403] 処理部 101は、時刻 kと最大データ数 Lとを比較する。処理部 1 01 は、時刻 kが最大データ数より大きければ処理を終了し、以下であれば次のステップに進む。（不要であれば条件文を取り除くことができる。または、再スタートする。 ) [ステップ S405] 処理部 1 01は、フィルタゲインに対応する項 K_{k+ 1}を式（27 )、 (3 1 )を用いて再帰的に計算する。

[ステップ S406] 処理部 1 01は、 R_e, _{k+ 1}を式（29)を用いて再帰的に計算する。

[ステップ S407] 処理部 1 01は、さらに K_s, _kを式（26)、（31 )を用いて計算する。

[ステップ S409] 処理部 1 01は、ここで、存在条件 EXC>0を判定し、存在条件を満たせばステップ S41 1に進む。

[ステップ S41 3] —方、処理部 1 01は、ステップ S409で存在条件を満たさなけれぱ r _fを増加し、ステップ S401に戻る。

[ステップ S41 1 ] 処理部 1 01は、式（25)の1~!_∞フィルタのフィルタ方程式を更新する。

[ステップ S41 5] 処理部 1 01は、 R_r. _{k+ 1}を式（30)を用いて再帰的に計算する。また、処理部 101は、 _k+1を式（28)、（31)を用いて再帰的に計算する。

[ステップ S419] 処理部 101は、時刻 kを進ませて（k = k+1)、ステップ S403に戻リ、データがある限り続ける。

なお、処理部 101は、 H_∞フィルタ計算ステップ S405~S415及び存在条件の計算ステップ S409等の各ステップで求めた適宜の中間値及び最終値を必要に応じて適宜記憶部 105に記憶し、また、記憶部 105から読み出すようにしてもよい。

6.エコーキャンセラ

つぎに、エコーキャンセリング問題の数理モデルを作成する。

まず、受信信号 iu_k}がエコーパスへの入力信号となることを考慮すれば、エコーパスの（時変）インパルス応答 {hjk]}により、ェコ一{d_k}の観測値 {y_k}は次式で表される。

N-1

Vk = d_k + v_k

+ v_k, 二 0， 1,2，··· (33) ここで、 u_k, y_kはそれぞれ時刻 t_k( = kT;Tは標本化周期)における受信信号とエコーを表し、 v_kは時刻 t_kにおける平均値 0の回線雑音とし、 hi[k], i = 0, N-1 は、時変インパルス応答であり、そのタップ数 Nは既知とする。このとき、インパルス応答の推定値 {h [k]}が時々刻々得られれば、それを用いて次のように疑似エコーが生成される。

ΛΓ-1

= ^hi[k}u_k- = 0,1,2，.·' (34)

i=0

これをエコーから差し引けば (y_k— cT_k=0)、エコーをキャンセルすることができる。ただし、 k— i<0のとき u_k— = 0とする。

以上より、問題は直接観測可能な受信信号 iu_k}とエコー iy_klからエコーパスのインパルス応答 {hi[k]}を逐次推定する問題に帰着できる。一般に、エコーキャンセラに H_∞フィルタを適用するには、まず式（32)を状態方程式と観測方程式からなる状態空間モデルで表現しなければならなし、。そこで、問題がィンパルス応答 {hi[k]}を推定することであるから、 ihi[k]}を状態変数 x_kとし、 w_k程度の変動を許容すれば、エコーパスに対して次の状態空間モデルを立てることができる。

, N

Xk+i Xk + G_kw_kl

Vk (36)

ただし、

1 IT

Xk = … , w - 1 [ w_k = (i),

H_k = , Uk-N+l) このような状態空間モデルに対するハイパーおよび高速 H_∞フィルタリングアルゴリズムは先に述べて通りである。また、インパルス応答の推定の際、送信信号の発生を検知するとその間推定を中止するのが一般的である。 7.インパルス応答に対する評価

(動作の確認）エコーパスのインパルス応答が時間的に不変であり（hjk] = h 、かつそのタップ数 Nが 48である場合について、シミュレーションを用いて、本高速アルゴリズムの動作を確認する。

47

t'=0 (38) なお、図 6は、ここでのインパルス応答 {h の値を示す図である ₍ ここで、インパルス応答 { 山₌₀ ²³は、図示の値を採用し、その他 {h山 ₌₂₄ ⁴⁷は 0とする。また、 _kは平均値 0、分散 σ_ν ² = 1.0 Χ ΐ (Γ⁶の定常なガウス白色雑音とし、標本化周期 Tを便宜上 1. 0とする。

また、受信信号 {ujは次のように 2次の ARモデルで近似する。

u_k= 一，+ a₂u_k - ₂+w_k (39) ただし、 0^=0. 7, ₂ = 0. 1とし、 w_k' は平均値 0、分散 CT_W' ² = 0, 04の定常なガウス白色雑音とする。

(インパルス応答の推定結果）

図 7に、定理 3の数値的に安定な高速アルゴリズムによるインパルス応答の推定結果を示す。ここで、図 7 (b)の縦軸は、

を表す。

これより、本高速アルゴリズムによって良好に推定出来ていることがわかる。ただし, P =1一 f)、ズ（r_f)= T_f一²、 Χ^Λ ₀|₀ = 0、∑ ！ ₀ = 201とし、計算は倍精度で行つた。また、存在条件を確認しつつ、 r_f=5.5 と設定とした _c

8. 定理の証明

8-1. 定理 2の証明

次の関係式

0

が成り立つとき、両辺の 2 X 2ブロック行列の各項を比較すれば次式が得られる。

(42)

∑k₊i\k + p-^lK_kR_e-lK^T _k = ρ-'Σ^ (43) ：れは定理 1の F_k=Iのときの式（1 3)のリカツチ方程式と一致する。ただし,

0 (44) 一方、 AJA^T= BJB^Tが成り立つとき、 Bは Jーュニタリ行列 Θ (k)を用いて Β = Α Θ (k)と表すことができる。よって、式（40)より定理 1のリカツチ方程式は次式と等価でめ W o

(45)

8

(

このとき、上式の両辺 J= ( ㊉一S)—ノル厶を比較すれば、左辺の X、 Υ、 Ζを以下のように決定することができる。ただし、 Sは対角成分が 1または一 1をとる対角行列とする。 (1, 1)一ブロック行列

― Re,k + (-Re,fc+1 - ¾+l)一 ( e,fc - Rk) = -Re.fe+l

1 = Ji {C^T _k =

(2, 1)—ブロック行列

YJ_XX τ

[

0

+1

0

+

0

：れより、 y である

(2, 2)—ブロック行列一 ZSZ^T + YJ_XY^T

T

0 T

p~^lL_kR~lL_k

Kk

0

+ p— [∑, k+l\k

A ρ~^ιΨ 一∑k\k-lj + P一¹

0

-, Τ

^Λ +1

= -^∑_k+llk^^T +∑k+2\k+l + R

0 'e, +1

0 これより、

、4 Z. T

-ZSZ^T = ¾₊₂|_ft+1― ^∑_k+Mk^^T = -L_k+1R~!₊₁SR; ₊₁Ll₊₁ となり、 Z = L_k+1R~I₊₁ を得る _c

8 -3. 定理 4の証明

観測行列 H_kがシフト特性をもち、かつ

J = (J, ® -S) のとき、定理 2と同様な方法によって次の関係式が得られる ₍

0

o

L e(k) (46) k+i I —

0 P ² L_k ただし、

Θ ( ( ん雨/?； _+lJX )

とし、 ∑(k)^T(R_eth θ -R_r>k)∑(k) = (R_eM1 Θ -R_r,k^) となるように R_r, _k+1を決定する。次に、式（46)の 3行目【こ R_r, _k+1の更新式を新たに追加すれば、最終的に次式が得られる。

,_1 iT ^T

Re,k一 C_k+iL_kn_r<kL_k C_k+l 0

' 0 · ■i n - ^Γ 0 "

K_k e, ^^ +1 も

^ j

0 ^r,k一 Lf. ^_k+1R_e lCk+l _k

(48) この両辺の 3 x 2ブロック行列の各項の対応から次のゲイン行列 K__kの更新式が得られる。

e,fc+l e,fc一

1

, ―

産業上の利用可能性

一般に、通常の民間の通信機器などでは、コストと速度の面から単精度で計算が行われる場合が多い。このため、本発明は実用的な状態推定アルゴリズムとして様々な産業分野にその効果をもたらすであろう。また、本発明は、通信システムや音響システムにおけるエコーキャンセラ、音場再生又は騒音制御などに適用することができる。

Claims

次式で表される状態空間モデルに対して、

xk+i = F_kx_k+G_kw_k

y_k= H_kx_k+v_k α ι

主Ε Π

ここで、

X_k :状態ベクトルまたは単に状態

w_{k :}システム雑音の v_{k :}観測雑音一一.

y_k :観測信号囲

:出力信号

F_k :システムのダイナミックス

G_k :駆動行列

評価基準として忘却係数Oで重み付けされた外乱からフィルタ誤差への最大エネルギーゲインを予め与えられた上限値 r_fに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数 /0の最適化を同時に行うためのシステム推定方法であって、処理部は、上限値 _f、フィルタの入力である観測信号 y_k、観測行列 H_kを含む値を記憶部又は入力部から入力するステップと、処理部は、前記上限値 r_fに従い、状態空間モデルに関連する忘却係数 Pを決定するステップと、

処理部は、記憶部から初期値又はある時刻の観測行列 H_kを含む値を読み取り、前記忘却係数 ;0を用いて次式で表されるハイパー H_∞フィルタを実行するステップ ^x k| k^=Fk-i^x k-1 | k- i +^Ks, k(yk— ^Hk^Fk - 1^X - 1 | k- 1)

ここで、

観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

K_{s k}:フィルタゲイン

処理部は、ハイパー H_∞フィルタに関する求められた値を記憶部に記憶するステップと、処理部は、求められた観測行列、又は、観測行列 Hjとフィルタゲイン K_s. iによリ、前記上限値 r_f及び前記忘却係数 Pに基づく存在条件を計算するステップと、処理部は、上限値 r_fを小さくしていき前記ハイパー H_∞フィルタを実行するステツプを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、

を含む前記システム推定方法。

2. 処理部は、前記存在条件を次式に従い計算する請求項 1に記載のシステム推定方法。

3. 処理部は、前記存在条件を次式に従い計算する請求項 1に記載のシステム推定方法。

-Q-i + Pi) > 0, = 0,···, (18) ここで、

(19)

4. 前記忘却係数 p及び前記上限値 χ _fは、次式の関係である請求項 1に記載のシステム推定方法。 <ρ=Λ-χ{γ,)≤Λ (ただし、（ r_f)は、（1 ) = 1、ズ（∞) =0を満たす r _f の単調減衰関数）

5. 前記ハイパー H_∞フィルタを実行するステップは、処理部は、前記フィルタゲイン K_s, _kを、次式により求めることを特徴とする請求項 1に記載のシステム推定方法。

T , „\-l

K_$ik = ！ ¾も— +p) (12)

-で、

P 0

/; A|fc— 1 fc， ¾ - (14)

0<ρ = 1- χ(_7/) <1, />l

(16) なお、式（16)の右辺はより一般化することもできる ₍

x_k:状態ベクトルまたは単に状態

y_k:観測信号

:出力信号

F_k:システムのダイナミックス H_k:観測行列

x _{| k}:観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

|_kの誤差の共分散行列に対応

K_{s k}:フィルタゲイン e_{f i:}フィルタ誤差

R_e,_k:補助変数

6. 前記ハイパー H_∞フィルタを実行するステップは、処理部は、初期条件に基づき、フィルタゲイン K_s. _kを前記式（12)を用いて計算するステップと、

処理部は、前記式（11)の1"1_∞フィルタのフィルタ方程式を更新するステップと、処理部は、 ∑ _|k、∑ _{+1 |k}を前記式（13)を用いて計算するステップと、処理部は、前記各ステップを、時刻 kを進ませて繰り返し実行するステップとを含む請求項 5に記載のシステム推定方法。

7. 前記ハイパー H_∞フィルタを実行するステップは、処理部は、前記フィルタゲイン K_s, _kを、ゲイン行列 K_kを用いて、次式により求めることを特徴とする請求項 1に記載のシステム推定方法。

た (20) だ

K _βし β_s6-,,k = K_k (：, 1)/¾(1, 1) K_k

(21)

(22)

(23)

でり、 Θ ( ）は j-ュニタリ行列、すなわち θ ( μθ (ん) ^T = Jを満たし、《7 = ^^)、ま単位行列である。また、 f (:,l)は行列 K_Aの 1列目の列ベクトルを表す。

なお、式（21)、 (22)において、 J ¹および J，は削除可能である。

ここで、

観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

y_k:観測信号

F_k:システムのダイナミックス

K_{s k}:フィルタゲイン

H_k:観測行列

_kの誤差の共分散行列に対応

〇（k):Jーュニタリ行列

R_{e k}:補助変数

8. 前記ハイパー H_∞フィルタを実行するステップは、

処理部は、 K_k、 1/2

k+1 I k を前記式（22)を用いて計算するステップと、処理部は、初期条件に基づき、フィルタゲイン K_s. _kを前記式（21 )を用いて計算するステップと、

処理部は、前記式（20)の!"^フィルタのフィルタ方程式を更新するステップと, 処理部は、前記各ステップを、時刻 kを進ませて繰り返し実行するステップとを含む請求項。つ 7に記載のシステム推定方法。

9. 前記ハイパー H_∞フィルタを実行するステップは、処理部は、前記フィルタゲイン K_s. _kを、ゲイン行列 K_kを用いて、次式により求めることを特徴とする請求項 1に記載のシステム推定方法。

Xk\k = Xk-i\k-i + K_s>k(y_k一 Hk&k^k^) (61)

(62)

0

ここで、 Θ ( ）は任意の J-ュニタリ行列であり、 0_fc = が成り立つ。

ただ

(23) ここで、

|_k:観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

y_k:観測信号

^Ks.k:フィルタゲイン

H_k:観測行列

Θ (k) :J—ュニタリ行列

R_e,_k:補助変数

0. 前記ハイパー H_∞フィルタを実行するステップは、

処理部は、 K - _kを前記式（63)を用いて計算するステップと、処理部は、初期条件に基づき、フィルタゲイン K - _s, _kを前記式（62)を用いて計

ん

算するステッ 11プと、

処理部は、前記式（61)の1~1_∞フィルタのフィルタ方程式を更新するステップと、処理部は、前記各ステップを、時刻 kを進ませて繰り返し実行するステップとを含む請求項 9に記載のシステム推定方法。 1. 前記ハイパー H_∞フィルタを実行するステップは、処理部は、フィルタゲイン K_s, 、ゲイン行列 K—_kを用いて、次式により求める請求項 1に記載のシステム推定方法。

Xk\k = Xk-i\k-i + K_s>k(yk - (25)

(26)

0

(28)

K_k

" T ^ T iリ ~

¾-,fc+l = ^r,k一 L_kC _k+1R_e>kCk+l k (30) ただし、

H i = ["HI u{k + l - N)] = [u(k + 1) u_k], H_x = [w(l)， 0, ... , 0]

P 0

∑ Ho = diag{ ², …， +²}, p= l - _χ( )

0 — PT'; ²

Ln一 p''K_k (31)

なお、上式は K—_kの代わりに K_kに関しても整理できる。ここで、

y_k:観測信号

F_k:システムのダイナミックス

H_k:観測行列

x _|k:観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

K_{Si k}：フィルタゲイン；ゲイン行列 K—_kから得られる。

R_e,_k、 L~_k:補助変数

12. 前記ハイパー H_∞フィルタを実行するステップは、

処理部は、予め定められた初期条件に基づき、 K - _k+1を前記式（27)を用いて再帰的に計算するステップと、処理部は、システムゲイン K_s, _kを前記式（26)を用いて計算するステップと、処理部は、存在条件を計算するステップと、

処理部は、前記存在条件を満たせば、前記式（25)の1~1_∞フィルタのフィルタ方程式を更新し、時刻 kを進ませて繰り返し各前記ステップを繰り返し実行するステップと、処理部は、前記存在条件を満たさなければ上限値 r_fを増加するステップとを含む請求項 11に記載のシステム推定方法。 13. さらに、次式により時刻 kの状態推定値 _{X |k}から出力信号の推定値 z^v _k|kを求めるようにした請求項 1に記載のシステム推定方法。

^Z k I k― k I k

14. 前記 H_∞フィルタ方程式を適用し、状態推定値 >T_k | _kを求め、

擬似エコーを次式のように推定し、求められた擬似エコーで実際のエコーを打ち消すことによりエコーキャンセラを実現する請求項 1に記載のシステム推定方法。 dk = 0,1,2,··· (34

15. 次式で表される状態空間モデルに対して、

y_k=H_kx_k+v_k

²k^=Hk^xk

ここで、

x_k:状態ベクトルまたは単に状態

w_k:システム雑音

v_k:観測雑音

y_k:観測信号

_Zk:出力信号

F_k:システムのダイナミックス

G_k:駆動行列

評価基準として忘却係数Oで重み付けされた外乱からフィルタ誤差への最大工ネルギーゲインを予め与えられた上限値 r_fに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数 ?の最適化を同時にコンピュータに実行させるためのシステム推定プログラムであつて、処理部は、上限値 r_f、フィルタの入力である観測信号 y_k、観測行列 H_kを含む値を記憶部又は入力部から入力するステップと、処理部は、前記上限値 r_fに従い、状態空間モデルに関連する忘却係数) 0を決定するステップと、

処理部は、記憶部から初期値又はある時刻の観測行列 H_Kを含む値を読み取リ、前記忘却係数 Pを用いて次式で表されるハイパー H_∞フィルタを実行するステップと、 X k I k ⁼ Fk-i^x k-1 I k- 1 +k_s, k (Yk— Hk' k-i^x k-1 I k- 1 )

ここで、

x i _k :観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

F_k：システムのダイナミックス

K_{s k} :フィルタゲイン

処理部は、ハイパー H_∞フィルタに関する求められた値を記憶部に記憶するステツプと、処理部は、求められた観測行列、又は、観測行列 Hiとフィルタゲイン K_s,；により、前記上限値 r_f及び前記忘却係数) 0に基づく存在条件を計算するステップ

処理部は、上限値 r_fを小さくしていき前記ハイパー H_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶するステップと、

をコンピュータに実行させるためのシステム推定プログラム。 1 6. 次式で表される状態空間モデルに対して、 y_k= H_kx_k+v_k

zk^{= H}k^xk

ここで、

x_k :状態ベクトルまたは単に状態 w_k:システム雑音

v_k:観測雑音

y_k:観測信号

:出力信号

F_k:システムのダイナミックス

G_k:駆動行列

評価基準として忘却係数 Pで重み付けされた外乱からフィルタ誤差への最大工ネルギーゲインを予め与えられた上限値 r_fに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数) 0の最適化を同時にコンピュータに実行させるためのシステム推定プログラムを記録したコンピュータ読み取り可能な記録媒体であって、処理部は、上限値 r _f、フィルタの入力である観測信号 y_k、観測行列 H_kを含む値を記憶部又は入力部から入力するステップと、処理部は、前記上限値 r_fに従い、状態空間モデルに関連する忘却係数 Pを決定するステップと、

処理部は、記憶部から初期値又はある時刻の観測行列 H_kを含む値を読み取リ、前記忘却係数Oを用いて次式で表されるハイパー H_∞フィルタを実行するステップと、

X k I k ⁼ '~k-1^X k- 1 I k- 1 + ， k(Vk一 ^k^k-^ k- I | k-1)

ここで、

χ^Λ _ι<:観測信号 y。〜y_kまでを用いた時刻 kの状態 x_kの推定値

F_k:システムのダイナミックス

K_{s k}:フィルタゲイン

処理部は、ハイパー H_∞フィルタに関する求められた値を記憶部に記憶するステツプと、処理部は、求められた観測行列、又は、観測行列 Hiとフィルタゲイン K_s, iにより、前記上限値及び前記忘却係数 0に基づく存在条件を計算するステップ

をコンピュータに実行させるためのシステム推定プログラムを記録したコンビユータ読み取り可能な記録媒体。 1 7. 次式で表される状態空間モデルに対して、

xk+i =F_kx_k+G_kw_k

y_k= H_kx_k+v_k

4一

ここで、

x_k:状態ベクトルまたは単に状態

w_k :システム雑音

v_{k :}観測雑音

y_k :観測信号

出力信号

F_k:システムのダイナミックス

G_{k :}駆動行列

評価基準として忘却係数Oで重み付けされた外乱からフィルタ誤差への最大工ネルギーゲインを予め与えられた上限値 r_fに対応する項より小さく抑えるように定めた、推定アルゴリズムにおいて、状態推定のロバスト化と忘却係数 Pの最適化を同時に行うためのシステム推定装置であって、推定アルゴリズムを実行する処理部と、

前記処理部により読み取り及び Z又は書き込みがなされ、状態空間モデルに関連する各観測値、設定値、推定値を記憶した記憶部と、

を備え、

前記処理部は、上限値 r _f、フィルタの入力である観測信号 y_k、観測行列 H_kを含む値を記憶部又は入力部から入力すること、前記処理部は、前記上限値 7 ^に従い、状態空間モデルに関連する忘却係数 Pを決定すること、

前記処理部は、記憶部から初期値又はある時刻の観測行列 H_kを含む値を読み取り、前記忘却係数 Pを用いて次式で表されるハイパー H_∞フィルタを実行すること、

X k | k~~ Fk—l X k - 1 | k- 1 + ^s, k (Yk一 Hk' k - 1 ^X k - 1 | k - ここで、

|_k :観測信号 y。~y_kまでを用いた時刻 kの状態 x_kの推定値

F_k_-,：システムのダイナミックス

K_{s k} :フィルタゲイン

前記処理部は、ハイパー H_∞フィルタに関する求められた値を記憶部に記憶すること、前記処理部は、求められた観測行列 Hj、又は、観測行列とフィルタゲイン K_s, ；により、前記上限値 r _f及び前記忘却係数 )0に基づく存在条件を計算すること、前記処理部は、上限値を小さくしていき前記ハイパー H_∞フィルタを実行するステップを繰り返すことで、各時刻で前記存在条件が満たされる範囲で上限値を小さく設定し、その値を記憶部に記憶すること

を備えた前記システム推定装置。