計測自動制御学会東北支部 第 297 回研究集会 (2015.10.23)
資料番号 297-5

動的システムモデルを用いた信号予測
Signal Prediction via Dynamical System Models○村松鋭一 ∗Eiichi Muramatsu∗*山形大学*Yamagata University
キーワード : 適応推定（adaptive estimation）, 離散時間信号（discrete-time signal）, 線形システム（linear system）, 予測（prediction），状態方程式（state equation）

連絡先 : 〒 992-8510 山形県米沢市城南 4–3–16 山形大学大学院理工学研究科応用生命システム工学専攻
村松鋭一，E-mail: muramatu@yz.yamagata-u.ac.jp

1.

はじめに

は線形サーボ系の理論においても目標値信号にシステム制御系において信号の未来値を予測することにより，制御系の特性を向上させることが期待できる．例えば，ビジュアルサーボ系
1, 2) において，追従対象物の動きを予測できれば，カメラの制御に予測値を活用して追従特性を高めることができる．この例に限らず，制御系における目標値信号の予測は，制御性能の向上につながる．
本稿では，離散時間信号の未来値を予測することについて考察する．過去から現在までの信号のデータは測定によって入手可能とし，それらを用いて現時点から未来の信号の予測値をオンラインで発生するシステムの設計法を提案する．

対して設けられいる．本稿で考える予測と推定の問題においては，予測対象の信号を発生するシステムは未知システムと仮定される．
未来値予測のためのシステムのオンライン推定において，本稿では適応的なシステムを導入する．過去から現在（予測を行う時刻）まで，測定データを用いて推定システムのパラメータと状態量を逐次更新する．こうして得られるパラメータと状態の推定値を利用することにより，信号の未来値を発生できる．本稿の基礎的な結果は著者による文献 3) において述べられているが，本稿ではさらに，パラメータと状態推定に用いるシステムを低次元化できることを示す．
なお本稿では伝達関数

信号の未来値を予測するには，まずその信号

G[z] = C (zI − A)−1 Bの特性を推定することが有効となるが，本稿で(1)

