DEIM Forum 2023 2b-1-2

モデル予測制御のための高速時系列予測アルゴリズム
藤原

廉†,††

松原 靖子†

木村輔†

櫻井 保志†† 大阪大学産業科学研究所 〒 567–0047 大阪府茨木市†† 大阪大学大学院情報科学研究科 〒 565–0871 大阪府吹田市
E-mail: †{r-fujiwr88,yasuko,tasuku,yasushi}@sanken.osaka-u.ac.jp
あらまし

本論文では，大規模制御応答時系列データストリームにおける制御量予測手法である C-C AST について述べる．C-C AST は，制御量（Controlled sequence），動作信号，操作量の三要素で構成される制御応答時系列データから，制御量の時系列パターンを捉えることで，パターン間の遷移に基づく高速な制御量予測を実現する．より具体的には，動作信号および操作量を考慮できるように動的システムを拡張し，制御応答時系列データを適応型動的システ
ムとしてモデル化することで，重要なパターンや複雑なパターンの遷移を柔軟に表現する．提案手法は，(a) 制御応答
時系列データストリームから重要な特徴を発見し，刻々と変化していく潜在的なパターンやパターン遷移を高速かつ自動的に認識し，(b) 将来的な制御量予測を実現する．さらに，提案手法は (c) データストリームの長さに依存しない．
実データを用いた実験では，提案手法が制御応答時系列データストリームの中から重要な時系列パターンを発見し，制御量予測を高精度に行うことを確認した．さらに，最新の既存手法と比較し大幅な精度向上を達成し，その計算速度はデータサイズに依存せず，高速に動作することを明らかにした．
キーワード 時系列データストリーム，制御量予測, ストリーム処理

また，一般に工業装置から得られる制御応答時系列データは，

1 まえがき宇宙ロケットの誘導制御，飛行機の姿勢制御，自動車の車両追従制御，工作機械・工業ロボットの自動制御など，時間発展する制御対象を動的システム（Dynamic System）としてモデル化するモデル予測制御（MPC：Model Predictive Control）[1, 2]
は，数多くの分野で導入され活躍している．この制御方法では，予測した制御量（Controlled sequence）に基づいて現在の操作量を調整することで，未来の制御量をシステムの目標と一致するように繰り返し最適化を行う．
しかし，MPC は制御対象に予測モデルが存在することを前提としており，モデル化されていない場合は MPC を使用できない．そのため，制御対象の内部挙動が複雑である場合や制御対象に任意に決定できる成分（動作内容や材料など）が含まれる場合，モデル化が困難もしくはモデル化にかかるコストが大きくなる．そのような場合において，制御対象のモデル化が不要である PID 制御 [3] が使用される．その簡易性から PID 制御の適用範囲は MPC よりも広い．一方で PID 制御は現在の出力をもとに操作量を決定するため，予測値を使用する場合と比べて制御の反映までに時間がかかり，結果としてその制御性能はMPC に劣る．このような制御対象において，データからの学習によって将来の制御量の予測が可能であれば，PID 制御ではなく制御性能の高い MPC が適用可能となる．そこで，近年，深層学習をはじめとした機械学習に基づくモデル予測制御の研究に期待が集まっている [4, 5]．計算機の能力が大幅に向上したことを背景に，サンプル時間が短い工業装置等のシステムの制御，さらには高度な論理計算を要する制御問題への応用が模索されている [6–9]．

装置の機差や稼働条件によって異なる時系列パターンを持つ．
また，装置に任意に設定可能な成分 (動作内容や材料等) を含む
場合，それらの時系列パターンに関する事前情報は利用できな
い．加えて，工業装置の作業内容の切り替わりに伴って動作信号や操作量で構成される外部信号が変化し，追うように制御量が急激に変化する．そのため，将来の制御量を正確に予測するためには，外部信号に基づいて時系列パターンの急激な遷移を捉える必要がある．そしてその際に事前情報を利用することはできない．したがって，現在のデータに応じて，制御応答時系
列データから時系列パターンを柔軟かつ動的に推定し，推定した時系列パターンを外部信号に基づいて切り替えることで制御量を予測する必要がある．
そこで，本論文ではモデル化が困難である工業装置を対象とした制御量予測手法である C-C AST を提案する．C-C AST は動的システムの概念を拡張し，外部信号を陽に考慮可能な適応型
動的システムを使用することで，制御応答データの時系列パターン（本研究では，“レジーム” と呼ぶ）やパターン遷移を柔軟に表現する．より具体的には，以下の問題を扱う．
制 御 量 X = {x(1), ..., x(tc )} お よ び 外 部 信 号 Z ={z(1), ..., z(tc ), z(tc + 1), ...z(tc + ls )} が与えられたとき，現
時刻 tc から ls ステップ先の制御量 x(tc + ls ) を予測する．このとき，Z は，動作信号 Sig = {sig(1), ..., sig(tc ), sig(tc +1), ...sig(tc + ls )} および操作量 U = {u(1), ..., u(tc ), u(tc +1), ...u(tc + ls )} から構成される．
1. 1 具 体 例
図 1 はある工業装置の制御応答時系列データにおける，C-C AST と既存手法の出力の様子を比較したものである．

Value

2.5
0.0
2.50

2000

4000

Time

6000

8000

10000

6000

8000

10000

8000

10000

4000

Time

(b) 動作信号 (Sig)

2.5
0.0
2.5 3600 3800 4000 4200
Time

2.5
0.0
2.5

Value

2.5
0.0
7200 2.56800 7000
Time

Value

(c) 操作量 (U )7600 7800 8000 8200
Time(d) 提案手法の予測結果

2.5
0.0
7200 2.56800 7000
Time

Value

2.5
0.0
2.5

6000

7600 7800 8000 8200
Time(e) SSDNet の予測結果

2.5
0.0
7200 2.56800 7000
Time

Value

2.5
0.0
2.5 3600 3800 4000 4200
Time

