統計数理（2023）
第 71 巻 第 1 号 47–64©2023 統計数理研究所

特集「諸科学における統計数理モデリングの拡がり」［研究ノート］

核融合プラズマ制御に向けたデータ同化システムの開発
森下 侑哉1 ・村上 定義1 ・横山 雅之2,5 ・上野 玄太3,4,6（受付 2022 年 6 月 30 日；改訂 10 月 23 日；採択 10 月 28 日）

要

旨

核融合プラズマの高精度な解析・制御を実現するため，核融合プラズマの統合シミュレー
ションコードをシステムモデルとするデータ同化システム ASTI の開発を進めている．ASTIの最終的な目的は核融合プラズマの制御にあるが，既存のデータ同化の枠組みは制御の要素を含んでいない．そのため本研究では，逐次ベイズフィルタを拡張し，観測データを用いたシス
テムモデルの最適化と目標状態を実現する制御入力推定とを統合したデータ同化フレームワー
クを新たに開発した．本稿では，ASTI の概要を紹介するとともに，開発した制御用データ同化フレームワークについて説明する．また，数値空間上の仮想プラズマを ASTI によって制御する数値実験の結果を紹介する．
キーワード：データ同化，モデル予測制御，核融合プラズマ．
1.

はじめに

次世代の発電方法として，核融合による発電が注目されている．核融合発電の方式の一つである磁場閉じ込め核融合では，磁力線により編み込まれたカゴにプラズマを閉じ込め，1 億度を超える高温かつ高密度な状態に維持する．こうした極限状態にある核融合プラズマは，様々な時空間スケールの物理現象が混在する複雑系であり，その全体挙動のシミュレーションには，各物理現象を表現するモジュールを結合した統合コードが用いられる．核融合プラズマ全体の予測や解析，運転シナリオ策定，炉設計などに利用するために世界各地で様々な統合コードの開発が進められている．また，将来の核融合炉の運転には，プラズマを常に監視し，加熱等を制御するシステムが必要となる．この制御には，シミュレーションによる挙動予測に基づいて制御を決定するモデル予測制御（足立, 2014）が必要になると考えられ，その予測モデルとしても統合コードが重要となる．
統合コードは様々な面で活用が期待されるが，計算結果の不確実性や計算コストといった問題も抱えている．統合コードを構成する各モジュールは少なからず不確実な部分を持ってお1

京都大学大学院 工学研究科：〒 615–8540 京都府京都市西京区京都大学桂2核融合科学研究所 六ヶ所研究センター：〒 039–3212 青森県上北郡六ヶ所村大字尾駮字表舘 2–166

3

統計数理研究所：〒 190–8562 東京都立川市緑町 10–3

4

データサイエンス共同利用基盤施設：〒 190–0014 東京都立川市緑町 10–3

5

総合研究大学院大学 物理科学研究科核融合科学専攻：〒 509–5292 岐阜県土岐市下石町 322–6

6

総合研究大学院大学 複合科学研究科統計科学専攻：〒 190–8562 東京都立川市緑町 10–3

48

統計数理

第 71 巻 第 1 号 2023

り，これらの結合による相互作用も加わって，シミュレーション結果が大きな不確実性を持つことがある．特に，乱流輸送はプラズマの挙動を大きく左右する重要な物理過程であるにも関わらず，精度の高いモデルは未だ開発できていない．また，多くのモジュールで計算コストの大きい複雑なモデルを解いていることや，考慮が不十分な要素（例えば，プラズマを囲う壁の状態など）がプラズマ性能に無視できない影響を与えることを考えると，現状の統合コードを用いて実時間でのモデル予測制御を行うのは困難であると考えざるを得ない．
こうした問題を解決し，核融合プラズマの複雑な挙動を予測・制御するシステムを構築する
ため，著者らは統合コードを用いたデータ同化システム ASTI（Assimilation System for Toroidal
plasma Integrated simulation）（Morishita et al., 2022）を開発している．データ同化（ 口 他,とは，観測データを用いてシミュレーションを最適化する統計的2011; Gettelman et al., 2022）
手法であり，気象予測や海洋解析をはじめとする多くの分野で活躍している．状態変数として定義した物理変数やモデルパラメータを観測情報に基づいて逐次的に最適化することで，シス
テムモデルの挙動を実システムに近づけることができ，予測精度を高めることができる．また，解析的な用途として，観測データに対して時空間的に整合性の高いモデルパラメータや非観測量の推定を行うことができる．図 1 は，データ同化システム ASTI のイメージ図である．
ASTI によって，統合コードによるシミュレーションと観測とが確率分布を介して繋がり，観測情報を用いた統合コードの最適化やモデル推定が可能となっている．ASTI の大きな特徴は，核融合プラズマの複雑な挙動の予測制御をデータ同化を応用して実現しようとしている点にある．そのため，データ同化システムから実システム（核融合炉）への働きかけを含んだ新たなデータ同化フレームワーク DACS（Data Assimilation and Control System）を開発している．また，ASTI は解析システムとしての側面も持っており，核融合分野に蓄積されてきた多種多様な計測データを用いて，より高性能なモデル開発を行うことができる．ASTI は，数値空間と実空間とをつなぐ総合的なシステムになることが期待される．
本稿では，核融合分野のデータ同化の取り組みについて，特に ASTI の制御システムとしての側面に重点を置いて紹介する．一般的なデータ同化は，限られた観測情報から状態変数を推定するオブザーバーの技術であり，制御推定といった実システムへの能動的な働きかけを含んでいない．そこで本研究では，ASTI の基盤となるデータ同化のフレームワーク（逐次ベイズ
フィルタ）を，図 2 に示すように，制御を含む形へと拡張した．この DACS フレームワークにより，観測を用いてシステムモデルの不確実性を抑えつつ，目標状態を実現する制御入力を推定することができる．制御の分野では適応モデル予測制御に該当し，また，デジタルツイン（Liu et al., 2021; Kapteyn et al., 2021）を実現する一手法となる．これまで解析システムとし

図 1．核融合分野におけるデータ同化システム ASTI．

核融合プラズマ制御に向けたデータ同化システムの開発

49

図 2．ASTI の拡張．

て開発してきた ASTI に，DACS フレームワークに基づいた制御システムとしての機能を追加し，数値空間上の仮想プラズマを制御する数値実験を通してその有効性を検証した．次節において逐次ベイズフィルタと ASTI の概要を説明した後，第 3 節において制御への取り組みを紹介する．
2.

データ同化システム ASTI

ASTI では，データ同化手法の一つである逐次ベイズフィルタの枠組みに核融合プラズマのを組み込んでいる．ここで，簡統合シミュレーションコード TASK3D（Murakami et al., 2015）
単に逐次ベイズフィルタについて説明する．まず以下のようなシステムモデル（2.1）と観測モ
デル（2.2）から成る状態空間モデルを考える．
(2.1)

xt = ft (xt−1 , vt )(2.2)