は信号が仮想的な動的システムから発生されていると仮定し，信号特性の推定を，動的システ
ムのパラメータと状態量の推定に置き換える．
予測対象の信号を発生する仮想的システムは線形時不変の自由システムと仮定する．この仮定で表されるシステムの入出力信号の関係を
]
[A B
u[k]
(2)
y[k] =C 0と書き表すことがある．

–1–

2.

問題設定

2.1

信号の予測

状態方程式

予測の対象となる離散時間信号を r[k] (k =0, 1, 2, · · · ) とする．これは連続時間信号 r(t) をサンプリング周期 h0 でサンプルしたものとす

xr [k + 1] = Ar xr [k](3)

r[k] = Cr xr [k](4)を用いる．ただし，xr ∈ Rn , r ∈ R であり，る．

ArCr0 0 ...
 1 0 ...
0 1 ...
= 
 .. .. . .
.
 . .
0 0 ...
[=0 0 ...


0 −a10 −a2 
0 −a3 ..
.. .
. 1 −an]0 1(5)(6)である．

(4) 式の出力 r[k] は測定可能であるが，(3) 式における状態量 xr [k]，および，(5) 式におけるパラメータ ai (i = 1, · · · , n) は未知であるとする．信号 r[k] の特性の推定は，システムにおける (5) 式におけるパラメータ ai (i = 1, · · · , n)，Fig. 1および (3) 式における状態量 xr [k] の推定に帰Prediction system着される．
本稿で考える問題を Fig.1 を用いて説明する．
信号の未来値予測を行う現在時刻を t̄ = k̄h0 と2.3

予測値を発生するシステム

する．それよりも過去の時間における信号の値

r[k] (k = k̄, k̄ − 1, k̄ − 2, · · · ) は測定によって入手可能であり，予測に用いることができる．信号の未来値 r[k] (k = k̄ + 1, k̄ + 2, k̄ + 3, · · · ) の予測対象の信号を発生するシステム (3) 式，(4) 式に対応し，時刻 k̄ に予測値を発生するシステムを

予測値を r̂[k] (k = k̄ + 1, k̄ + 2, k̄ + 3, · · · ) と表すとする．考える問題は，予測値 r̂[k] (k =k̄ + 1, k̄ + 2, k̄ + 3, · · · ) を生成する予測システ
ムの設計法を与えることである．

2.2

仮想的な信号発生システム

信号 r[k] を予測するにあたっては，その特性を推定することが有効である．本稿では信号特性の推定のため，仮想的な動的システムを導入

xp [k + 1] = Âr xp [k](7)r̂[k] = Cr xp [k](8)と定義する．ただし，
0 0
 1 0

0 1
Âr = 
 .. ..
 . .
0 0


. . . 0 −â1. . . 0 −â2 
. . . 0 −â3 .. . . ..
. .
. . . . 1 −ân(9)し，信号 r[k] がそのシステムから発生しているである．(7) 式と (8) 式より，予測値 r̂[k] (k =と考える．仮想的な信号発生システムは線形時

k̄ + 1, k̄ + 2, k̄ + 3, · · · ) は，時刻 k̄ での状態不変自由システムと仮定する．本稿ではつぎのxp [k̄] とパラメータ âi (i = 1, · · · , n) によって決定づけられる．したがって，予測システムの–2–

設計問題においては，時刻 k̄ までに xp [k̄] と âi(i = 1, · · · , n) の値をいかに求めるかという問題を考えることになる．

適応的なオンライン推定と予測

3.

前節で述べた xp [k̄] と âi (i = 1, · · · , n) の決定のためには，(5) 式の ai の推定値と，(3) 式の状態 xr の推定値が必要となる．これらを得るために適応推定機構を導入する．

3.1



パラメータの推定未知ベクトルを[]T(10)a1 a2 . . . anと定義する．この推定値として，時間とともに..
.

0 0 ...
 0 0 ...
0 0 ...
An := 
 .. .. . .
.
 . .
0 0 ...

未知パラメータ ai (i = 1, · · · , n) をまとめたθ :=
0 −10 0 


0 0 

..
.. .
. 0 0
0 00 −1 


0 0 

..
.. .
. 
0 0 ... 0 00 0 ...
 0 0 ...


0 0 ...
A1 := 
 .. .. . .
.
 . .
0 0 ...
0 0 ...
 0 0 ...
0 0 ...
A2 := 
 .. .. . .
.
 . .

000..
.

000..
.








0 −1である．このような Ai (i = 1, · · · , n) と (6) 式の Cr に対して，調整可能なベクトルを[θ̃[k] :=Ai = Li Cr , i = 1, · · · , n]Tã1 [k] ã2 [k] . . . ãn [k]とする．このベクトルを適応的に調整する仕組

を満たす Li (i = 1, · · · , n) が存在し，それらはつぎのように表される．

みを与えるため，(3) 式，(4) 式を線形パラメータ[表現に書き換えて適応システムの設計を考える．

L1 =ここでいう適応システムとは，(3)，(4) 式の[L2 =出力 r[k] に対して r̃[k] → r[k] となるような出力 r̃[k] を発生するシステムである．そのシステ
ムの構成にあたり，(3) 式の状態方程式をつぎの[Ln =ように変形する．

xr [k + 1] = (A0 +n∑ai Ai )xr [k]i=1(13)(11)(12)−1 0 0 0 . . . 00 −1 0 0 . . . 0..
.
0 0 0 0 . . . −1]T]T]T.

また，行列 L0 を A0 − L0 Cr が安定であるような行列として定義する．
つぎの捕題は，θ̃[k] を適応調整則で調整する

ただし，

0
 1


0A0 :=  ..
 .
0001..
.

... 0 0... 0 0 
... 0 0 
. . .. .. 
. . . 0 ... 1 0

ことにより，r̃[k] → r[k] を達成するシステムの構成法を述べている．
補題１ 3) 離散時間信号 r[k] の推定値をr̃[k] := θ̃T [k − 1]ξ[k] + ξ0 [k] ,–3–(14)とする．ただし，[]Tξ[k] := ξ1 [k] ξ2 [k] . . . ξn [k](15)[]
A 0 − L 0 Cr L i
ξi [k] :=r[k], i = 0, 1, · · · , n.
Cr0(16)である．(14) 式におけるパラメータの推定値ベクトル θ̃[k] を，r[k] に対する推定誤差

e[k] := r̃[k] − r[k](17)を用いた適応調整則によって更新することにより，信号の推定値 r̃[k] とパラメータ推定値ベクトル θ̃[k] は，r̃[k] → r[k], θ̃[k] → θ̄(18)を満たす．ただし，θ̄ は定数ベクトルである．
証明は文献 3) にあるので，ここではその概略Fig. 2

Adaptive system

補題における (14) 式の r̃[k]，(16) 式の ξi (i =を述べる．(4) 式と (13) 式を用いると，(12) 式0, · · · , n) と，(17) 式の e[k] を用いた適応調整の状態方程式は

xr [k + 1] = (A0 − L0 Cr )xr [k]n∑+ai Li r[k] + L0 r[k] (19)システムは Fig.2 で表される．このシステムによって，パラメータ推定値ベクトル θ̃[k] を入手することが可能になる．