2.5
0.0
2.5

Time

Value

2.5
0.0
2.5 3600 3800 4000 4200
Time

4000

Value

2000

Value

2.5
0.0
2.5
0

2000

Value

2.50

Value

Value

Value

Value

Value

Value

(a) 制御量 (X)0.07600 7800 8000 8200
Time

2.5
0.0
2.5
2.5
0.0
2.5
2.5
0.0
2.5

9750 10000 10250
Time9750 10000 10250
Time9750 10000 10250
Time(f) OrbitMap の予測結果

図 1: 実データにおける 5 ステップ先の予測結果：(a)–(c) オリジナルデータ，(d)–(f) 各時刻における提案手法および比較手法の予測
結果．

図 1(a)–図 1(c) はそれぞれオリジナルデータの制御量（赤色

1. 2 本論文の貢献の線），動作信号（青色の線），操作量（緑色の線）を示してお本研究では，大規模制御応答時系列データストリームのためり，図 1(d)–図 1(f) は，図 1(a) の時刻 3400 以降の 4 つの異なの制御量予測手法である C-C AST を提案する．C-C AST は次のる時刻における，C-C AST および既存手法が予測した 5 ステッ
プ先の制御量を示している．
ここで，橙色の線は各手法の予測制御量を，灰色の線はオリジナルデータの制御量を示す．
図 1(d) が示すように，提案手法は新たなパターンへ瞬時に適応し，外部信号を考慮してパターン遷移を予測することで，パターンの遷移に伴う制御値の急上昇および急降下を含む，将来特長を持つ．
（ 1 ） データ内に含まれる制御量の時系列パターン（レジーム）
に関する事前知識を必要とせず，パターンの潜在的な振舞
いを時間遷移と他の入力から把握する．
（ 2 ） 把握した時系列パターン（レジーム）を用いて最適な制御
量予測を行う．

の制御量を正確に予測できている．一方，図 1(e) より，SSDNet （ 3 ） モデル推定と予測に必要な計算コストはデータサイズに依存しない．
は，学習データに存在しない未知のパターン 3 の制御量を予測できず，精度が大きく低下していることが確認できる．同様に，
図 1(f) より，OrbitMap は，パターン遷移のタイミングを誤って予測したため，予測結果に複数のスパイクが含まれることがわかる．

2 関連研究時系列データの解析および予測に関する研究は多岐にわたり，自己回帰モデル（AR ：autoregressive model），線形動的システム（LDS ：linear dynamical systems）およびカルマンフィルタ（KF ：Kalman ﬁlters） などの古典的時系列解析を活用した手法

や，近年増加している Deep Neural Network（DNN）に基づく

非線形動的システムにおけるパラメータ集合 {s0 , p, Q, A, u, V }手法など数多く提案されている．

をレジーム θbase と定義する．

古典的時系列解析を拡張した時系列予測手法として，AW-SOM [10]，TBATS [11]，PLiF [12] をはじめ，数多く提案され3. 2 動的システム
動的システムとは，現在の出力が同時刻の操作量や過去の操ている．RegimeCast [13] は大量に生成され続ける多次元セン
サデータから潜在的なパターンをリアルタイムに推定し，適応的に将来を予測し続ける．また OrbitMap [14] は時系列データの潜在状態および状態遷移をリアルタイムに推定し，長期的な予測を可能とする．しかしながら，これらの手法は外部信号を考慮した状態遷移やモデリングに対応しておらず，制御応答時

作量，および，内部状態に依存するシステムのことである．動的システムは，さまざまな分野，たとえば，宇宙工学，自動車，
化学プラント，製薬，あるいは近年では，経済分野などにおいて応用されている．また，操作量を陽に考慮したモデルを使用することで，モデル化された現象そのものを制御できる．
動的システムにおいて，現在の操作量を u(t)，現在の内部状

系列の表現には適していない．また，大規模時系列におけるパターン発見のための手法も数多く提案されている [15–19]．これらの手法は時系列上の重要なパターンを発見することが可能

態を s(t) とすると，これらは式 (3) のように表現される．この
とき過去の入力は現在の内部状態に内包されている． また，センサーなどによって現時刻の状態が観測されるとき，センサーであるが，一方で，将来予測を目的としていない．
外部入力を考慮した DNN に基づく時系列の解析に関する研究もさかんである [20–23]．たとえば，DeepAR [23] は RNN をによる雑音を d とすると，観測値 e(t) は式 (4) のように表される．まとめると，動的システムは図 2 のように表現される．

用いて予測対象となる時系列データ以外の時系列データから特徴を抽出し，それらを用いて対象となる時系列データの将来の確率分布を予測している．また，SSDNet [20] は状態空間モデルと Transformer [24] を組み合わせた時系列予測手法である．
過去の時系列データと予測先の時刻までの外部入力を使用する
ことで，外的な要因によるトレンドや季節性を利用した時系列
データの将来予測を可能としている．しかし，これらの DNNに基づく手法は，モデルの学習において膨大な計算量が必要とされ，刻々と取得されるデータからリアルタイムにモデルを更新し，最適な予測値を推定し続けることは困難である．

ds(t)= As(t) + Bu(t)dt(3)e(t) = Cs(t) + d(4)

動的システムにおいて制御量に直接影響を与えるのは操作量とされている．したがって動的システムは操作量を考慮することは可能だが，動作信号についてはその限りではない．そのため，次節では，動的システムの概念を操作量だけでなく，動作信号も考慮できるように拡張した提案手法である C-C AST について述べる．

3 提案モデル3. 3 C-C AST

本章では提案手法のための基本的な概念と提案モデルについて述べる．

本研究の目的は，制御応答データストリームが与えられたとき，その中から重要かつ潜在的な時系列パターンを発見し，現在の潜在的なパターンおよび外部信号から将来の制御量を予測

3. 1 OrbitMap におけるレジーム