yt = ht (xt , wt )ここで，xt は時刻 t における状態ベクトルであり，yt は観測ベクトルである．システムモデ（システムノイズ）vt
ル（2.1）は状態ベクトルの時間発展を表現しており，時間発展に伴うノイズが考慮される．式（2.1）
中の時間発展演算子 ft が統合コード TASK3D に対応している．観測モ
デル（2.2）は，観測ベクトルと状態ベクトルの関係性を表現しており，観測ノイズ wt が考慮される．観測機器に起因する測定誤差に加え，シミュレーションモデルが考慮していない（表現
し切れなかった）成分による誤差もこの観測ノイズに含まれる. システムモデルが状態ベクトル
を観測の時刻まで運び，観測モデルが状態ベクトルと観測の比較を可能にする．
逐次ベイズフィルタでは，状態ベクトルの確率分布（状態分布）の時間発展を追う．今，時刻 ti の状態ベクトルの確率分布 p(xi |y0:i ) が所与であるとする．ここで，時刻 ti までの観測情
報 y0:i ≡ {y0 , y1 , y2 , . . . , yi } はこの状態分布に反映されているものとする．システムモデルを用いることで次の観測時点 ti+1 における予測分布 p(xi+1 |y0:i ) を得ることができる．観測時点
ti+1 において観測 yi+1 が得られると，p(xi+1 |y0:i ) を事前分布とし，観測モデルを用いて計算される尤度分布を用いることでベイズの定理が適用できる．こうして状態ベクトルの確率分布
は，新たな観測情報 yi+1 を取り込んだ事後分布（フィルタ分布）p(xi+1 |y0:i+1 ) へと更新される.
この最適化操作はフィルタと呼ばれ，この事後分布を用いてさらに次時点への予測を行う. 逐
次ベイズフィルタは，この予測とフィルタの繰り返しであり，各ステップで状態分布は以下のように変化する．
予測：(2.3)

p(xi |y0:i ) → p(xi+1 |y0:i ).

50

統計数理

第 71 巻 第 1 号 2023

フィルタ：
(2.4)

p(xi+1 |y0:i ) → p(xi+1 |y0:i+1 ).

ここでは，具体的な数式の計算に立ち入らずに逐次ベイズフィルタを説明した．詳細は，（ 口
他, 2011）などを参照されたい．
一般にデータ同化を必要とするのは複雑で非線形なシステムであるため，逐次ベイズフィル
タを実現するためには工夫が必要である．非線形なシステムでも実行できる強力なデータ同（Evensen, 2003）や粒子フィルタ（PF）化手法として，アンサンブルカルマンフィルタ（EnKF）（Kitagawa, 1996）などの状態分布を有限個のアンサンブルメンバーで近似する手法がある．
個々のアンサンブルメンバーをシミュレーションモデルを用いて時間発展させることで予測分布を近似するアンサンブルを計算し，観測データに基づいてアンサンブルメンバーを再構成することでフィルタを実現する．EnKF は，状態分布をガウス分布だと仮定して（すなわち，アンフィルタを行う手法であり，行列計算によりフィサンブルの 2 次のモーメントまでを用いて）ルタを実行できる．その反面，システムの非線形性が強く，変数間の関係を線形の範囲で捉えきれない場合（分布がガウス分布から大きく離れる場合）
や観測モデルが非線形な場合には，必ずしも推定がうまくいくとは限らない. PF は分布形状に一切の仮定を置かない反面，分布を有限個のアンサンブルメンバーで近似することの影響を受けやすいため，EnKF に比べ多くのアンサンブルメンバーを必要とする（上野, 2019）．ASTI には EnKF と PF を実装しているが，多くの場合，EnKF で十分対応できることを確認している．
ASTI では，システムモデルとして統合コード TASK3D が用いられている．磁場閉じ込め核
融合プラズマは，磁力線がトーラスを稠密に周回することで作られる面（磁気面と呼ばれる）を利用して閉じ込められる（図 3）．磁気面は入れ子構造を作っており，荷電粒子が磁力線に巻き付く性質を利用してプラズマを閉じ込めている．荷電粒子は磁力線方向に素早く動くため，磁気面上で物理量が一定であるとする仮定を置くことができ，図 3 のような複雑なトーラス形状のプラズマでも磁気面についての一次元流体に近似して取り扱うことができる．これは核融合
プラズマのマクロな取り扱いの一例であり，この一次元流体を特徴付けるパラメータとして，（磁気軸），1 がプラズマの最も外側の磁気面規格化小半径 ρ を用いる．ρ は，0 が中心の磁気面
（径方
に対応する
（図 3）．TASK3D もこの一次元流体近似を用いており，磁気面を横切る方向

図 3．磁力線によって作られる磁気面の例（左）とその入れ子構造（右）．

核融合プラズマ制御に向けたデータ同化システムの開発

51向）の粒子と熱の輸送現象に問題の焦点を絞っている．様々な物理過程から統合的に評価される拡散係数や対流速度，ソース項を用いて，粒子と熱の一次元輸送方程式を連立して解き，各粒子種の密度と温度の時間発展を計算することができる．
（LHD）（Takeiri et al., 2017）における粒子・熱これまでに，ASTI を用いて大型ヘリカル装置．LHD は，岐阜県土岐市にある世輸送のデータ同化を行ってきた（Morishita et al., 2022, 2020）
界最大級の超伝導磁場閉じ込め核融合実験装置であり，多数の放電データが蓄積されている．
密度と電子温度，イオン温度の径方向分布の時系列観測データを 40–500 ms 周期で同化することで，乱流モデルパラメータ等の推定を行っており，高い推定精度と妥当性が確認されている．
ASTI は，状態ベクトルや観測ベクトルを調整することで様々な推定問題に柔軟に対応することができるため，核融合プラズマの強力な解析システムとなることが期待される．
3.

データ同化を用いた核融合プラズマの制御

本節では，データ同化を用いた核融合プラズマ制御についての取り組みを紹介する．第 1 節
でも述べた通り，一般的なデータ同化は観測データから状態ベクトルを推定する技術であり，制御の要素は含んでいない．しかしながら，核融合分野での最終的な目標は核融合プラズマの制御にあるため，データ同化システムによる核融合プラズマの制御を実現するには，データ同化の枠組みを予測制御を含んだ形に拡張する必要がある．そこで，著者らは観測情報によるシ
ステムモデルの更新と目標状態を実現する制御入力の推定とを統合した新たなデータ同化のフを開発した．
レームワーク（DACS）（逐次ベイズフィルタ）
予測制御を含んだデータ同化は，特に以下の 2 点が通常のデータ同化と異なる．
（1）設定した目標状態に実システムの挙動が近づくように制御入力の推定を行う．
（2）システムの挙動を予測し，制御入力を推定し続けなければならないので，数値空間上の状態の時刻が常に実空間の時刻よりも先に進んでいる．
また，システムモデルは，状態の時間発展を実時間よりは速く計算できるものの，制御を推定するために何度も計算をやり直せるほど軽量ではないことを想定している．そのため，制御ア
ルゴリズム全体で，できる限り予測計算の回数を少なくしなければならない．そこで，DACS
フレームワークでは，後述する 3 つのフィルタ：z フィルタ，u フィルタ，y フィルタを導入すが実現され，y フィルタによって上記（2）の状況における観測
る．z フィルタによって上記（1）の同化が可能となる．また，u フィルタおよび y フィルタを用いることで，計算区間に重なりのない（一回の予測計算のみの）制御アルゴリズムを構築することが可能となる．
3.1 予測制御を含んだデータ同化フレームワーク
時間間隔 Δty 毎に観測データ y が得られ，時間間隔 Δtz おきに目標状態 z を実現する制御
入力 u∗ を推定する状況を考える．簡単のため，Δty = nΔtz (n ∈ N ) であるとする．まず初めに，時間の記載が煩雑になるのを防ぐため，以下のような時間の記法を導入する．
(3.1)

