👌

Density Controlによる最適輸送(1)

2024/07/27に公開

最適輸送問題はいかに輸送コストを小さくできるかという古くから知られた問題です。また最適制御問題はシステムをどのように制御すれば望む状態にできるかを考える工学などにおいて重要なテーマです。

最適輸送問題を最適制御問題として扱うことで動的システムの視点を導入することができ、各時刻でどのように輸送を行えばよいのかが分かります[1]。また最適制御問題の解法として確立されたツールであるハミルトン・ヤコビ・ベルマン(HJB)方程式と類似の方程式が動的最適輸送問題の解法において現れるため、最適制御の分野で確立された効率的な解法を使用して最適輸送問題を解くことができます[2]。従って最適制御の分野をよく知る方にとっても最適輸送問題が興味ある対象になると思います。最適輸送問題と最適制御問題の関連するテーマはこのような動機により興味ある分野となっています。

最適輸送問題を最適制御問題として動的に解く方法として流体による定式化が提案されていて[3]、Density Controlなどと呼ばれたりしています。この手法により最適輸送問題を時間的な逐次処理で最適化していくことができます。

この記事から始まる一連の記事は3部構成の予定で最適輸送問題をDensity Controlの手法で解くことについて解説します。パート1であるこの記事ではDensity Controlの概要を説明します。パート2では数値的に解くための具体的なアルゴリズムを解説して、パート3では実際に数値計算をする予定です。

以下では初めに最適輸送問題について簡単に述べます。次に最適制御問題について簡単な例を用いながら説明します。そして最後に流体による定式化を説明します。

最適輸送問題

最適輸送問題の説明をします。例えば与えられた2点に対して、一方からもう一方へ質点を低コストで運ぶ方法を考えましょう。色々なコストが考えられますが単純に距離の2乗だとすれば、最適解は測地線(最短曲線)に沿って運べばよいことになります。

では運ぶべき質点が n 個あって初期配置と最終配置が指定されている場合を考えます。ただし質点は互いに区別できず与えられた最終配置のどこにどの質点を置いてもよいものとします。このような問題はある分布を別の分布へコストが最小となるように移動させる問題であり、どのように動かすのが最適なのかはすぐには分かりません。これが、分布が離散的な場合の最適輸送問題です。

分布が連続な場合にコストを距離の2乗で定めるときの最適輸送問題はMonge-Kantorovich質量輸送問題と呼ばれ、正確な定式化は以下のようになります。
\mathbb{R}^dにおいて2つの密度分布\rho_0,\rho_T
\int_{\mathbb{R}^d}\rho_0(x)dx^d=\int_{\mathbb{R}^d}\rho_T(x)dx^d=1
を満たすものが与えられているとき、滑らかな写像M:\mathbb{R}^d\to \mathbb{R}^dで以下の2つを満たすものを求める問題です。
(i) 任意の可測集合A\subset\mathbb{R}^dに対して
\int_{x\in A}\rho_T(x)dx^d=\int_{M(x)\in A}\rho_0(x)dx^d
を満たす。
(ii)
\int_{\mathbb{R}^d}|M(x)-x|^p\rho_0(x)dx^d
を最小にする。
特に p=2 のとき、L^2-MKP (Monge-Kantorovich Problem)と呼ばれます。写像 M を求める方法は多くの研究により提案されています。

最適制御問題

最適制御問題とは、何らかの時間変化するシステムにおいて与えられたときに、例えばコストなどの性能指標を最小化(または最大化)するために、外部からそのシステムにどのように影響を与えていけばよいかを考える問題です。外部からシステムに与える影響を制御入力と呼ばれ、以下では時刻 t における制御入力を u(t) と表すことにします。

最適制御問題の単純な例として、線形システムにおける最適制御問題を見てみます。A,\ Bn\times n の行列とするとき、状態変数 x: [0,T]\to\mathbb{R}^n と、制御入力 u: [0,T]\to\mathbb{R}^n を用いて制御対象となる系が

\dot{x}(t)=A x(t)+ B u(t),\ x(0)=x_0\tag{1}