することである．そこで本研究では，3. 1 節で述べたレジーム

本節では大規模時系列予測の既存手法である OrbitMap [14]を拡張し，制御応答データストリームにおけるレジームを定義において導入されているレジーム θbase について述べる．まず、する．そして外部信号に基づいてレジームの遷移を捉えることレジーム θbase は次の 2 種類の状態から構成される．

で制御量予測を実現する．このような条件下において，提案手法は次の 2 つの特長を持つ必要がある．

•

s(t)：時刻 t におけるレジーム θbase の潜在値．

•

e(t)：時刻 t における制御量の推定値．潜在値 s(t) を用いて生成される．

(P1) 動的システムで表現されるレジーム(P2) 外部信号を考慮したレジーム遷移

OrbitMap において，複数の潜在的パターンを含んでいる時系

3. 3. 1 動的システムで表現されるレジーム (P1)

列イベントストリームを次のような潜在的非線形微分方程式として表現する．

まず，本節では単一のレジームの表現法について述べる．つまり，ここではレジームの遷移は存在しないものとする．また，

3. 1 節で定義したレジーム θbase は外部からの影響を考慮してds(t)= p + Qs(t) + AS(t)dte(t) = u + V s(t)(1)

いないため，3. 2 節にて述べた動的システムに基づくレジームθ を定義する．
(2)

ここで，初期条件を s(0) = s0 ，ds(t)/dt を時刻 t の導関数とし，S(t) を s(t) の 2 次形式の行列とする．また，それぞれの潜在的パターンを式 1，式 2 のように表現し，単一の潜在的なレジーム θ は次の３種類の状態から構成される．

•

s(t)：時刻 t におけるレジーム θ の潜在値．

•

u(t)：時刻 t における操作量．

かつ s(ts ) が soutに十分に近い場合，レジーム θj への遷
i

noise

移が発生する．具体的にはレジーム制御部が sj (ts ) = v inj

System

および，操作量 u(ts ) をレジーム θj に入力し，時刻 ts 以

Sensor降の推定値をレジーム θj を使用して推定する．
図 2: 動的システムにおける入出力の様子

ここでは，レジーム θi ，レジーム θj 間の遷移は，潜在値

si (t) または外部信号 z(t) が遷移ベクトルに十分に接近し•e(t)：時刻 t における制御量の推定値．潜在値 s(t) およびた時のみ発生する（i.e.，Ds (t) = ||si (t) − souti || < ρ かつ

操作量 u(t) を用いて生成される．

Dz (t) = ||z(t) − z out
i || < ρz ）．ここで，ρ は遷移強度，ρz は制御強度とする．

ここで，レジーム制御部（Controller）を導入する．レジー

まとめると，以下のように定式化される．

ム制御部は時刻 t における外部信号 z(t) および潜在値 s(t) を

C-S WITCHER. sF B (t) を時刻 t にレジーム制御部が出力する帰入力とし，操作量 u(t) および帰還入力値 sF B (t) を出力する．

還入力値とする．このとき，レジーム θi から θj への遷移は次

レジーム制御部および帰還入力値 sF B (t) の詳細は次節にて述のように sF B (t) が切り替わることによって表現される．

べる．



si (t) (1 <= t < ts )insF B (t) = sj
(t = ts )sj (t) (ts < t)

単一のレジームおよびレジーム制御部は図 3(a) のように表される．レジームはレジーム制御部から時刻 t の操作量 u(t) および帰還入力値 sF B (t) を与えられることで，次の時刻の潜在値

s(t + 1) および推定値 e(t + 1) を生成する．つまり，式 (3) を(6)

次の式 (5) のように拡張し，式 (5) によって潜在値 s(t + 1) をout
ここで，Ds (t) = ||si (t) − souti ||，Dz (t) = ||z(t) − z i || とし

生成し，式 (4) で e(t + 1) を生成する．

たとき，ts =arg min{t|Ds (t)<=ρ,Dz (t)<= ρz }Ds (t) となる．

ここでの si (t) はレジーム θi が時刻 t において生成した潜在ds(t)= AsF B (t) + Bu(t)dt(5)

まとめると，以下の定義を得る．
［定義 1］（レジーム θ） θ を単一のレジームにおけるパラメー
タ集合とする：θ = {A, B, C, d}．

3. 3. 2 外部信号を考慮したレジームの遷移 (P2)
次に本研究では，レジームの遷移方法として C-S WITCHER

outin
値であり，また，{souti , z i , sj } はレジーム θi とレジーム θj

における遷移ベクトルである．
まとめると，提案モデルは次の要素で構成される．
［定義 2］（C-C AST の全パラメータ集合 M ） 全レジームの集合を Θ = θ1 , ...θr ，V を全遷移ベクトル集合とするとき（i.e., {sout , z out , sin } ∈ V ），提案モデルの全パラメータ
集合 M は {Θ, V } ∈ M となる．

を提案する．C-S WITCHER において，レジーム制御部は状態遷
outin
移ベクトル (souti , z i , si ) を持ち，図 3(b) で示すように，レジーム制御部が現時刻の潜在値 s(t) および外部信号 z(t) を元に，適切なレジームに帰還入力値 sF B (t) および操作量 u(t) を入力することで各時刻の推定値 e(t) を生成する．そして，入力を受けたレジームは潜在値 s(t + 1) をレジーム制御部にフィードバックする．そして，フィードバックされた潜在値 s(t + 1)
および外部入力が遷移条件を満たしたとき，図 3(c) に示すよ
うにレジーム制御部は入力経路を遷移先のレジームに切り替え

4 アルゴリズム
本章では，制御応答時系列の予測手法である C-C AST のアル
ゴリズムについて述べる．

4. 1 問 題 定 義
ここで本手法で必要な概念について定義を行う．また表 1 に
主な記号と定義を示す．

る．このように入力先を外部信号とフィードバックされた潜在［定義 3］（制御量データストリーム X ） X を制御量から構成