i=1と書き表すことができる．これより，(3) 式, (4)式のシステムの出力として，信号 r[k] は[]n∑
A0 − L0 Cr Lir[k]r[k] =aiCr0i=1[]A0 − L0 Cr L0+r[k] (20)
Cr0という形式で書ける．(10), (15), (16) 式の定義

3.2r[k] を発生するシステムの表現

信号の予測を行うには，r[k] を発生するシステムのパラメータと状態量が必要となることは前にも述べた．ここで，次節で述べる状態推定とそれに続く信号予測の準備として，r[k] を発生するシステムを，パラメータ推定値 θ̃[k] を用いて表現しておく．

を用いれば，

補題１における (18) 式より，時間の経過とと

r[k] =n∑i=1
T

もにai ξi [k] + ξ0 [k]= θ ξ[k] + ξ0 [k]r̃[k] ≃ r[k], θ̃[k] ≃ θ̄(21)となる．この段階では，(14) 式はとなる．この形式においては，未知パラメータr[k] = θ̄T ξ[k] + ξ0 [k]

(22)

(23)

ベクトル θ が線形に現れているため，適応調整則 5) を適用することができ，補題で述べたの特となる．(3) 式，(4) 式のシステムの出力 r[k] を性が得られる．

(21) 式で書き表せることと同様に，(23) 式で書–4–0, · · · , n) を用いて求められることがつぎに補き表した r[k] は，

x̃r [k + 1] = Ār x̃r [k](24)

r[k] = Cr x̃r [k] ,(25)

0 0 ...
 1 0 ...

0 1 ...
Ār = 
 .. .. . .
.
 . .
0 0 ...

0 −ā10 −ā2 n∑
0 −ā3  = A0 +āi Ai .
.. ..
i=1.
. 
1 −ān(26)によって定義されるシステムの出力と見なすことができる．(26) 式における āi は θ̄ の要素[θ̄ =ā1 ā2 · · ·題によって示される．
補題２ 3)n∑āi ηi [k] + η0 [k](31)i=1とすると，x̂r [k] → x̃r [k](32)が成り立つ．
この補題によって，ξi [k] を発生するシステムの状態をパラメータ推定値 āi を係数とする線形結合によって，(24) 式の状態量 x̃r [k] の推定]Tānx̂r [k] :=(27)から得られるものである．

値が得られることがわかる．

3.4

信号の未来値予測

r[k] を出力するシステムの状態方程式は (3)，(4) 式でも与えられているが，ai ̸= āi ，xr [k] ̸=

前節までで信号予測システムに必要なパラメーx̃r [k] であり，(24) ∼ (27) 式は r[k] を出力する

タと状態量の入手方法が示されたので，ここでシステムの別表現となっている．その特徴は入は時刻 k̄ における信号予測システムの設定をま手可能なパラメータ推定値 āi (i = 1, · · · , N ) をとめておく．予測を行う時刻 k̄ において，用いているところである．あとは (24) 式の状態r̃[k̄] ≃ r[k̄], θ̃[k̄] ≃ θ̄ ,量 x̃r [k] の推定値を入手できれば，r[k] の未来(33)が満たされているとする．この時刻において，値を予測することが可能になる．