ti,j = t0,0 + iΔty + jΔtz

ここで，t0,0 は初期時刻，i, j は整数であるとする．Δty = nΔtz を仮定しているため，ti,n = ti+1,0である．また，A(i,j) と表記した場合，時刻 ti,j における変数 A を表すものとする．特に，yiのように j を省略して書いた場合は，j = 0（ti,0 ）であると約束する．また，複数時刻をまとめて，A(0,0):(i,j) ≡ {A(0,0) , A(0,1) , A(0,2) , . . . , A(1,0) , . . . , A(i,j) } のように表記する．

統計数理

52

第 71 巻 第 1 号 2023

時刻 ti,j における状態ベクトルを以下の形で定義する．

x̃(i,j)(3.2)
x(i,j) =u(i,j)
ここで，x̃(i,j) は系の状態やモデルパラメータなどが含まれる部分であり，通常のデータ同化
において状態ベクトルとして定義されるものである．u(i,j) は，制御入力となる変数を成分とするベクトルであり，時刻 ti,j−1 から時刻 ti,j = ti,j−1 + Δtz の間の系の時間発展を決める．また，制御入力の推定時には目標状態 z(i,j) が与えられているものとする．目標状態は，予め時系列データとして設定しておいてもよいし，制御の進行に合わせて決定してもよい．目標状
態 z(i,j) を実現するために推定された制御入力ベクトルを，確率変数 u(i,j) と区別して u∗(i,j) とする．
制御推定を取り入れるために，以下の 4 つの式から構成される状態空間モデルを考える．
(3.3)

x(i,j+1) = f(i,j+1) (x(i,j) , v(i,j+1) )

(3.4)

z
z(i,j) = hz(i,j) (x(i,j) , w(i,j)
)

(3.5)

u
u∗(i,j) = H u x(i,j) + w(i,j)

(3.6)

yi

= hyi (x(i,0) , wiy )式（3.3）は，システムモデルであり，システムノイズ v(i,j+1) を考慮した x(i,j) → x(i,j+1) の時
間発展を記述する．ここで，この時間発展の間，u は時間変化しないものとする．すなわち，
式（3.3）を x̃ と u について分けて書くと，(3.7)

x̃(i,j+1) = f˜(i,j+1) (x̃(i,j) , u(i,j+1) , ṽ(i,j+1) ),

(3.8)

u
u(i,j+1) = u(i,j) + v(i,j+1)となる．ṽ ，v u はそれぞれ，x̃，u に対するシステムノイズである．また，u(i,j) へのシステム
uは，時間発展計算の前に付与するものとした．式（3.4）は，状態ベクトル x と目
ノイズ v(i,j+1)zを考える．
標状態ベクトル z との関係を表現したモデルである．ここで，目標状態ノイズ w(i,j)∗
式（3.5）は，状態ベクトル x と制御入力ベクトル u との関係を表現したモデルであり，状態ベu
クトルの定義（3.2）
から，変換行列 H u を用いて書かれる．ここでも，制御入力ノイズ w(i,j)
をと（3.5）
をそれぞれ，目標状態モデル，制御入力モデルと呼ぶことにする．式
考える．式（3.4）と同様の観測モデルであり，観測ノイズ wiy を考慮する．この状態空間モデル（3.6）は，式（2.2）
は，状態ベクトルを制御入力を含んだ形に拡張し，第 2 節で紹介した逐次ベイズフィルタの状（，2.2）に目標状態モデル（3.4）と制御入力モデル（3.5）を加えたモデルである．
態空間モデル（2.1）
開発した DACS フレームワークは，状態分布に対する 4 つの操作：予測，z フィルタ，u フィ
ルタ，y フィルタで構成される．各操作における状態分布の変化を前もって以下に示しておく．
予測：
p(x(i,j) |y0:k , u∗(0,1):(i,j) ) → p(x(i,j+1) |y0:k , u∗(0,1):(i,j) ).

(3.9)
z フィルタ：(3.10)

p(x(i,j+1) |y0:k , u∗(0,1):(i,j) ) → u∗(i,j+1) = E(u(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) ).

u フィルタ：(3.11)

p(x(i,j+1) |y0:k , u∗(0,1):(i,j) ) → p(x(i,j+1) |y0:k , u∗(0,1):(i,j+1) ).

核融合プラズマ制御に向けたデータ同化システムの開発

53

y フィルタ：(3.12)

p(x(i,j) |y0:k , u∗(0,1):(i,j) ) → p(x(i,j) |y0:k+1 , u∗(0,1):(i,j) ).

ここで，式（3.12）
では，後述する予測ホライズン h(≥ Δtz ) を用いて，ti,j = tk+1 + h の時間関係が成立しているものとする．以下，各操作について詳しく説明する．
今，時刻 ti,j の分布 p(x(i,j) |y0:k , u∗(0,1):(i,j) ) が所与であるとする．ここで，最新の観測時
刻 tk (< ti,j ) までの観測情報 y0:k と時刻 ti,j に至るまでの制御入力 u∗(0,1):(i,j) の情報が，この分布に反映されているものとする．予測（3.9）は，p(x(i,j) |y0:k , u∗(0,1):(i,j) ) から Δtz 先の分布
p(x(i,j+1) |y0:k , u∗(0,1):(i,j) ) を計算する操作であり，予測ステップ（2.3）と同じである．予測分布p(x(i,j+1) |y0:k , u∗(0,1):(i,j) ) には，時刻 ti,j における系の状態とモデルパラメータ（x̃(i,j) ）が持つ不確実性，時間発展計算に伴う不確実性，制御入力 u(i,j+1) による不確実性が含まれる．
z フィルタ（3.10）
は，予測分布と目標状態 z(i,j+1) から最適制御入力 u∗(i,j+1) を推定する操作である．本稿では，目標状態 z(i,j+1) を実現する最適制御入力 u∗(i,j+1) を条件付き確率分布p(u(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) ) による期待値で定義する．すなわち，
∗(3.13)
u(i,j+1) = u(i,j+1) p(u(i,j+1) | y0:k , u∗(0,1):(i,j) , z(i,j+1) )du(i,j+1)である．確率分布 p (u(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) ) は，式（3.2）より，(3.14)p (u(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) )= p(u(i,j+1) , x̃(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) )dx̃(i,j+1)
= p(x(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) )dx̃(i,j+1)