値 s(t + 1) を使用して切り替えることで，外部信号を考慮したされるデータストリーム X = {x(1), ..., x(tc )} とし，現レジームの動的な遷移を実現する．すなわち，C-S WITCHER で

時刻を tc とする．

は次のような規則に従ってレジーム遷移が発生する．
［定義 4］（外部信号データストリーム Z ） z(t) を時刻 t におけ•レジーム θi は各時刻 t においてレジーム制御部から帰還入力値 sF B (t) および操作量 u(t) を受け取り s(t + 1) を式(5) に従って生成し，レジーム制御部に入力する．レジー
ム遷移が発生するまでこの操作は繰り返される．

•

レジーム制御部が時刻 ts における外部信号 z(ts )，および潜在値 si (ts ) を受け取ったとき，z(ts ) が z outに十分に近い，i

る動作信号と制御量の組 z(t) = {u(t), sig(t)} とするとき，

Z を動作信号および操作信号から構成される外部信号デー
タストリーム Z = {z(1), ..., z(tc ), z(tc + 1), ...z(tc + ls )}とする．このとき，現時刻を tc とし，ls ステップ先の予測を行うものとする．

Regime θ1

Regime θ1

noise(

Regime θ(Controller))

SystemSystemController

Controller

noise
SystemSensorSensorRegime θ2

Regime θ2

noise

Sensor

System(a) 単一レジームとレジーム制御部

Sensor

System(b) 時刻 t における制御量推定

Sensor(c) レジーム遷移

図 3: C-C AST 概要：(a) 単一のレジームはレジーム制御部から入力 (u(t), sF B (t)) を受け取ることで推定値 e(t + 1) を出力する．(b)
各時刻においてレジーム制御部は１つのレジームのみに入力を与え，推定値を出力する．(c) レジーム遷移条件を満たしたとき，レジーム制御部は入力を与えるレジームを切り替える．

Algorithm 1 C-C AST

表 1: 主な記号と定義

1: Input: (a) New value x(tc ) at time points tc

記号

定義

tc

現在の時刻

X

制御量のデータストリームx(t)

時刻 t における制御量

U

操作量量のデータストリーム(b) Updated candidate C ′u(t)

時刻 t における操作量(c) Estimated variables e(tc + ls )sig

動作信号のデータストリーム(b) New value z(tc + ls ) at time points ts + ls(c) Parameter set M = {Θ, V }
out
in
out
(d) Candidate C = {θc , sout}p , z p , sc , sc

2: Output: (a) Updated parameter set M ′3: /* Estimate C and Update M */4: {M ′ , C ′ } = C-E STIMATOR(x(tc ), z(tc + ls ), M, C)

sig(t) 時刻 t における動作信号

5: /* Forecast ls step ahead estimated value */Z

外部信号のデータストリーム：{U , Sig} ∈ Zz(t)時刻 t における外部信号：z(t) = {u(t), sig(t)}s(t)

時刻 t における潜在値e(t)

時刻 t における推定値

Xc

カレントウィンドウ：X c = X[tm : tc ]tc + lc ] を長さ lc + ls のリファレンスウィンドウとする．こZ ref

リファレンスウィンドウ：Z ref = Z[tm : tc + ls ]こで，Z ref は外部信号データストリーム Z の時刻 ts からΘ

全レジームのパラメータ集合：Θ = {θ1 , ...θr }全遷移ベクトルの集合：{sout , z out , sin } ∈ V

時刻 tc + lc までの部分シーケンスを表す．

V
M

C-C AST の全パラメータ集合：{Θ, V } ∈ M

6: e(tc + ls ) = C-G ENERATOR(z(tc + ls ), M ′ , C ′ )7: return {M ′ , C ′ , e(tc + ls )}まとめると，本論文で取り組む問題を以下のように定義する．

ここで，毎時刻において新たなエントリー x(tc ) が発生し，［問題 1］ 制御量データストリーム X および外部信号データ
ストリーム Z が与えられたとき，ls ステップ先の制御量tc が増加するものとする．そこで，最新の時刻において観測さを出力し続ける．より具体的には，各時刻 tc において，れた制御量の集合のうち，現在のレジーム θc の開始時点 ts か(a) カレントウィンドウ X c およびリファレンスウィンド

ら現時刻までの部分集合をカレントウィンドウと呼び，次のよウ Z ref を用いてパラメータ集合 M を更新する．

うに定義する．

(b) 更新後のパラメータ集合 M およびリファレンスウィ［定義 5］（カレントウィンドウ X c ） X c = X[ts : tc ] を長さ

ンドウ Z ref を用いて ls ステップ先の制御量 x(tc + ls )lc のカレントウィンドウとする．ここで，X c は制御量デーを予測する．

タストリーム X の時刻 ts から時刻 tc までの部分シーケン
これらを逐次的，かつ高速に行うことが求められる．ここで，スを表す．

現時刻 tc におけるレジームを, θc ，θc の直前に存在したレジー
カレントウィンドウ X c が与えられたとき，次の目標は，パ

ムを θp とする．このとき，(a) では，現時刻のレジーム θc ，レ

ラメータ集合 Θ から最適なレジームを発見し，式 (5)，式 (4)

out
in
out
ジーム間の遷移ベクトル {soutp , z p , sc , sc } を更新していに基づき ls ステップ先の 制御量を推定することである．

る．そこで，現在のレジーム，遷移ベクトルをまとめて，候補

また，ls ステップ先の予測をするにあたり，時刻 tc + ls まで

out
in
out
モデル C = {θc , soutp , z p , sc , sc } として保持する．

の動作信号，操作量は既知とする．そこで，現時刻で既知の動作信号および操作信号のうち，現在のレジーム θc の開始時点

4. 2 提案アルゴリズム
C-C AST は次のアルゴリズムで構成される (Algorithm 1)．

