👌

Density Controlによる最適輸送

2024/07/27に公開

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

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

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

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

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

最適輸送問題

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

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

分布が連続な場合にコストを距離の２乗で定めるときの最適輸送問題は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$ で以下の２つを満たすものを求める問題です。
(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,\ B$ を $n\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$ 上の１パラメータ変換群 $\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など）に関していくらかの準備が必要ですのでこの記事では詳細に踏み込まずアイデアの紹介までとしたいと思います。

計算アルゴリズム

上記までに最適輸送を実現する流体の流れ $(\rho(t,x), v(t,x))$ は

L(\phi,\rho,v):=\int_0^Tdt\int_Ddx \{\rho(t,x)|v(t,x)|^2 + \phi(t,x)(\partial_t\rho+\textrm{div}(\rho v))\}dtdx

と置くとき、

\inf_{\rho, v} \sup_{\phi} L(\phi, \rho, v)

を実現することを論じた。この記事では具体的にこれを求めるアルゴリズムを解説する。

拡張ラグランジュ関数の導入

新たに運動量 $m=\rho v$ を定義すると、 $L(\phi,\rho,v)$ は

L(\phi, \rho, m) = \int_0^T \int_D \left( \frac{|m|^2}{2\rho} - \partial_t \phi \rho - \nabla_x \phi \cdot m \right)dxdt- \int_D \left( \phi(0, \cdot) \rho_0 - \phi(T, \cdot) \rho_T \right)dx

となるから

\inf_{\rho, m} \sup_{\phi} L(\phi, \rho, m) \tag{1}

を考えればよい。変分すると

\partial_t \phi + \frac{|m|^2}{2\rho} = 0,\ \frac{m}{\rho} = \nabla_x \phi,\ \partial_t \rho + \textrm{div}m = 0

\rho(0, \cdot) = \rho_0,\ \rho(T, \cdot) = \rho_T

が得られる。

これまで空間 $\mathbb{R}^n\ (n=2,3)$ を考えてきたが、時間方向の次元を加えて全体として $M:=\mathbb{R}\times \mathbb{R}^n$ を考え、これを拡張空間と呼ぶことにする（数学における時空とは別の概念なのでこのように呼び分けることにする）。拡張空間 $M$ のベクトル $\xi$ は時間成分が $\tau$ で空間成分( $\mathbb{R}^n$ のベクトル)が $s$ のとき、

\xi = (\tau, s)

と書くことにする。

このとき

\mu= (\rho, m)\\ q = (a,b)

と置く（ $a, b$ は $\rho, m$ の双対変数として導入する）。また拡張空間のベクトル同士の内積 $\cdot$ はEuclid内積（成分ごとの積の和）とする。さらに $L^2$ 内積を

\langle\mu,q\rangle:=\int_0^T\int_D\mu\cdot q\ dtdx

と定義する。 $D\subset\mathbb{R}^n$ は被積分関数の台である。以下では登場する場は全て $\mathbb{R}^n$ においてコンパクト台を持つことを仮定し、 $D$ は関連する全ての場の台を含むとする。

ここで

K =\{ q: \mathbb{R} \times \mathbb{R}^3 \rightarrow \mathbb{R} \times \mathbb{R}^3;\ \ a + \frac{|b|^2}{2} \leq 0 \}

と定義すると、

\frac{|m(t, x)|^2}{2\rho(t, x)} = \sup_{(a, b) \in K} [a(t, x) \rho(t, x) + b(t, x) \cdot m(t, x)]

となる。

F(q) = \begin{cases}0 \ \ (q\in K)\\ \infty \ \ (q\notin K)\end{cases}\\ G(\phi) = \int_D \left[ \phi(0, \cdot) \rho_0 - \phi(T, \cdot) \rho_T \right] \, dx

と定義すると、(1)は

\sup_\mu \inf_{\phi, q} \left[ F(q) + G(\phi) + \langle \mu, \nabla_{t,x} \phi - q \rangle \right] \tag{2}

と同値である。ここで拡張空間のgradientを

\nabla_{t,x} \phi(t,x) := (\partial_t\phi, \textrm{grad}\phi)

と書いた。

さらに（２）に対して、 $\mu$ を制限 $\nabla_{t,x} \phi(t,x)-q=0$ に関するLagrange乗数と読み替えることにより、拡張ラグランジュ関数を

L_r(\phi, q, \mu) := F(q) + G(\phi) + \langle \mu, \nabla_{t,x} \phi - q \rangle + \frac{r}{2} \langle \nabla_{t,x} \phi - q, \nabla_{t,x} \phi - q \rangle

と定義すると、

\sup_{\mu} \inf_{\phi, q} L_r(\phi, q, \mu) \tag{3}

を考えてもよいことになる。

極値問題（３）を解くにはALG2アルゴリズムを利用する。ALG2アルゴリズムは以下である。

$(\phi^{n-1}, q^{n-1}, \mu^n)$ が与えられたとする。
Step A : 以下を満たす $\phi^n$ を見つける。

L_r(\phi^n, q^{n-1}, \mu^n) \leq L_r(\phi, q^{n-1}, \mu^n), \quad \forall \phi.

Step B : 以下を満たす $q^n$ を見つける。

L_r(\phi^n, q^n, \mu^n) \leq L_r(\phi^n, q, \mu^n), \quad \forall q.

Step C : 以下のように $\mu^n$ を更新する。

\mu^{n+1} = \mu^n + r (\nabla_{t,x} \phi^n - q^n)

Step A に戻る。
以下では各stepについて詳細に論じる。

Step A

\phi ^n= \phi + \varphi, \quad \forall \phi \in C(\mathbb{R}^n), \quad L_r(\phi^n, q^{n-1}, \mu^n) \leq L_r(\phi, q^{n-1}, \mu^n)

であるとすると、 $L_r(\phi^n, q^{n-1}, \mu^n)$ を $\varphi$ で変分することで必要条件として

G(\varphi) + \langle \mu^n, \nabla_{t,x} \varphi \rangle + r \langle \nabla_{t,x} \phi^n - q^{n-1}, \nabla_{t,x} \varphi \rangle = 0 \tag{A1}

が得られる。

r \langle \nabla_{t,x} \phi^n - q^{n-1}, \nabla_{t,x} \varphi \rangle = r \langle (\partial_t \phi^n - a^{n-1}, \textrm{grad} \phi^n - b^{n-1}), (\partial_t\varphi, \textrm{grad} \varphi) \rangle \\ = r\int_0^Tdt\int_Ddx\left( \partial_t\varphi(\partial_t\phi^n-a^{n-1})+(\textrm{grad}\phi^n-b^{n-1})\cdot\textrm{grad}\varphi \right)\\ \langle \mu^n, \nabla_{t,x} \varphi \rangle = \langle (\rho^n, m^n), (\partial_t \varphi, \textrm{grad} \varphi) \rangle =\int_0^Tdf\int_Ddx( \partial_t\varphi\rho^n + m^n\cdot \textrm{grad}\varphi )

であるから、

\langle \mu^n, \nabla_{t,x} \varphi \rangle + r \langle \nabla_{t,x} \phi^n - q^{n-1}, \nabla_{t,x} \varphi \rangle\\ =\int_0^Tdt\int_Ddx\left( \partial_t\varphi (\rho^n + r(\partial_t\phi^n-a^{n-1})) + (m^n + r(\textrm{grad}\phi^n-b^{n-1}))\cdot\textrm{grad}\varphi \right)\\ =\int_0^Tdt\int_Ddx \frac{d}{dt}(\varphi (\rho^n + r(\partial_t\phi^n-a^{n-1}))) - \int_0^Tdt\int_Ddx( \varphi (\partial_t\rho^n + r\partial_t(\partial_t\phi^n-a^{n-1}))\\ + \int_0^T\int_Ddx(m^n + r(\textrm{grad}\phi^n-b^{n-1}))\cdot\textrm{grad}\varphi )\\ =\int_Ddx \{\varphi(T,x) (\rho^n(T,x) + r(\partial_t\phi^n(T,x)-a^{n-1}(T,x)))\} - \int_Ddx \{\varphi(0,x) (\rho^n(0,x) + r(\partial_t\phi^n(0,x)-a^{n-1}(0,x)))\}\\ - \int_0^Tdt\int_Ddx \varphi \{\partial_t\rho^n + r\partial_t(\partial_t\phi^n-a^{n-1}) + \textrm{div}(m^n + r(\textrm{grad}\phi^n-b^{n-1}))\}\\ + \int_0^Tdt\int_Ddx \textrm{div}(\varphi(m^n + r(\textrm{grad}\phi^n-b^{n-1})))

ここで、場がコンパクト台を持つと仮定しているから、

\langle \mu^n, \nabla_{t,x} \varphi \rangle + r \langle \nabla_{t,x} \phi^n - q^{n-1}, \nabla_{t,x} \varphi \rangle\\ =\int_Ddx \{\varphi(T,x) (\rho^n(T,x) + r(\partial_t\phi^n(T,x)-a^{n-1}(T,x)))\} - \int_Ddx \{\varphi(0,x) (\rho^n(0,x) + r(\partial_t\phi^n(0,x)-a^{n-1}(0,x)))\}\\ - \int_0^Tdt\int_Ddx \varphi \{\partial_t\rho^n + \textrm{div}(m^n) + r(\partial^2_t\phi^n + \Delta\phi^n-\partial_ta^{n-1} - \textrm{div} b^{n-1})\}

となる。従って、(A1)は

G(\varphi) + \langle \mu^n, \nabla_{t,x} \varphi \rangle + r \langle \nabla_{t,x} \phi^n - q^{n-1}, \nabla_{t,x} \varphi \rangle\\ =\int_Ddx \{\varphi(T,x) (\rho^n(T,x) -\rho(T,x) + r(\partial_t\phi^n(T,x)-a^{n-1}(T,x)))\}\\ - \int_Ddx \{\varphi(0,x) (\rho^n(0,x) -\rho(0,x) + r(\partial_t\phi^n(0,x)-a^{n-1}(0,x)))\}\\ - \int_0^Tdt\int_Ddx \varphi \{\partial_t\rho^n + \textrm{div}(m^n) + r(\partial^2_t\phi^n + \Delta\phi^n-\partial_ta^{n-1} - \textrm{div} b^{n-1})\}=0

となるから、結局

-r\Delta_{t,x}\phi^n=\nabla_{t,x}\cdot(\mu^n-rq^{n-1}) \tag{A2}\\ r\partial_t\phi^n(T,x)=\rho_T -\rho^n(T,x) + ra^{n-1}(T,x)\\ r\partial_t\phi^n(0,x)=\rho_0 -\rho^n(0,x) + ra^{n-1}(0,x)

が得られる。
$\mu^n, q^{n-1}$ が与えられているとすると、(A2)は $\phi^n$ に関するPoisson方程式である。

Step B

$\phi^n,\mu^n$ が与えられたとき、

\inf_q\{F(q)+ \frac{r}{2} \langle \nabla_{t,x} \phi^n - q, \nabla_{t,x} \phi^n - q \rangle + \langle \mu^n, \nabla_{t,x} \phi^n - q \rangle\}\\ =\frac{r}{2}\inf_{q\in K}\{ \langle \nabla_{t,x} \phi^n - q, \nabla_{t,x} \phi^n - q \rangle + \frac{2}{r}\langle \mu^n, \nabla_{t,x} \phi^n - q \rangle\}\\ =\frac{r}{2}\inf_{q\in K}\{ \langle \nabla_{t,x} \phi^n - q + \frac{1}{r}\mu^n, \nabla_{t,x} \phi^n - q + \frac{1}{r}\mu^n\rangle - \frac{1}{r^2}||\mu^n||^2\}

を解けばよい(最後の式は平方完成した)。
よって

\inf_{q\in K}\{ \langle \nabla_{t,x} \phi^n - q + \frac{1}{r}\mu^n, \nabla_{t,x} \phi^n - q + \frac{1}{r}\mu^n\rangle\}

を考えればよい。

\alpha^n(t,x):=\partial_x\phi^n(t,x) + \frac{\rho^n(t,x)}{r}\\ \beta^n(t,x):=\textrm{grad}\phi^n(t,x) + \frac{m^n(t,x)}{r}

と置くと、 $q=(a,b)$ は

\inf_{(a,b)}\{(a-\alpha^n)^2+|b-\beta^n|^2;\ a+\frac{|b|}{2}\le0\} \tag{B1}

を求めればよい。これは各 $(t,x)$ 毎に求めることができる。

step C

最後に $\mu$ を $L_r$ が大きくなるように更新するためには、

\mu^{n+1} = \mu^n + r (\nabla_{t,x} \phi^n - q^n) \tag{C1}

とすればよい。

最適輸送を実現する流れを求めるアルゴリズム

これまでの結果をまとめるとアルゴリズムは以下のようになる。

始めに $\rho_0(x),\rho_T(x)$ が与えられている。

$(\phi^{n-1}, q^{n-1}, \mu^n)$ が与えられたとする。
Step A : 以下を満たす $\phi^n$ を見つける。

-r\Delta_{t,x}\phi^n=\textrm{grad}_{t,x}(\mu^n-rq^{n-1})\\ r\frac{\partial}{\partial t}\phi^n(T,x)=\rho_T(x) -\rho^n(T,x) + ra^{n-1}(T,x)\\ r\frac{\partial}{\partial t}\phi^n(0,x)=\rho_0(x) -\rho^n(0,x) + ra^{n-1}(0,x)

Step B : 以下を満たす $q^n=(a^n,b^n)$ を見つける。

\alpha^n(t,x):=\partial_x\phi^n(t,x) + \frac{\rho^n(t,x)}{r}\\ \beta^n(t,x):=\textrm{grad}\phi^n(t,x) + \frac{m^n(t,x)}{r}\\ \inf_{(a,b)}\{(a^n-\alpha^n)^2+|b^n-\beta^n|^2;\ a+\frac{|b|}{2}\le0\}