で あ る た め ，状 態 分 布 p(x(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) ) の x̃(i,j+1) に つ い て の 周 辺 化 によ り 求 め る こ と が で き る ．こ の 状 態 分 布 p(x(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) ) は ，予 測 分 布p(x(i,j+1) |y0:k , u∗(0,1):(i,j) ) に，目標状態モデル（3.4）を用いて目標状態 z(i,j+1) を同化した分布である．すなわち，z フィルタは，時刻 ti,j+1 に目標状態 z(i,j+1) が実現したとして，その際に行われた制御を推定していると解釈できる．また，この制御推定は，時刻 ti,j における系の状態やモデルパラメータが持つ不確実性に対してロバストな制御推定であると言える．アンサ
ンブル近似による手法を用いる場合には，予測分布 p(x(i,j+1) |y0:k , u∗(0,1):(i,j) ) のアンサンブルに z(i,j+1) を同化し，できた p(x(i,j+1) |y0:k , u∗(0,1):(i,j) , z(i,j+1) ) のアンサンブルの u の平均を計算すれば u∗(i,j+1) が求まる．そのため，計算コストはフィルタ計算のコストのみである．
u フィルタ（3.11）は ，推 定 し た 制 御 入 力 u∗(i,j+1) に よ り 制 御 を 行 っ た 場 合 に 予 測 さ れる 状 態 分 布 p(x(i,j+1) |y0:k , u∗(0,1):(i,j+1) ) を 求 め る 操 作 で あ る ．u フ ィ ル タ は ，予 測 分 布p(x(i,j+1) |y0:k , u∗(0,1):(i,j) ) に，制御入力モデル（3.5）を用いて u∗(i,j+1) を同化することで実現∗できる．フィルタ後の分布 p(x(i,j+1) |y0:k , u(0,1):(i,j+1) ) を u フィルタ分布と呼ぶことにする．u
フィルタ分布は予測分布に u∗(i,j+1) の情報を付与した分布であり，予測分布に含まれる不確実の不確実性が取り除かれた（低減された）分布である．すなわち，性のうち，制御入力（u(i,j+1) ）u フィルタ分布は，主に状態 x̃(i,j+1) の不確実性を表現した分布であり，Δtz 前の状態 x̃(i,j)
および時間発展計算が持つ不確実性が含まれている．この不確実性は，y フィルタ時に観測情
報を未来に伝える働きをする．アンサンブル近似を用いる場合に，予測計算を行う前の分布
p(x(i,j) |y0:k , u∗(0,1):(i,j) ) の各アンサンブルメンバーの u(i,j+1) 部に u∗(i,j+1) を代入して予測計算をやり直しても，p(x(i,j+1) |y0:k , u∗(0,1):(i,j+1) ) を推定できる．こちらの方法の方が正確である
が，同じ時間区間の予測計算を 2 回行うことになる．u フィルタは，この 2 度目の予測計算を

統計数理

54

第 71 巻 第 1 号 2023

省く操作であり，計算コストを無視できないモデルでモデル予測制御を実現するためには重要な操作である．予測 →z フィルタ →u フィルタによって，目標状態 z(i,j+1) を実現する制御
入力 u∗(i,j+1) を推定するとともにその制御を行った場合の状態分布を推定できる．この操作のセットによって，さらに先の制御推定に進むことができる．
y フィルタ（3.12）
は，新たな観測 yk+1 を u フィルタ分布に反映させるためのフィルタである．制御を考える場合，観測時点 tk+1 からの制御 u∗(k+1,1) が時刻 tk+1 になる前に決定できていなければならないため，数値空間の時刻は実空間の時刻よりも先に進んでいる．y フィルタの役割は，時間を隔てた最新の u フィルタ分布へ観測情報を反映させることである．ここで，この時間の隔たりを予測ホライズン h と定義する．すなわち，
h = ti,j − tk+1

(3.15)

である. ここで，ti,j は最新の u フィルタ分布の時刻とする．制御の都合から h ≥ Δtz
で な く て は な ら な い ．こ の 時 間 を 隔 て た 最 適 化 操 作 は ，観 測 時 点 と 最 新 の u フ ィ ル
タ 分 布 の 時 点 に お け る 状 態 ベ ク ト ル の 同 時 確 率 分 布 p(x(i,j) , x(k+1,0) |y0:k , u∗(0,1):(i,j) ) を
考 え る こ と で 実 現 で き る ．す な わ ち ，こ の 同 時 分 布 p(x(i,j) , x(k+1,0) |y0:k , u∗(0,1):(i,j) ) に 観
測 yk+1 を同化することで p(x(i,j) , x(k+1,0) |y0:k+1 , u∗(0,1):(i,j) ) を得て，更に xk+1 について周 辺 化 す る こ と で y フ ィ ル タ 分 布 p(x(i,j) |y0:k+1 , u∗(0,1):(i,j) ) を 得 る ．ア ン サ ン ブ ル 近 似を 用 い る 場 合 に は ，p(x(k+1,0) |y0:k , u∗(0,1):(k+1,0) ) の ア ン サ ン ブ ル を 保 存 し て お き ，最 新の u フ ィ ル タ 分 布 p(x(i,j) |y0:k , u∗(0,1):(i,j) ) の ア ン サ ン ブ ル と 連 結 す る こ と で 同 時 分 布p(x(i,j) , x(k+1,0) |y0:k , u∗(0,1):(i,j) ) のアンサンブルを構成することができる．ここに yk+1 を同 化 し ，時 刻 tk+1 に 対 応 す る 部 分 を 切 り 離 す こ と で 周 辺 化 す れ ば ，y フ ィ ル タ 分 布p(x(i,j) |y0:k+1 , u∗(0,1):(i,j) ) のアンサンブルを手に入れることができる．y フィルタについても，y フィルタ後の観測時刻の分布 p(x(k+1,0) |y0:k+1 , u∗(0,1):(k+1,0) ) のアンサンブルから予測計
算をやり直して p(x(i,j) |y0:k+1 , u∗(0,1):(i,j) ) のアンサンブルを求めた方が正確である．y フィル
タは，観測を状態分布に反映させると同時にこの予測計算をやり直すコストを省いている．u
フィルタと y フィルタにより，予測計算区間の重なりを完全に排除することができ，各時間区間について一度のみの予測計算で制御アルゴリズムを構築できる．予測計算のコストが無視できないためにイタレーションを用いるモデル予測制御に使用できなかったシミュレーションモ
デルも，DACS フレームワークではシステムモデルとして用いることができる可能性がある．
それに加えて，観測の同化による最適化を前提として，より簡約なモデルを採用することができるため，予測計算の更なる高速化が期待できる．
3.2 制御アルゴリズム
パラメータ n と予測ホライズン h を設定すると，具体的な制御アルゴリズムを構築することができる．もちろん，h = Δtz とする場合が最も素早く観測を制御推定に反映できるが，計算コストの面で余裕がなくなる可能性がある．このため，n や h は予測計算やフィルタ計算，並列計算の通信速度や観測情報の処理速度などを考慮して，慎重に決定する必要がある．本稿では，n = 3 とし，h = Δty = 3Δtz の場合を検証する．この場合の時間関係を図 4 に示す．
n = 3，h = Δty を考えると，k = i − 1 であり，式（3.9）（3.12）は以下のようになる．
–
予測：
p(x(i,j) | y0:i−1 , u∗(0,1):(i,j) ) → p(x(i,j+1) | y0:i−1 , u∗(0,1):(i,j) ).