で与えられているとします。このとき以下の汎関数

J(x,u):=\int_0^T({}^tx(t)Qx(t)+{}^tu(t)Ru(t))dt \tag{2}

を最小にするために u(t) をどのような関数にすればよいかを考えることになります。最適制御問題のポイントは制約条件(1)が制御入力 u を含む状態方程式であるということです。

流体の方法による最適制御

始めに述べた最適輸送問題の定式化における最適写像 M の存在証明や最適写像を求める方法はこれまでにいくらか提案されてきました。最適写像 M を具体的に作る方法として、恒等写像から逐次的に少しづつ作りあげていく方法は制御理論との相性もよさそうで、また理解しやすいだろうと思われます。

2000年に[3]において提案された L^2-MKPを流体を利用した最適制御問題として捉える方法について概説します。そのアイデアは、質量分布を時間に依存する流体の密度分布 \rho(t,x) と見なし、その流体の流れ場 v(t,x) を適切に制御することで、その流れに乗って流体が移動することが最適輸送写像になっているようにするという考えです。この方法はDensity controlと呼ばれることがあります。

すなわち、\mathbb{R}^d上の1パラメータ変換群 \varphi:[0,T]\times\mathbb{R}^d\to\mathbb{R}^d\varphi(0,x)=x,\varphi(T,x)=M(x) となるような \varphi を求めようというわけです。
流速場を v(t,x)=\partial_t\varphi(t,x) とするとき、密度 \rho(t,x) が連続の式と境界条件
\partial_t\rho+{\rm div}(\rho v)=0 \cdots (3)
\rho(0,x)=\rho_0,\rho(T,x)=\rho_T \cdots (4)
を満たすことを要請します。また(3),(4)が満たされるときMKP(i)が満たされることが分かります。

さらに
T\int_{\mathbb{R}^d}\int^T_0\rho(t,x)|v(t,x)|^2dx^ddt\geq \int_{\mathbb{R}^d}|M(x)-x|^2\rho_0(x)dx^d
となることが少し計算すると分かるので、MKP(ii)の条件の代わりに
S(\rho,v):=T\int_{\mathbb{R}^d}\int^T_0\rho(t,x)|v(t,x)|^2dx^ddt
を最小にする問題を考えてもよいことになります。
従って L^2-MKPは(3),(4)を状態方程式、制御入力を v(t,x)、状態変数を \rho(t,x) と見なしたとき、S を最小化する最適制御問題ということになります。

(3)式は任意の t,x に対して成り立たないといけないので、t,x に依存するLagrange multiplier \phi(t,x) を使うと
S(\rho,v,\phi):=\int_{\mathbb{R}^d}\int^T_0\left(\frac{1}{2}\rho|v|^2+\phi(\partial_t\rho+{\rm div}(\rho v))\right)dx^ddt
の極値問題を解くことが必要条件であることが分かります。
\rho,v,\phi で変分して整理すると
\partial_t\phi+\frac{1}{2}|{\rm grad}\phi|^2=0 \cdots(5)
v={\rm grad}\phi \cdots(6)
が得られます。
よってL^2-MKPを解くには連立系(3),(4),(5),(6)を解き、速度場 v(t,x) に沿って T 秒移動すれば求める写像が得られることになります。
時刻 t において制御入力 v(t,x) をどのように決定していけばよいのかは技術的には拡張ラグランジュ法([4]のchap2 Application to the Stokes and Navier-Stokes equationなど)に関していくらかの準備が必要ですのでこの記事では詳細に踏み込まずアイデアの紹介までとしたいと思います。

参考文献

[1] Charlotte Bunne, "Optimal Transport in Learning, Control, and Dynamical Systems"
[2] N.Ghoussoub, Y-H. Kim, A.Z. Palmer, "OPTIMAL TRANSPORT WITH CONTROLLED DYNAMICS AND FREE END TIMES"
[3] J.D. Benamou, Y. Brenier, "A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem"
[4] Fortin, M., Glowinski, R. " Augmented Lagrangian methods. Applications to the numerical solution of boundary value problems"

Discussion