(9) 式の予測値発生システムのパラメータ âi (i =3.3

1, · · · , n) を状態推定âi = ãi [k̄], i = 1, · · · , n(24) 式の状態量 x̃r [k] の推定値を得る方法をこの節で述べる．ここでは，ξi (i = 0, · · · , n)を発生している (16) 式のシステムの状態量が関係する．(16) 式を状態方程式を用いて書くと{
η0 [k + 1] = (A0 − L0 Cr )η0 [k] + L0 r[k](28)
ξ0 [k] = Cr η0 [k]{η1 [k + 1] = (A0 − L0 Cr )η1 [k] + L1 r[k](29)
ξ1 [k] = Cr η1 [k]..
.
{ηn [k + 1] = (A0 − L0 Cr )ηn [k] + Ln r[k](30)
ξn [k] = Cr ηn [k](34)と設定する．ここで，ãi [k̄] (i = 1, · · · , n) は θ̃[k̄]の要素として入手できるものである．また補題

2 の (31) 式に基づき，(7) 式の予測値発生シス
テムの状態量 xp をxp [k̄] =n∑ãi [k̄]ηi [k̄] + η0 [k̄](35)i=1と設定する．こうして得られる (7)，(8) 式のシステムは，信号の予測値 r̂[k] (k = k̄ + 1, k̄ +2, k̄ + 3, · · · ) を生成する．となる．(24) 式の状態 x̃r [k] の推定値は，(28)信号 r[k] に基づくシステムの適応推定と，時刻 ¯[k] における未来値予測は Fig.3 のようにま∼ (30) 式のシステムにおける状態 ηi [k] (i =とめられる．図中の (a) は Fig.2 で表されるシ–5–質を利用することにより，各システムで同じ状態を共有させることができ，低次元化が可能となる．後の記述を簡潔にするため，

AL := A0 − L0 Cr(36)を定義しておく．
(28) ∼ (30) 式のシステムが，r[k] を入力とし，ηi [k], ξi [k] (i = 1, · · · , n) を出力していることをより明確にするため，(28) ∼ (30) 式をまとめて
  AL0···η0 [k + 1]..
 η1 [k + 1]  .
0AL
 


= .
..
.
.


..
..
 ..
.
ηn [k + 1]
0···0
L0
 L1 
+  .  r[k]
 .. 

Fig. 3

Adaptive prediction system

ステムの動作を表す．図中の (b) は時刻 k̄ における (34) 式と (35) 式による設定を表している．
図中の (c) は (7) 式と (8) 式による予測値 r̂[k](k = k̄ + 1, k̄ + 2, k̄ + 3, · · · ) の生成を表している．

4.

システムの低次元化

信号予測のためには，(28) ∼ (30) 式のシス
テムによって，η0 [k], · · · , ηn [k]，および，ξ0 [k],· · · , ξn [k] を生成しなければならない．予測対
象の信号が n 次の仮想的システムから発生してIη0 [k]
 η1 [k]   0 
  .
..
  .
.
  .
 
 ηn [k]   0
= ξ0 [k]   Cr  ξ1 [k]    0

..
 
  .
.
 ..
ξn [k]
0

Ln0I..
.
···
0Cr..
.
······..
.
..
.
0···
..
.
..
.
0



0
η0 [k]..   η [k] 1. 

..


.
0 
ηn [k]AL
0.. . 0  η0 [k]I   η1 [k] 

.
.
0 
.
..  η [k]n. 
0 Ce(37)(38)と書く．このシステムの低次元化を考える．上式より，r から ηi (i = 0, · · · , n)，あるいは rシステムを組まなければならない．オンラインから ξi (i = 0, · · · , n) への伝達関数行列は，
−1
 (zI − AL ) Li, i = 0, · · · , n (39)あるいは C (zI − A )−1 LriL

での推定と予測を考えると，システムの低次元と表される．これらはすべて (zI − AL )−1 を含化が望まれる．

み，n + 1 個のシステムの特性多項式はすべていると仮定した場合，(28) ∼ (30) 式のシステ
ムが n + 1 個あることから，それらの総次数はn(n + 1) となり，n の値が大きい場合，高次の

この節では，(28) ∼ (30) 式のシステムの低次元化，すなわち，η0 [k], · · · , ηn [k], ξ0 [k], · · · ,ξn [k] をより次数の低いシステムから発生させることを考える．(28) ∼ (30) 式のシステムの特徴は，状態 ηi [k] (i = 0, · · · , n) に掛かる行列がどのシステムにおいても A0 − L0 Cr と Cr で同じになっており，r[k] に掛かる行列が各システ
ムにおいて異なっているところである．この性

det(zI − AL ) = z n + αn z n−1 + · · · + α2 z + α1(40)となる．そこで，この係数 α1 , · · · , αn を用いて，

010 ···0 .
.. ..
..
..
 ..
.
.
.
. 
 .
..
..
Ac :=  ..
 (41).
.
0··· ··· 01  0−α1 −α2 · · · · · · −αn–6–を定義し，これと[Bc := 0 · · ·0 1

数値例5.

]T(42)

連続時間信号(π)+ 0.2r(t) = e−0.3t sin 10t +4を用いて可制御正準形式を構成することを考える．これらを用いて (39) 式で表されているシス
テムは，
−1

 Ti (zI − Ac ) Bc
あるいは
 C T (zI − A )−1 Br icc(47)をサンプリング周期 0.1 [sec] でサンプルした離散時間信号を r[k] とし，この信号の予測をシ, i = 0, · · · , n (43)ミュレーションする．
仮想的な信号発生システムを 3 次の自由シスとも表すことができる．ただし，正方行列 Ti はテム

adj(sI − AL )Li
v1n sn−1 + · · · + v12 s + v11

..
=

.

0

xr [k + 1] =  1
0[r[k] =0vnn sn−1 + · · · + vn2 s + vn1


0 −a1

0 −a2  xr [k] (48)
1 −a3](49)
0 1 xr [k]から得られる係数 vkℓ (k, ℓ = 1, · · · , n) を用いて
v11 v12 · · · v1n
 .
..
..
.. 
Ti =  ..
(44).
.
. vn1 vn2 · · · vnn生成しながらパラメータ推定を行う．その結果と定義される行列である．

を Fig.4 ∼ Fig.6 に示す．Fig.4 は，補題 1 で述とする．ai (i = 1, 2, 3) は未知パラメータであり，状態 xr も未知とする．(45) 式，(46) 式のシステムを構成し，ηi [k], ξi [k] (i = 1, · · · , n) を(43) 式の伝達関数行列を用いると，r[k] を入べた r̃[k] → r[k] の特性を表しており，破線の力，ηi [k], ξi [k] (i = 1, · · · , n) を出力とするシr̃[k] が実線の r[k] に漸近していることが確かめ

ステムの状態方程式表現は，られる．Fig.5 は適応調整則によるパラメータのz[k + 1] = Ac z[k] + Bc r[k] 
η0 [k]T0  η1 [k]   T1  .   .  .   .  .   .   ηn [k]   Tn   ξ [k]  =  CT  z[k]0  0   ξ1 [k]   CT1  .   . 
 .   . 
 .   . 
ξn [k](45)

収束を表しており，θ̃[k] の各要素 ãi (i = 1, 2, 3)が補題 1 で述べたように定数に収束していることが確認できる．
シミュレーションでの経過時間 t = 3[sec] で(46)は，Fig.4， Fig.5 に見られるように，r̃ ≃ r となり，また θ̃ が一定値になっている．この時刻を現在時刻 k̄ と見なして，信号の予測値を発生させた．r[k] の予測値 r̂[k] (k = k̄ + 1, k̄ + 2, k̄ +CTn

3, · · · ) が Fig.6 において○印でプロットされてとなる．

おり，r(t) のサンプル値 r[k] を予測できている(28) ∼ (30) 式のシステムが生成する η0 [k],ことが確認できる．

· · · , ηn [k], ξ0 [k], · · · , ξn [k] は，状態の初期値の影響が無視できる程度に時間が経過すれば，(45)1.5
1式，(46) 式のシステムでも生成することができる．(28) ∼ (30) 式のシステムの次数が n(n + 1)次であるのに対し，(45) 式，(46) 式のシステム

0.5
0−0.5
−1
0

0.5

1

1.5

Fig. 4–7–

2.5

3

3.5

4

4.5

5

5.5

6

6.5

7

7.5

8

8.5

time[sec]は n 次であり，より低次元なシステムとして

ηi [k], ξi [k] (i = 1, · · · , n) を生成できる．

2

Convergence of r̃[k] to r[k]99.510参考文献0−1−21) F. Chaumette and S. Hutchinson, “Visual
00.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

Servo Control, Part I, Basic Approaches”,

21IEEE Robotics & Automation Magazine,

000.5

1

1.5

2

2.5

3

3.5

4

4.5

55.5

6

vol. 13, 2006.

0−12) D. Navarro-Alarcon and Y. Liu, “A Dy-−200.5

11.5

2

2.5

33.544.555.5

6namic and Uncalibrated Method to Vi-time[sec]sually Servo-control Elastic Deformations
Fig. 5

Estimates of the unknown parameters

by Fully-constrained Robotic Grippers”,
Proc. of IEEE International Conferenceon Robotics & Automation pp. 4457-4462,

1.5

2014.

10.503) E. Muramatsu, “Prediction of Discrete-−0.5Time Signals via Adaptive Estimation”,−1−1.5
00.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6time[sec]Fig. 6

Proc. 10th Asian Control Conference,2015.

Prediction of reference signal4) P. Ioannou and J. Sun, Robust Adaptive
Control, Prentice Hall, 1996.

6.

おわりに5) 新中新二, 適応アルゴリズム，産業図書，

本稿では離散時間信号の予測システムの設計法とその低次元化を考察した．信号を仮想的なシステムからの出力であると見なすことにより，信号の特性の推定を，動的システムのパラメータと状態の推定に置き換え，推定によって得られるパラメータと状態を用いて信号を予測する方法を述べ，そのシステムの次数を低くする方法を与えた．

n(n + 1) 次から n 次に低次元化が可能なことから，次数 n が大きな値であっても，計算の負荷は従来よりも重くならない．今後の研究課題として，予測対象の信号が数多くの周波数成分を含むなど，より複雑な信号の場合を想定し，次数 n の値を大きく見積もる場合の予測精度の検証が考えられる．

–8–

1990.