(3.16)z フィルタ：(3.17)

p(x(i,j+1) | y0:i−1 , u∗(0,1):(i,j) ) → u∗(i,j+1) = E(u(i,j+1) | y0:i−1 , u∗(0,1):(i,j) , z(i,j+1) ).

核融合プラズマ制御に向けたデータ同化システムの開発

55

図 4．制御アルゴリズムにおける時間関係（Δty = 3Δtz ，h = Δty の場合）．

図 5．状態分布が情報を取り込む流れ（Δty = 3Δtz ，h = Δty の場合）．

u フィルタ：(3.18)

p(x(i,j+1) | y0:i−1 , u∗(0,1):(i,j) ) → p(x(i,j+1) | y0:i−1 , u∗(0,1):(i,j+1) ).

y フィルタ：(3.19)

p(x(i+1,0) | y0:i−1 , u∗(0,1):(i+1,0) ) → p(x(i+1,0) | y0:i , u∗(0,1):(i+1,0) ).

図 5 は，この制御アルゴリズムにおいて状態分布が目標状態と観測情報を取り込んでいく流れ)が所与であるとすると，を表したものである．状態分布 p(x(i,0) |y0:i−1 , u∗(0,1):(i,0)（図中ハート）∗
予測 →z フィルタ →u フィルタの繰り返しにより，制御入力 u(i,1):(i+1,0) を推定しながら Δty
先の u フィルタ分布 p(x(i+1,0) |y0:i−1 , u∗(0,1):(i+1,0)（図中ダイヤ）)まで計算を行う．観測 yi が得

56

統計数理

第 71 巻 第 1 号 2023

られると，最新の u フィルタ分布は y フィルタにより p(x(i+1,0) |y0:i , u∗(0,1):(i+1,0)（図中スペー))ド）へと更新される．ここで，観測時刻 ti の y フィルタ分布 p(x(i,0) |y0:i−1 , u∗(0,1):(i,0)（所与とした分布，図中ハート）が観測情報を未来（ti+1 ）に伝えるために用いられる．y フィルタ分布p(x(i+1,0) |y0:i , u∗(0,1):(i+1,0) ) から更に ti+2,0 へと制御推定が進んでいく．図 5 の 3 次元的な流（x と y で時刻がずれている点を考えなければ）通常の逐次ベイズ
れを x-y 平面に射影すると，フィルタのデータ同化の流れが現れる（x-y 平面上実線）．このことから，DACS フレームワー
クは，通常の逐次ベイズフィルタの予測ステップを予測＋ z フィルタ＋ u フィルタの一連の操作に拡張したものであると言える．この拡張によって，データ同化の枠組みに目標状態時系列
を用いて実システムに働きかける仕組みが取り込まれる．
状態分布のアンサンブル近似により上記のアルゴリズムを実現する．時刻 t0,0 から観測が取得できるとすると，制御の計算手続きは以下のようになる．
• 初期アンサンブルメンバーを生成し，時刻 t0,0 の y フィルタ分布として保存しておく．
• For i = 0, . . . , imax（観測周期）：（1）制御入力 u∗(i,1):(i+1,0) と時刻 t(i+1,0) までの u フィルタ分布の推定：
For j = 0, 1, 2（制御推定周期）：（ti,3 = ti+1,0 ）．
（a）時刻 ti,j から ti,j+1 まで予測計算を行う（a）で得た予測分布を用いて z フィルタを実行し，制御入力（b）目標状態 z(i,j+1) と 1u∗(i,j+1) を推定する．
（a）で得た予測分布を用いて u フィルタを実行する．
（c）制御入力 u∗(i,j+1) と 1（2）時刻 t(i,0) に観測を行い，観測ベクトル yi を取得する．
（3）観測 yi を用いた時刻 t(i+1,0) の u フィルタ分布の更新：（a）時刻 ti+1,0 の u フィルタ分布と保存しておいた時刻 ti,0 の y フィルタ分布を連結する．
y フィルタにより，3（a）で生成した連結アンサンブルに観測 yi を同化する．
（b ）（c）連結アンサンブルから時刻 ti+1,0 の部分を y フィルタ分布のアンサンブルとして取り出す．このアンサンブルは次の y フィルタに利用するために保存しておく．
全てのフィルタを EnKF により実装する場合，目標状態モデルおよび観測モデルが線形・ガ
ウス分布となる以下のような状態空間モデルを考える．
(3.20)

x(i,j+1) = f(i,j+1) (x(i,j) , v(i,j+1) ), v(i,j+1) ∼ N (0, Q(i,j+1) )(3.21)zzzz(i,j) = H z x(i,j) + w(i,j), w(i,j)∼ N (0, R(i,j))(3.22)uuuu∗(i,j) = H u x(i,j) + w(i,j), w(i,j)∼ N (0, R(i,j))(3.23)

yi = H y x(i,0) + wiy , wiy ∼ N (0, Riy )ここで，目標状態モデル（3.21）と観測モデル（3.23）を行列 H z と H y を用いて書いた．また，ノzu，R(i,j)，Riy は，それイズは全て平均 0 のガウス分布に従うものとした．行列 Q(i,j+1) ，R(i,j)
ぞれのノイズの共分散行列である．これらの共分散行列は，制御性能を決める重要なハイパー
パラメータである．本稿では簡単のため，これらの行列は全て対角行列であると仮定する．
システムノイズの共分散行列 Q(i,j+1) は，状態分布が持つ不確実性を制御する．特に制御入力 u に対するシステムノイズは，u∗ を重点的に探す領域の広さを決めるため重要である．目標状態の変化が緩やかな場合や車のハンドルのように有限の速度で制御入力が変化する場合には，比較的小さな分散で良い．一方で，制御入力に大きな変化（ジャンプ）
が必要となるような制御には，大きな分散が必要となる．行列 Q(i,j+1) は，システムモデルの洗練度や目標状態時
系列の変化の度合い，制御機器の制約などに応じて設定する必要がある．制御入力 u に対する

核融合プラズマ制御に向けたデータ同化システムの開発

57

システムノイズは毎予測計算前に，x̃ に対するシステムノイズは毎 y フィルタ計算後に付与するものとする．
目標状態ノイズは，一回の制御推定における目標状態の実現精度を表しており，共分散行列zz
R(i,j) により z フィルタの動作を調整することができる．R(i,j)の分散成分が小さいと，z フィ
ルタ時に予測分布が目標状態により強く引っ張られる．逆に分散成分が大きい場合には，予測分布を目標状態方向に引っ張る力が弱まり，目標状態時系列に対して，状態変数の変化が緩やかになることが期待される．システムモデルの不確実性が大きい場合には，目標状態ノイズを大きく取ることで状態変化の速度を遅くし，システムモデルと実システムの差異から生じるzを予測分布の共分散行列に比例す
制御推定のズレを緩和することができる．ここでは，R(i,j)
る形(3.24)z(R(i,j))ll = rz2 (H z V(i,j) (H z )T )ll

で与える．rz は比例係数であり，V(i,j) は予測分布のアンサンブルの共分散行列である．また，( )ll は，行列の第 l 対角成分を表し，T は行列の転置を表す．
制御入力ノイズは，制御入力時の誤差を表現しているが，多くの場合十分小さく取って問題ないと考えられる．制御に対する影響として，制御入力 u の変化速度を z フィルタとは独立に調整できるが，結果的な影響としては目標状態ノイズが持つ影響とほとんど同じである．本稿
では，制御入力ノイズは十分小さく取るものとする．EnKF の場合，u の次元がアンサンブル
メンバー数に比べて十分に小さいならば，分散を 0 としても問題なく u フィルタ計算を行うことができる．PF の場合には，制御入力ノイズを小さく取りすぎると粒子の退化が顕著になるため，慎重に設定する必要がある．
観測ノイズの役割は，通常の逐次ベイズフィルタと同じである．ASTI では，観測ノイズを観測値と予測値の差に比例する形
2(3.25)(Riy )ll = ry2 yi − H y x̂(i,0) l ,で与えると，安定的にデータ同化を行えることがわかっている（Morishita et al., 2022, 2020）．
ここで，ry は比例係数であり，x̂(i,0) は観測時点に対応する状態分布 p(x(i,0) |y0:i−1 , u∗(0,1):(i,0) )のアンサンブル平均である．また，( )l はベクトルの第 l 成分を表す．
3.3 仮想プラズマを用いた制御数値実験
3.2 節で構築した制御アルゴリズムを ASTI に実装し，TASK3D によって数値的に生成した
仮想プラズマを制御する数値実験によりその有効性を検証した．全てのフィルタは EnKF により実装した．目標状態としてプラズマ中心（ρ = 0）における密度と電子温度を取り，密度，電子温度，イオン温度の径方向分布が Δty = 0.6 秒おきに観測できると仮定する．本稿では，電子密度とイオン密度は等しいものとし，合わせて “密度” と表記する．仮想プラズマは軽水素
プラズマを仮定し，LHD の磁場配位（磁気軸大半径 3.6m，磁場強度 2.85T）を使用する．表 1に，本数値実験での状態ベクトル，目標状態ベクトル，観測ベクトルを構成する変数をまとめた．Mi は各変数がベクトル中で持つ次元であり，σinit は初期分布の標準偏差，σQ はシステ
（ρ = 0.0, 0.1, 0.2, . . . , 1.0)
ムノイズの標準偏差である．Mi が 11 である状態変数は，径方向分布を持っている．標準偏差が % 単位で表記されているものは，標準偏差をその時の状態分布の平均値に比例する形で決定する．また，表中で日本語の説明に ∗ が付いているものは，その変数にかかる数因子を状態変数とすることを意味する．例えば，乱流粒子拡散係数（DTB ）であれ
TB
ば，TASK3D 内で D = dmodel × d のように状態変数 d による調整を受ける．ここで，dmodelは通常のシミュレーションで用いられるモデルである．x̃ は主に乱流による拡散係数の不確実