ts から時刻 tc + ls までの部分集合をリファレンスウィンドウと呼び，次のように定義する．
［定義 6］（リファレンスウィンドウ Z ref ） Z ref•= Z[ts :C-E STIMATOR ：カレントウィンドウ X c ，リファレンス
ウィンドウ Z ref ，パラメータ集合 M ，および候補モデル

Algorithm 2 C-E STIMATOR

Algorithm 3 C-G ENERATOR

1: Input: (a) New value x(tc ) at time points tc

1: Input: (a) New value z(tc + ls ) at time points ts + ls(b) Parameter set M = {Θ, V }(b) New value z(tc + ls ) at time points ts + ls(c) Parameter set M = {Θ, V }

out
in
out
(c) Candidate C = {θc , sout}p , z p , sc , sc

out
in
out
(d) Candidate C = {θc , sout}p , z p , sc , sc

2: Output: (a) Updated parameter set M ′2: Output: (a) Estimated variables e(tc + ls )3: /* Initialize */(b) Updated candidate C ′4: s(tc ) = sin
c ; ts = tc ; Z ref = Z[ts : tc + ls ]3: X c = X[ts : tc ]; Z ref = Z[ts : tc + ls ]5: while ts < tc + ls do

4: /* Update current regime */6:/* Generate latent variables and estimated variables */out
5: {θc , sin}=c , sc7:(s(ts ), ..., s(tc + ls ), e(tc + ls )) = generate(θc , Z ref , s(ts ))8:/* Search transition */9:

tbest = tc + ls

10:

for i = 1 to r doarg min f (θc , sout, sincc , X c , Z ref )out ,sin
θc ,sccin
6: if f (θc , sout
,
s
,
X
,
Z ref ) > ρ then
c
c
c

7:

/* Find a better regime in Θ */8:out{θc , sin}=c , scarg minf (θc , sout , sin , X c , Z ref )θ∈Θ,sout ,sin

if f (θc , sout, sincc , X c , Z ref ) > ρ then9:10:/* Create new regime */

11:

out
{θc , sin
} = RegimeCreation(X c , Z ref )c , scΘ = Θ ∪ θc

12:
13:

end if11:/* Find any transition vector that is to close trajectory */12:

if (sout, z out, sincci ) ∈ V then

13:

14:

arg min||s(t) − sout||c||s(t)−sout||<ρ,cts <t<tc +ls
if ||z(t) − z out|| < ρz and t < tbest thenct=tbest = t; θf = θi // Select shortest transition15:

14: end if

16:

15: /* Detect segment Xc end and Insert new transition vector */17:

end if

16: if f (θc , sout, sincc , X c , Z ref ) > ρ then18:

end for

17:

out
in
V = V ∪ {soutp , z p , sc }19:/* Transition θc → θf */

18:

out
out
(sout
, z(tc ))p , z p ) = (sc20:θc = θf ; ts = tbest ; s(ts ) = sin
f

19:

out
θc = sin=∅c = sc21: end while

20: end if

end if

22: return e(tc + ls )21: M ′ = {Θ, V }
out
in
out
22: C ′ = {θc , sout}p , z p , sc , sc

23: return {M ′ , C ′ }Regime Creation. カレントウィンドウ Xc に未知の時系列パターンが現れたとき，新たなレジーム θ = {A, B, C, d} を推定

C を与えたとき，パラメータ集合 M ，および候補モデルする必要がある．OrbitMap [14] においてパラメータ推定アルゴ

C を更新する (Algorithm 2)．

リズムに expectation–maximization（EM）アルゴリズムを採用しているが，このアルゴリズムにおいて外部入力は白色ノイズ•

C-G ENERATOR ：候補モデル C ，およびパラメータ集合 M ，として扱われるため，今回の提案手法には適さない．そこで，本

リファレンスウィンドウ Z ref を用いて，ls ステップ先の研究では Adjoint Sensitivity Method [25,26] を拡張したパラメー

推定値 e(tc + ls ) 予測する (Algorithm 3)．

タ推定アルゴリズムを提案する．具体的には，Adjoint Sensitivity

4. 2. 1 C-E STIMATOR

Method を使用し，損失関数に対するパラメータ {A, B} の勾カレントウィンドウ X c ，リファレンスウィンドウ Z ref ，パ配を求め，それを用いた勾配法で {A, B} を推定し，パラメー

ラメータ集合 M ，候補モデル C を与えたときを考える．こタ {C, d} を Levenberg-Marquardt（LM）アルゴリズム [27] をin
こで，損失関数を f (θc , soutc , sc , X c , Z ref ) と定義する．f (·)

用いて推定する．これを誤差が収束するまで繰り返す．

は推定値とオリジナルデータのフィッテング精度を示す（i.e.,∑ tc
inf (θc , soutc , sc , X c , Z ref ) =t=ts ||x(t) − e(t)||）．
関数 C-E STIMATOR は次の手順で行われる．
in
out
in
(I) (θc , soutc , sc ) を f (θc , sc , sc , X c , Z ref ) が最小になる

ように更新する．
in(II) f (θc , soutc , sc , X c , Z ref ) > ρ となるとき，パラメータ集
in合 M から，f (θc , soutc , sc , X c , Z ref ) が最小になるパラメータを探索する．
in(III) f (θc , soutc , sc , X c , Z ref ) > ρ となるとき，新たなレジー

ムを生成し，パラメータ集合 Θ に追加する．

(IV) 遷移ベクトルの生成（もし存在すれば）：これまでの手順のなかで，遷移ベクトル集合 V に存在しない遷移が確認され4. 2. 2 C-G ENERATOR
次にモデル候補 C ，パラメータ集合 M およびリファレン
ンスウィンドウ Zref が与えられたときを考える．関数 C--

G ENERATOR は次の手順で行われる．
(I) s(tc ) を初期値として，時刻 tc + ls までの推定値をレジー
ム θc を使用して生成する．

(II) 生成された潜在値 {s(tc ), ..., s(tc + ls )} および，外部信号{z(tc ), ...z(tc +ls )} について，Ds (ts ) = ||si (ts )−souti || <
ρ かつ Dz (ts ) = ||z(ts ) − z out
c || < ρz となる時刻 ts が存在すれば，潜在値を {s(tc ), ..., s(ts − 1), sin
f , ..., sf (tc + ls )}inと更新する．このとき，{sinf , ..., sf (tc + ls )} は，sf を初期値として，時刻 tc + ls までの推定値をレジーム θf を使用して生成する．

たとき，新たな遷移ベクトルを生成し，遷移ベクトル集合

V に追加する．

(III) もし ts < tc + ls ならば，tc = ts として手順 (II) に戻る．

そうでないなら，推定値を出力して終了する．

パターン遷移を予測できていないため，制御量の急激な変動を予測できていない．

5 評価実験
本論文では C-C AST の有効性を検証するため，実データを用いた実験を行った．本章では以下の項目について検証する．

Q1 制御量予測に対する提案手法の有効性

5. 2 Q2：提案手法の精度
本節では，C-C AST の予測精度を検証するため，既存手法である SSDNet および OrbitMap と精度比較を行った．
図 5 は，工業制御応答データストリームにおける C-C AST のQ2 リアルタイム予測に対する提案手法の精度の検証

予測結果の精度を示している．具体的には，オリジナルデータ

Q3 データストリームに対する計算時間の検証の制御量と 5 ステップ先の予測制御量の推定値の二乗平均誤差

また，実験は 512GB のメモリ，AMD EPYC 7502 2.5GHz の（RMSE：root mean square error）および平均絶対誤差（MAE：

32 コア CPU および NVIDIA RTX A6000 48GB GPU を搭載しmean absolute error）を示している．図に示す通り，提案手法はた Linux マシン上で実施した．

最新の既存手法である SSDNet および OrbitMap と比較し，2 つの評価指標のどちらにおいても高い予測精度を持つ．SSDNet

比較手法. 有効性を検証するために，以下の最新の既存手法とはモデルをオンラインに更新できず，テスト区間に含まれる未比較を行った．

知のパターンに対応できないため精度が低下している．また，OrbitMap は外部信号を考慮した時系列パターンの表現やパター•

OrbitMap [14]：大規模時系列イベントストリームのリアル
タイム予測手法．モデル推定時の閾値は学習時に最もフィッ
テング性能が良い値を使用した．

•

ン遷移の予測ができないため，適切に制御量を予測することができない．
また，状態の遷移がさらに複雑になった場合，初めて現れる

SSDNet [20]：Transformer に基づく深層モデルを用いた時遷移に関しては提案手法を含め全手法で予測精度が低下する．

系列予測手法であり，将来の外部信号を考慮して予測を行

しかし，その後に再び同じ状態遷移が現れた際，提案手法はそう．Transformer のエンコーダ層，デコーダー層のユニット

れに適応可能である．一方で，SSDNet はモデルをオンライン

数はそれぞれ 3 とした．また，最適化アルゴリズムとしてに更新できないため，初めて現れた時と同様に予測精度が低下

Adam [28] を使用し，300 エポック学習した中で最もフィッ

する．そして OrbitMap は状態遷移を適切に予測できないため，テング性能が良いモデルを使用した．

初めて現れた時と同様に予測精度が低下する．つまり，状態遷移が複雑になるほど各手法の予測精度は低下するが，提案手法

データセット. ある工業装置の制御応答時系列データストリー

は比較手法と比較して精度の低下が小さいと考えられる．

ムを使用した．データセットは平均値と分散値で正規化 (z-normalization) して使用した．また，データストリームの時刻 0

5. 3 Q3：提案手法の計算時間

から時刻 3400 までを学習データとして使用した．具体的には，続いて，提案手法の計算コストについて検証する．図 6(a) は，

SSDNet では学習データをモデル学習のみに使用し，提案手法

各時刻 tc における計算コストを，提案手法の C-C AST と既存手

および OrbitMap は初期パラメータの推定およびモデル学習に法である SSDNet および OrbitMap と比較したものである．ま用いた．また，SSDNet のウィンドウサイズについて，10 からた，図 6(b) では，提案手法のデータストリーム全体の計算時間

50 まで 10 刻みでそれぞれ実験を行い，学習データにおいて誤の分布を示している．モデルの実行において，提案手法および差が最小となったウィンドウサイズを用いた．

OrbitMap は CPU のみを用い，深層学習である SSDNet は CPU

5. 1 Q1：提案手法の有効性
本節では，制御応答時系列データに対する C-C AST の予測能力を検証する．図 1 および図 4 は実際の制御応答時系列データに対する提案手法および比較手法の予測結果である．そして，
図 4 では各パターン（レジーム A-C）の開始時における予測結果を比較する．すでに 1 章の図 1 においても示したように，提案手法は外部信号を利用することで動作信号に基づいてパター
ン遷移を予測可能であるため，遷移の開始時に発生する制御量の急激な変動を予測できている．一方，図 4(b) および図 4(c)より，比較手法は制御量の急激な変動に対応できていないことが分かる．特に，学習データに含まれないレジーム C の変化において，制御量の立ち上がりが遅れていることを確認できる．

SSDNet は，C-C AST と同様に外部信号を考慮できるため，部
分的に急激な変動に対応できているが提案手法ほど精度は高くない．一方，OrbitMap は，外部信号を扱えず動作信号に基づく

および GPU を併用した．なお，図中の y 軸は対数スケールで示している．
図 6(a) に示す通り，大規模時系列の高速予測手法である Or-

bitMap と同様に, C-C AST はデータストリームの長さに依存せず，高速に動作することが確認できる．さらに，提案手法は，深層学習である SSDNet と比較して大幅な性能向上を達成した．
具体的には，SSDNet と比較し 136,700 倍の高速化を実現している．
また，本実験で使用した装置において MCMPC [29] を適用する場合，予測値を求めた後に操作量を最適化する必要がある．そのため，データサイズの大きさに対する計算時間だけでなく，計算時間の分布も非常に重要である．図 6(b) は，デー
タセットの各時刻の処理にかかる時間の度数分布をプロットしたものである．本実験で使用したデータセットのサンプリング
レートは 100ms であり，これはモデル予測制御のためには少

5200

5200

5190
Time

7720

RegimeB

2.5
0.0
2.5

5200

7730
Time

7740

RegimeC

2.5
0.07020 2.5

7010

7015
Time

7720

(b) SSDNet の予測結果

RegimeA

2.5
0.0
2.5

7015
Time

Value

2.5
0.0
2.5

Value
5190
Time

7010

(a) 提案手法の予測結果

RegimeA

2.5
0.0
2.5

RegimeC

2.5
0.07020 2.5

RegimeB

7730
Time

7740

RegimeC

2.5
0.07020 2.5
Value

5190
Time

RegimeB

Value

Value

2.5
0.0
2.5

Value

Value
Value
Value

RegimeA

2.5
0.0
2.5

7010

7015
Time

7720

7730
Time

7740

(c) OrbitMap の予測結果

図 4: 各パターン (レジーム) 開始時における予測精度：(a) 提案手法および (b)–(c) 比較手法の予測結果

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

0

6 む す び
本論文では，制御応答時系列データストリームにおける制御量予測手法である C-C AST について述べた．C-C AST は，外部信号である動作信号および操作量を考慮できるように動的システムを拡張し，制御応答時系列データを適応型動的システムとしてモデル化することで，制御量の重要なパターンや外部信号

0
SSDNet OrbitMap C-Cast

SSDNet OrbitMap C-Cast(a) RMSE

によるパターンの遷移を柔軟に表現し，最適な制御量予測を実現する．

(b) MAE

図 5: 工業制御応答ストリームにおける制御量の予測精度

実データを用いた実験では，C-C AST が制御応答時系列デー

C-CastOrbitMap
SSDNet

102
100
10 2
7000

8000
9000
Sequence length(a) 各時刻における計算速度

10000

Number of data points(in thousands)

Wall clock time (sec.)

タに対し，制御量の時系列パターンやパターン遷移を高速かつ
C-Cast3K

今後の課題として，様々な制御応答時系列データをより柔軟に表現するための高度なモデル学習や，予測モデルを利用し，最2K

適な操作量を決定する手法について検討していく予定である．

1K
0

謝辞

本 研 究 の 一 部 は JSPS 科 研 費 ，JP20H00585，

JP21H03446，JP22K17896, 国 立 研 究 開 発 法 人 情 報 通 信 研 究<50 50 - 100 100<
Running times (msec.)(b) 平均計算速度

図 6: 各時刻 tc における計算コスト（a）と平均値（b）

なくとも 100ms 以内に計算を完了させる必要があることを意味する．3.7K データポイント（全データポイントの約 97%）
では 100ms 以内に，3.1K データポイント（全データポイン
トの約 83%）では 50ms 以内に計算が完了する．このように、

C-C AST は制御部分に操作量を最適化するための十分な時間を与えることができる。

継続的に発見し，制御量予測を高精度に行うことを確認した．

機構委託研究 NICT 03501, 総務省 SCOPE JP192107004，JS-

TAIP 加速課題 JPMJCR21U4，ERCA 環境研究総合推進費 JPMEERF20201R02, の助成を受けたものです．
文

献[1] Eduardo F Camacho and Carlos Bordons Alba. Model predictive
control. Springer science & business media, 2013.
[2] David Q. Mayne, James B. Rawlings, Christopher V. Rao, and Pierre
O. M. Scokaert. Constrained model predictive control: Stability andoptimality. Autom., 36(6):789–814, 2000.
[3] Karl Johan Åström and Tore Hägglund. PID Controllers: Theory,
Design, and Tuning. ISA - The Instrumentation, Systems and Automation Society, 1995.
[4] Lukas Hewing, Kim P Wabersich, Marcel Menner, and Melanie N

[5]

[6]

[7]

[8]

[9]
[10]

[11]

[12]

[13]

[14]

[15]
[16]

[17]

[18]

[19]

Zeilinger. Learning-based model predictive control: Toward safe
learning in control. Annual Review of Control, Robotics, and Autonomous Systems, 3:269–296, 2020.
Benjamin Karg and Sergio Lucia. Efﬁcient representation and approximation of model predictive control laws via deep learning. IEEE
Trans. Cybern., 50(9):3866–3878, 2020.
D.Q. Mayne and H. Michalska. Receding horizon control of nonlinear systems. IEEE Transactions on Automatic Control, 35(7):814–
824, 1990.
Toshiyuki Ohtsuka. A continuation/gmres method for fast computation of nonlinear receding horizon control. Autom., 40(4):563–574,
2004.
Grady Williams, Paul Drews, Brian Goldfain, James M. Rehg, andEvangelos A. Theodorou. Aggressive driving with model predictive path integral control. In Danica Kragic, Antonio Bicchi, andAlessandro De Luca, editors, 2016 IEEE International Conferenceon Robotics and Automation, ICRA 2016, Stockholm, Sweden, May
16-21, 2016, pages 1433–1440. IEEE, 2016.
Peter Ortner and Luigi del Re. Predictive control of a diesel engine
air path. IEEE Trans. Control. Syst. Technol., 15(3):449–456, 2007.
Spiros Papadimitriou, Anthony Brockwell, and Christos Faloutsos.
Adaptive, hands-off stream mining. In Johann Christoph Freytag,
Peter C. Lockemann, Serge Abiteboul, Michael J. Carey, Patricia G.
Selinger, and Andreas Heuer, editors, Proceedings of 29th International Conference on Very Large Data Bases, VLDB 2003, Berlin,
Germany, September 9-12, 2003, pages 560–571. Morgan Kaufmann, 2003.
Alysha M. De Livera, Rob J. Hyndman, and Ralph D. Snyder. Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American Statistical Association,
106(496):1513–1527, 2011.
Lei Li, B. Aditya Prakash, and Christos Faloutsos. Parsimonious linear ﬁngerprinting for time series. Proc. VLDB Endow., 3(1):385–396,
2010.
Yasuko Matsubara and Yasushi Sakurai. Regime shifts in streams:
Real-time forecasting of co-evolving time sequences. In Balaji Krishnapuram, Mohak Shah, Alexander J. Smola, Charu C. Aggarwal,
Dou Shen, and Rajeev Rastogi, editors, Proceedings of the 22nd
ACM SIGKDD International Conference on Knowledge Discovery
and Data Mining, San Francisco, CA, USA, August 13-17, 2016,
pages 1045–1054. ACM, 2016.
Yasuko Matsubara and Yasushi Sakurai. Dynamic modeling andforecasting of time-evolving data streams. In Ankur Teredesai,
Vipin Kumar, Ying Li, Rómer Rosales, Evimaria Terzi, and George
Karypis, editors, Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD
2019, Anchorage, AK, USA, August 4-8, 2019, pages 458–468. ACM,
2019.
Varun Chandola, Arindam Banerjee, and Vipin Kumar. Anomaly detection: A survey. ACM Comput. Surv., 41(3):15:1–15:58, 2009.
David Hallac, Sagar Vare, Stephen P. Boyd, and Jure Leskovec.
Toeplitz inverse covariance-based clustering of multivariate time series data. In Proceedings of the 23rd ACM SIGKDD International
Conference on Knowledge Discovery and Data Mining, Halifax, NS,
Canada, August 13 - 17, 2017, pages 215–223. ACM, 2017.
Takato Honda, Yasuko Matsubara, Ryo Neyama, Mutsumi Abe, andYasushi Sakurai. Multi-aspect mining of complex sensor sequences.
In Jianyong Wang, Kyuseok Shim, and Xindong Wu, editors, 2019
IEEE International Conference on Data Mining, ICDM 2019, Beijing, China, November 8-11, 2019, pages 299–308. IEEE, 2019.
Kouki Kawabata, Yasuko Matsubara, and Yasushi Sakurai. Automatic sequential pattern mining in data streams. In Wenwu Zhu,
Dacheng Tao, Xueqi Cheng, Peng Cui, Elke A. Rundensteiner, David
Carmel, Qi He, and Jeffrey Xu Yu, editors, Proceedings of the 28th
ACM International Conference on Information and Knowledge Management, CIKM 2019, Beijing, China, November 3-7, 2019, pages
1733–1742. ACM, 2019.
Yasuko Matsubara, Yasushi Sakurai, and Christos Faloutsos. Au-

[20]

[21]

[22]

[23]

[24]

[25]
[26]

[27]
[28]

[29]

toplait: automatic mining of co-evolving time sequences. In Curtis E. Dyreson, Feifei Li, and M. Tamer Özsu, editors, International
Conference on Management of Data, SIGMOD 2014, Snowbird, UT,
USA, June 22-27, 2014, pages 193–204. ACM, 2014.
Yang Lin, Irena Koprinska, and Mashud Rana. Ssdnet: State spacedecomposition neural network for time series forecasting. In James
Bailey, Pauli Miettinen, Yun Sing Koh, Dacheng Tao, and Xindong
Wu, editors, IEEE International Conference on Data Mining, ICDM
2021, Auckland, New Zealand, December 7-10, 2021, pages 370–
378. IEEE, 2021.
Yuan Xue, Denny Zhou, Nan Du, Andrew M. Dai, Zhen Xu, Kun
Zhang, and Claire Cui. Deep state-space generative model for correlated time-to-event predictions. In Rajesh Gupta, Yan Liu, Jiliang Tang, and B. Aditya Prakash, editors, KDD ’20: The 26th
ACM SIGKDD Conference on Knowledge Discovery and Data Mining, Virtual Event, CA, USA, August 23-27, 2020, pages 1552–1562.
ACM, 2020.
Yuyang Wang, Alex Smola, Danielle C. Maddix, Jan Gasthaus, Dean
Foster, and Tim Januschowski. Deep factors for forecasting. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedingsof the 36th International Conference on Machine Learning, ICML
2019, 9-15 June 2019, Long Beach, California, USA, volume 97
of Proceedings of Machine Learning Research, pages 6607–6617.
PMLR, 2019.
David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim
Januschowski. Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting,36(3):1181–1191, 2020.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion
Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Isabelle Guyon, Ulrike von Luxburg, Samy
Bengio, Hanna M. Wallach, Rob Fergus, S. V. N. Vishwanathan, andRoman Garnett, editors, Advances in Neural Information Processing
Systems 30: Annual Conference on Neural Information Processing
Systems 2017, December 4-9, 2017, Long Beach, CA, USA, pages
5998–6008, 2017.
Lev Semenovich Pontryagin. Mathematical theory of optimal processes. 1987.
Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Samy Bengio, Hanna M. Wallach, Hugo Larochelle, Kristen Grauman, NicolòCesa-Bianchi, and Roman Garnett, editors, Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8,
2018, Montréal, Canada, pages 6572–6583, 2018.
Jorge J. Moré. The levenberg-marquardt algorithm: Implementationand theory. In Numerical Analysis, pages 105–116, 1978.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In Yoshua Bengio and Yann LeCun, editors, 3rdInternational Conference on Learning Representations, ICLR 2015,
San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings,
2015.
Nikolas Kantas, JM Maciejowski, and A Lecchini-Visintini. Sequential monte carlo for model predictive control. In Nonlinear model
predictive control, pages 263–273. Springer, 2009.