58

統計数理

第 71 巻 第 1 号 2023

表 1．状態ベクトル (x̃, u)，目標状態ベクトル z，観測ベクトル y の構成．

図 6．各状態変数の役割．

性を考慮して構成した．制御入力ベクトル u は，電磁波加熱の一つである電子サイクロトロンの入力パワー PECH とプラズマの端での中性粒子密度を調整する数因子 ξn
周波数加熱（ECH）
により定義する．PECH は電子の熱輸送における加熱項を決定し，主に電子温度分布に影響する．ξn は粒子輸送のソース項を決定し，主に密度分布に影響する．これらの各状態変数の役割を図 6 にまとめる．PECH は 0.5–5.0 MW の間を 0.5 MW 刻みで取る離散値であるとし，ξn は連続値であるとする．z フィルタ時には全ての状態変数が連続変数として取り扱われるため，∗
PECH については，式（3.13）の右辺で計算される期待値に最も近い離散値を PECHとする．制御1推定の周期は Δtz = 3 Δty = 0.2 秒である．
仮想プラズマの初期状態は，初期制御入力 u(0,0) = (0.5 MW, 1.0) に対する定常状態であるとする．初期状態分布の平均は，密度，温度については仮想プラズマの定常状態の値を用い，∗付き変数については 1.0（通常のシミュレーションの値）であるとする．目標状態の変化として，
密度と電子温度を (1 × 1019 m−3 , 2 keV) から同時に上昇させていき，(2 × 1019 m−3 , 5 keV) まで（a）（b）の波線）．ここで，eV（電子ボル
上昇したところで状態を維持する時系列を考える（図 7◦ト）はプラズマの温度を表す単位であり，1eV は約 11300 C，keV はその 103 倍で約一千万 ◦ Cである．数値実験として，システムモデルに不正確な部分がある状況下で，観測情報を用いてシステムモデルの挙動を実システムに近づけながら制御を行う問題を考える．そのためこの数

核融合プラズマ制御に向けたデータ同化システムの開発

59

図 7．仮想プラズマを制御する数値実験の結果．

値実験では，仮想プラズマの乱流粒子拡散係数 DTB と電子の乱流熱拡散係数 χTBに 1/n に比e
TB
TB
例する項を追加する．D と χe はそれぞれ，密度勾配と電子温度勾配にかかり，径方向外向きの乱流粒子束と電子乱流熱流束を生む．粒子拡散係数 DTB について，ASTI では，定数モ
デル(3.26)DTB = dconst × dを仮定するが，仮想プラズマ内では，1/n の項を加えた形

1(3.27)DTB = dAlc 1 +nを仮定する．ここで，dconst = 1.4，dAlc = 0.4 とする．電子の乱流熱拡散係数 χTBについては，e
ASTI では，gyro-Bohm モデルと呼ばれるモデル(3.28)

χTB
= CgB
e

Te ρi
× ce
eB a

を仮定するが，仮想プラズマ内では，1/n の項を加えた形


Te ρi
1
1
(3.29)χTB=
CgB+ CAlce2eB anを用いる．ここで，CgB = 1.5，CAlc = 1.61 である．B ，ρi ，a はそれぞれ，磁場強度，イオン
ラーモア半径，プラズマ小半径である．また，dconst ，dAlc ，CgB ，CAlc の値は，複数の実験結果

60

統計数理

第 71 巻 第 1 号 2023に基づいて行われた先行研究（Sakai et al., 2015; Wakasa et al., 2010）から妥当な値を設定した．
図 7 に，数値実験の結果を示す．この数値実験では，rz = 0.6，ry = 0.8 とした．また，アンサンブルメンバー数は 420 とした．この状況設定においては，アンサンブルメンバー数を 200程度まで少なくしても制御精度は大きく低下しないことを確認している．図 7 の（a）と（b）はそ
れぞれ，プラズマ中心における電子温度 Te と密度 n の時間変化を表しており，仮想プラズマ（実線）が目標状態（破線）を追随できているのがわかる．また，（c）と（d）はその際に推定された（b）における，“u-filtered prediction” は u フィルタ
制御入力の値を示している．ここで，（a），の平均と 1 標準偏差であり，“Obs filtered” は対応する時分布（p(x(i,j+1) | y0:i−1 , u∗(0,1):(i,j+1) )）刻の観測まで同化した分布 p(x(i,0) | y0:i , u∗(0,1):(i,0) ) の平均である．（c）における，“Prediction”
は予測ステップ後の制御入力分布であり，“z-filter estimate” は z フィルタにより推定された制である．また，実線 “Actual control input” は，実際に模擬プラズマの制御に
御入力の値（u∗ ）（b）
どちらにおいても，変化区間（2.4–4.8 秒）では，目標状態と用いられた入力値である．（a），
仮想プラズマの差が広がっていることがわかる．これは，システムモデルと仮想プラズマのモ
デルの差異による影響であり，観測情報を用いてモデルの不足分を押さえ込みながら制御がと（f）
はモデルの不足分を補うために調整された乱流輸送のモデルパラメータ進められる．（e）
を示している．
（g）と（h）は，仮想プラズマ内で式（3.29），（3.27）を用いて計算（径方向 2 点分）と DTB の値
（点）と，ASTI 内で ce および d によって調整を受けた χTBと DTB のされる χTBee
アンサンブルにおける平均値（実線および破線）を比較した図である．序盤では，仮想プラズマと ASTI で値が大きく異なっている場合でも，観測の同化による ce および d の調整によって，両者の値が近づいているのがわかる．
図 8 は，時刻 ti = 4.2 秒の状態分布を表すアンサンブルメンバーが各ステップでどのような変化を受けるかを示した図である．ここで用いたアルゴリズムでは，h = Δty としているため，観測の時刻におけるアンサンブルは，予測分布 p(x(i,0) | y0:i−2 , u∗(0,1):(i−1,2) ) から
p (x(i,0) | y0:i−2 , u∗(0,1):(i−1,2) )(3.30)→ p(x(i,0) | y0:i−2 , u∗(0,1):(i,0) )(3.31)→ p(x(i,0) | y0:i−1 , u∗(0,1):(i,0) )(3.32)→ p(x(i,0) | y0:i , u∗(0,1):(i,0) )の様に 3 回の変化を受ける．
（a）と（b）
は，変化（3.30）の図であり，予測分布（“Prediction”）から
u フィルタ分布（“u-filtered”）
への変化を表している．（a）はプラズマ中心における電子温度と密度，（b）は制御入力変数についての散布図である．（a）の散布図における “Target” は目標状態であり，z フィルタ＋ u フィルタにより u フィルタ分布が目標状態に近づいているのが分かる．また，制御入力ノイズを小さく取ることで，（b）における u フィルタ分布のアンサンブルのばらつきを実用上問題ない程度まで小さくしている．Δty 前（ti−1 = 3.6 秒）の観測が得られで，u フィルタ分布（“u-filtered”）は y フィルタ分布（“y-filtered (yi−1 )”）ると，y フィルタ（3.31）と（d）はこの変化を表した図であり，（c）はプラズマ中心における電子温へと更新される．（c）
度と密度，（d）は乱流輸送のモデルパラメータ (ce (ρ = 0.2)，d(ρ = 0.5)) に注目している．ま（電た，（e）はこの y フィルタ時に利用した時刻 ti−1 = 3.6 秒における状態分布のアンサンブル
子温度と密度）である．時刻 ti = 4.2 秒における観測が得られると，（c）（d）における y フィルに観測情報を伝えるために用いられる．この時の状態
タ分布は最新の u フィルタ分布（4.8 秒）
であり，図 8（f），（g）
のように，y フィルタ分布（“y-filtered (yi−1 )”）は対応
分布の変化が（3.32）
する時刻の観測を同化した分布（“y-filtered (yi )”）へと更新される．また，観測 yi で xi が更新（h）はこの y フィルタで更新された時されるのにともなって時刻 4.8 秒に観測情報が伝わる．

核融合プラズマ制御に向けたデータ同化システムの開発

61

図 8．図 7 の数値実験における時刻 ti = 4.2 s の状態分布を表すアンサンブルメンバーの変化の様子．

刻 ti+1 = 4.8 秒における状態分布のアンサンブル（中心の電子温度と密度）である．
この例では，制御入力ノイズが u に対するシステムノイズよりも十分に小さいため，予測の状態分布における u と x̃ の相関はほステップ前（u に対するシステムノイズを付与した後）
とんど 0 になっている．そのため，y フィルタにより未来の u フィルタ分布を修正する際に uへの影響はほとんどない．しかし，何らかの理由で制御入力ノイズを大きく取る必要がある時には，yi を同化する y フィルタ時に u の分布 p(u(i+1,0) | y0:i−1 , u∗(0,1):(i+1,0) ) も修正を受け，yフィルタ分布の平均が u∗(i+1,0) からずれてしまう可能性がある．粒子フィルタによるアンサン
ブルの退化を防ぎたい時などがこの状況として考えられる．この時には，予測ステップ前にアンサンブルの u 部分を毎回新たにガウシアンから生成し，u と x̃ の相関を断っておくなどの対策が考えられる．

統計数理

62

第 71 巻 第 1 号 2023

最後にノイズを決めるパラメータ rz と ry について言及する．定常区間（4.8 秒-）では，
0.3 ≤ rz , ry ≤ 1.0 の幅広い範囲で高い精度（観測と目標状態の平均平方二乗誤差率（RMSPE）が5% 以内）となることを確認している．一方で，変化区間（2.4–4.8 秒）での制御精度は rz に影響されやすい．モデルの不確実性が大きい場合には制御推定の精度が落ちるため，rz を大きく取り，状態の変化速度を遅くした方が安定した制御を行うことができる．また，ry を小さく取りすぎると，計測誤差が大きいような観測値の情報も状態分布に強く反映されてしまい（過適合），システムモデルの予測能力を低下させてしまうことが考えられる．そのため，システム
モデルや観測値の不確実性が大きい場合や未知の場合には，rz , ry をある程度大きい値（≥ 0.5）に選ぶのが安全である．システムモデルの適応能力は観測ノイズ（ry ）の他にも x̃ に対するシ
ステムノイズによって，調整することができる．ただし，x̃ に対する大きなシステムノイズ ṽ
は，x̃ に含まれる状態変数間の相関を薄め，y フィルタ時の時間を隔てた最適化を弱める．また，取れるアンサンブルメンバーが限られる状況下では制御入力推定にも悪影響する．そのため，ṽ に過度に大きな分散を設定しない方が安全である．
4.

まとめと今後の課題

本稿では，核融合分野でのデータ同化システムの構築に関して，制御に重きを置いて説明した．制御を含んだデータ同化のフレームワーク（DACS）を開発しており，核融合プラズマに限らず多くの制御問題に適用可能である．システムモデルが持つ不確実性を観測による逐次最適化により押さえ込むことができるため，不確実性の大きいシステムモデルを用いたモデル予測制御が実現できる．また，DACS フレームワークを用いると，予測区間に被りのない並列性に優れた制御アルゴリズムを構築することができるため，比較的重いシミュレーションコードでもシステムモデルとして用いることができる．
この手法の課題として，システムモデルの適応の限界が挙げられる．これは，観測の同化周期に対して実システムが速く変化する場合に，実システムの変化にモデルパラメータの最適化が追いつかなくなる状況である．システムモデルによる実システムの近似精度が大きく低下す
ると，誤った制御入力を推定することになり，制御精度が低下する．この問題は，表現能力の高いモデルを用いて柔軟に実システムの変化に適応できるようにするか，観測周期を上げて状態ベクトルの最適化頻度を上げることなどで解決することができる．ただし，どのような状況に陥るか想定できない場合も多いので，状況に合わせて観測周期やノイズのパラメータを調整できる仕組みを導入することが理想的である．
また，実際の制御の場面では，系の状態の安定性などを考慮に入れて制御を行わなければならない．特に核融合プラズマでは，多くの不安定状態が存在し，制御を誤ればプラズマが消えてしまう状況も起こり得る．このようなプラズマの持続性能に直結する課題に対するデータ駆動的な研究も進められており（Yokoyama et al., 2020），実際の制御にはこうした研究との連携が重要である．現在，本稿で言及した制御アルゴリズムを利用して，実際の核融合プラズマを制御する実証実験を計画中である．
謝

辞，核融合科学研
本研究は統計数理研究所共同研究（2021-ISMCRP-2005，2022-ISMCRP-2026）
究所共同研究（NIFS20KLPT007，NIFS22KAPT008），科研費（JP21J14260）の助成を受けたものです．また，本研究で紹介した計算結果は国際核融合エネルギー研究センター（IFERC）のスーパーコンピュータシステム JFRS-1 により計算されたものです．

核融合プラズマ制御に向けたデータ同化システムの開発

参

考

文

63

献

足立修一 (2014). モデル予測制御の基礎, 日本ロボット学会誌, 32(6), 499–502.
Evensen, G. (2003). The ensemble Kalman filter: Theoretical formulation and practical implementation, Ocean Dynamics, 53(4), 343–367, http://dx.doi.org/10.1007/s10236-003-0036-9.
Gettelman, A., Geer, A. J., Forbes, R. M., Carmichael, G. R., Feingold, G., Posselt, D. J.,
Stephens, G. L., van den Heever, S. C., Varble, A. C. and Zuidema, P. (2022). The future ofEarth system prediction: Advances in model-data fusion, Science Advances, 8(14), p. eabn3488,
http://dx.doi.org/10.1126/sciadv.abn3488.
口知之, 上野玄太, 中野慎也, 中村和幸, 吉田 亮 (2011).『データ同化入門』, 朝倉書店, 東京.
Kapteyn, M. G., Pretorius, J. V. and Willcox, K. E. (2021). A probabilistic graphical model foundation for enabling predictive digital twins at scale, Nature Computational Science, 1(5), 337–347,
http://dx.doi.org/10.1038/s43588-021-00069-0.
Kitagawa, G. (1996). Monte Carlo filter and smoother for non-Gaussian nonlinear state space models, Journal of Computational and Graphical Statistics, 5(1), 1–25, http://dx.doi.org/10.1080/
10618600.1996.10474692.
Liu, M., Fang, S., Dong, H. and Xu, C. (2021). Review of digital twin about concepts, technologies, and industrial applications, Journal of Manufacturing Systems, 58, 346–361, http://dx.doi.org/10.1016/j.jmsy.2020.06.017.
Morishita, Y., Murakami, S., Yokoyama, M. and Ueno, G. (2020). Data assimilation system based on
integrated transport simulation of Large Helical Device plasma, Nuclear Fusion, 60(5), 056001,
http://dx.doi.org/10.1088/1741-4326/ab7596.
Morishita, Y., Murakami, S., Yokoyama, M. and Ueno, G. (2022). ASTI: Data assimilation systemfor particle and heat transport in toroidal plasmas, Computer Physics Communications, 274,
108287, http://dx.doi.org/10.1016/j.cpc.2022.108287.
Murakami, S., Yamaguchi, H., Sakai, A., Wakasa, A., Fukuyama, A., Nagaoka, K., Takahashi,
H., Nakano, H., Osakabe, M., Ida, K., et al. (2015). Integrated transport simulations ofhigh ion temperature plasmas of LHD, Plasma Physics and Controlled Fusion, 57(5), 054009,
http://dx.doi.org/10.1088/0741-3335/57/5/054009.
Sakai, A., Murakami, S., Yamaguchi, H., Wakasa, A., Fukuyama, A., Nagaoka, K., Takahashi, H.,
Nakano, H., Osakabe, M., et al. (2015). Integrated particle transport simulation of NBI plasmas
in LHD, Plasma and Fusion Research, 10, 3403048–3403048, http://dx.doi.org/10.1585/pfr.10.
3403048.
Takeiri, Y., Morisaki, T., Osakabe, M., Yokoyama, M., Sakakibara, S., Takahashi, H., Nakamura, Y., Oishi, T., Motojima, G., Murakami, S., et al. (2017). Extension of the operational regime of the LHD towards a deuterium experiment, Nuclear Fusion, 57(10), 102023,
http://dx.doi.org/10.1088/1741-4326/aa7fc2.
上野玄太 (2019). 粒子フィルタとデータ同化, 統計数理, 67, 241–253.
Wakasa, A., Fukuyama, A., Murakami, S., Miki, M., Yokoyama, M., Sato, M., Toda, S., Funaba,
H., Tanaka, K., Ida, K., Yamada, H., Honda, M. and Nakajima, N. (2010). Integrated transport
simulation of LHD plasmas using TASK3D, Proceedings of 23rd IAEA Fusion Energy Conference(Daejon), P4.029.
Yokoyama, T., Yamada, H., Masuzaki, S., Miyazawa, J., Mukai, K., Peterson, B. J., Tamura, N.,
Sakamoto, R., Motojima, G., Ida, K., et al. (2020). Prediction of radiative collapse in Large
Helical Device using feature extraction by exhaustive search, Journal of Fusion Energy, 39(6),500–511.

64

Proceedings of the Institute of Statistical Mathematics Vol. 71, No. 1, 47–64 (2023)

Development of Data Assimilation System for Fusion Plasma Control
Yuya Morishita1 , Sadayoshi Murakami1 , Masayuki Yokoyama2,5 and Genta Ueno3,4,6
1 Department of Nuclear Engineering, Kyoto University
2 Rokkasho Research Center, National Institute for Fusion Science
3 The Institute of Statistical Mathematics
4 Joint Support-Center for Data Science Research
5 Department of Fusion Science, School of Physical Sciences, The Graduate University for Advanced

Studies, SOKENDAI
6 Department of Statistical Science, School of Multidisciplinary Sciences, The Graduate University for

Advanced Studies, SOKENDAI

We developed a data assimilation system, named ASTI, to analyze and control fusion plasma behavior. Although the goal of ASTI is to control fusion plasmas with highaccuracy, existing data assimilation frameworks do not include control processes. We
therefore developed a data assimilation framework that integrates system model updates
and optimal control-input estimation. The proposed framework provides model predictive
control even when the system model has large uncertainties. In this article, we provide an
overview of ASTI and show the proposed data assimilation framework. We demonstrate
the eﬀectiveness of the framework through a numerical experiment to control virtual fusion
plasma.

Key words: Data assimilation, model predictive control, fusion plasma.