随伴モデルを
作ってみよう

榎本剛

Lotka-Volterraモデル

捕食・被食モデル
化学反応 (Lotka 1920)や魚種交替 (Volterra 1926)モデル
海洋生態系（NPZ）モデル (Franks et al. 1986)
積雲対流の自己組織化のモデル (Nober and Graf 2005)

支配方程式系

二元連立非線型常微分方程式系

\[ \begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}t} &= x(a_1 + a_2x + a_3y) \\ \frac{\mathrm{d}y}{\mathrm{d}t} &= y(a_4 + a_5y + a_6x) \end{aligned} \]

変数の定義

変数	定義	単位	標準値
\(x,\,y\)	密度	数\(\,\mathrm{m}^{-2}\)	計算
\(a_1,\,a_4\)	比増加率	\(\mathrm{d}^{-1}\)	\(4,-6\)
\(a_2,\,a_5\)	自己依存増加率	\((\)数\(\,\mathrm{m^{-2}})^{-1}\mathrm{d}^{-1}\)	\(-2, 2\)
\(a_3,\,a_6\)	相手依存増加率	\((\)数\(\,\mathrm{m^{-2}})^{-1}\mathrm{d}^{-1}\)	\(-4, 4\)
\(x_1,\, y_1\)	初期密度	数\(\,\mathrm{m}^{-2}\)	\(1, 1\)

離散化

\[ \begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}t} &= x(a_1 + a_2x + a_3y) \\ \frac{\mathrm{d}y}{\mathrm{d}t} &= y(a_4 + a_5y + a_6x) \end{aligned} \]

Fortran
Python

do n = 1, nmax-1
  x(n+1) = x(n) + dt * (x(n) * (a(1) + a(2) * x(n) + a(3) * y(n)))
  y(n+1) = y(n) + dt * (y(n) * (a(4) + a(5) * y(n) + a(6) * x(n)))
end do

for i in range(nmax-1):
    x[n+1] = x[n] + dt * (x[n] * (a[0] + a[1] * x[n] + a[2] * y[n]))
    y[n+1] = y[n] + dt * (y[n] * (a[3] + a[4] * y[n] + a[5] * x[n]))

時間発展

15日間の時間発展
dt = 0.001, nmax = 15001

位相平面

おさらい

変分法データ同化では入力（初期値やパラメタ）を推定
出力（予測値）と観測との乖離（コスト函数）の最小化
コスト函数の入力についての勾配を用いて最適化
随伴を求めてから離散化か、離散化してから随伴か。
直接勾配が計算できないので連鎖律を利用する。
Lagrangeの未定乗数を用いる方法 (Lawson et al. 1995)

Lagrangeの未定乗数法

モデル\(F(X, Z, ...)\)、制御変数 \(X\)、状態変数 \(Z\)
ステップ \(n\) での \(\lambda_Z\) は \(Z\) に関するモデルの微分と
ステップ \(n+1\) での \(\lambda_Z\)との積
観測がある時刻はコスト函数の \(Z\) 微分を加える。
\(\lambda_X\)は \(\lambda_Z\) と \(X\) に関するモデルの微分の積の総和

\[ \lambda_X = \lambda_X + \lambda_Z\frac{\partial F}{\partial X} \]

制御変数と状態変数、随伴変数

制御変数は \(X \equiv (x_1,\,y_1,\,a_1,\,a_2,\,a_3,\,a_4,\,a_5,\,a_6)^\mathrm{T}\)
状態変数は \(Z \equiv (x, y)^\mathrm{T}\)
パラメタ数の長さ6の1次元配列 a
ステップ数の長さnmaxの1次元配列 x
a の随伴変数は a をつけた長さ6の1次元配列 aa
x と y の随伴変数は、それぞれ ax と ay

随伴モデル作成のルール

時間を逆行するので、コードを逆順に辿る。
ループは逆に回す。
順行コードの右辺を制御変数 a6 で微分する。
微分に状態変数の随伴 ay を掛ける。
これを制御変数の随伴 a6 に足し込む。
制御変数 y が複数回出てきたら、それぞれ微分する。
全ての制御変数について行う。
随伴変数は0に初期化する。

捕食者パラメタの随伴

\[ \frac{\mathrm{d}y}{\mathrm{d}t} = y(a_4 + a_5y + a_6x) \]

Fortran
Python

! y(n+1) = y(n) + dt * (y(n) * (a(4) + a(5) * y(n) + a(6) * x(n)))
aa(6) = aa(6) + dt * x(n) * y(n) * ay(n+1)
aa(5) = aa(5) + dt * y(n) * y(n) * ay(n+1)
aa(4) = aa(4) + dt * y(n) * ay(n+1)

# y[n+1] = y[n] + dt * (y[n] * (a[3] + a[4] * y[n] + a[5] * x[n]))
aa[5] = aa[5] + dt * x[n] * y[n] * ay[n+1]
aa[4] = aa[4] + dt * y[n] * y[n] * ay[n+1]
aa[3] = aa[3] + dt * y[n] * ay[n+1]

状態変数の随伴

\[ \frac{\mathrm{d}y}{\mathrm{d}t} = y(a_4 + a_5y + a_6x) \]

Fortran
Python

! y(n+1) = y(n) + dt * (y(n) * (a(4) + a(5) * y(n) + a(6) * x(n)))
ax(n) = ax(n) + dt * a(6) * y(n) * ay(n+1)
ay(n) = ay(n) + dt * a(5) * y(n) * ay(n+1)
ay(n) = ay(n) + (1 + dt * (a(4) + a(5) * y(n) + a(6) * x(n))) * ay(n+1)

# y[n+1] = y[n] + dt * (y[n] * (a[3] + a[4] * y[n] + a[5] * x[n]))
ax[n] <- ax[n] + dt * a[5] * y[n] * ay[n+1]
ay[n] <- ay[n] + dt * a[4] * y[n] * ay[n+1]
ay[n] <- ay[n] + (1 + dt * (a[3] + a[4] * y[n] + a[5] * x[n])) * ay[n+1]

被食者パラメタの随伴

\[ \frac{\mathrm{d}x}{\mathrm{d}t} = x(a_1 + a_2x + a_3y) \]

Fortran
Python

! x(n+1) = x(n) + dt * (x(n) * (a(1) + a(2) * x(n) + a(3) * y(n)))
aa(3) = aa(3) + dt * y(n) * x(n) * ax(n+1)
aa(2) = aa(2) + dt * x(n) * x(n) * ax(n+1)
aa(1) = aa(1) + dt * x(n) * ax(n+1)

# x[n+1] = x[n] + dt * (x[n] * (a[0] + a[1] * x[n] + a[2] * y[n]))
aa[2] = aa[2] + dt * y[n] * x[n] * ax[n+1]
aa[1] = aa[1] + dt * x[n] * x[n] * ax[n+1]
aa[0] = aa[0] + dt * x[n] * ax[n+1]

状態変数の随伴

\[ \frac{\mathrm{d}x}{\mathrm{d}t} = x(a_1 + a_2x + a_3y) \]

Fortran
Python

! x(n+1) = x(n) + dt * (x(n) * (a(1) + a(2) * x(n) + a(3) * y(n)))
ay(n) = ay(n) + dt * a(3) * x(n) * ax(n+1)
ax(n) = ax(n) + dt * a(2) * x(n) * ax(n+1)
ax(n) = ax(n) + (1 + dt * (a(1) + a(2) * x(n) + a(3) * y(n))) * ax(n+1)

# x[n+1] = x[n] + dt * (x[n] * (a[0] + a[1] * x[n] + a[2] * y[n]))
ay[n] = ay[n] + dt * a[2] * x[n] * ax[n+1]
ax[n] = ax[n] + dt * a[1] * x[n] * ax[n+1]
ax[n] = ax[n] + (1 + dt * (a[0] + a[1] * x[n] + a[2] * y[n])) * ax[n+1]

aa(:) = 0d0
ax(:) = 0d0
ay(:) = 0d0
do nmax-1, 1, -1
  aa(6) = aa(6) + dt * x(n) * y(n) * ay(n+1)
  aa(5) = aa(5) + dt * y(n) * y(n) * ay(n+1)
  aa(4) = aa(4) + dt * y(n) * ay(n+1)
  ax(n) = ax(n) + dt * a(6) * y(n) * ay(n+1)
  ay(n) = ay(n) + dt * a(5) * y(n) * ay(n+1)
  ay(n) = ay(n) + (1 + dt * (a(4) + a(5) * y(n) + a(6) * x(n))) * ay(n+1)
  aa(3) = aa(3) + dt * y(n) * x(n) * ax(n+1)
  aa(2) = aa(2) + dt * x(n) * x(n) * ax(n+1)
  aa(1) = aa(1) + dt * x(n) * ax(n+1)
  ay(n) = ay(n) + dt * a(3) * x(n) * ax(n+1)
  ax(n) = ax(n) + dt * a(2) * x(n) * ax(n+1)
  ax(n) = ax(n) + (1 + dt * (a(1) + a(2) * x(n) + a(3) * y(n))) * ax(n+1)
end do

aa = np.zeros(6)
ax = np.zeros(nmax)
ay = np.zeros(nmax)
for i in reversed(range(nmax-1)):
    aa[5] = aa[5] + dt * x[n] * y[n] * ay[n+1]
    aa[4] = aa[4] + dt * y[n] * y[n] * ay[n+1]
    aa[3] = aa[3] + dt * y[n] * ay[n+1]
    ax[n] = ax[n] + dt * a[5] * y[n] * ay[n+1]
    ay[n] = ay[n] + dt * a[4] * y[n] * ay[n+1]
    ay[n] = ay[n] + (1 + dt * (a[3] + a[4] * y[n] + a[5] * x[n])) * ay[n+1]
    aa[2] = aa[2] + dt * y[n] * x[n] * ax[n+1]
    aa[1] = aa[1] + dt * x[n] * x[n] * ax[n+1]
    aa[0] = aa[0] + dt * x[n] * ax[n+1]
    ay[n] = ay[n] + dt * a[2] * x[n] * ax[n+1]
    ax[n] = ax[n] + dt * a[1] * x[n] * ax[n+1]
    ax[n] = ax[n] + (1 + dt * (a[0] + a[1] * x[n] + a[2] * y[n])) * ax[n+1]

数値最適化

得られた aa や ax、ayはコスト函数の初期勾配
最急降下法 \[ X^{(i+1)} = X^{(i)} - \alpha\nabla_{X(0)}J \]
勾配を用いる共軛勾配法や準ニュートン法
scipy.optimize
Nocedal L-BFGS CG+
NLopt

課題

捕食・被食モデルの時間発展（スライド6）を調べる。
- 初期状態やパラメタを変えてみよう。
- 位相空間にプロットする（スライド7）。
随伴モデルを作成する。
- 異なる初期値とから真の初期値を復元する。
- 異なる初期値とパラメタから真値を復元する。
- 最適化手法を変えてみる。
- 背景誤差や観測誤差を考慮する。

付録

質問と回答

Q: 随伴変数を0に初期化したら，最後まで0
A: 観測時刻にはコスト函数 \(J\) の勾配が加わる．

未定乗数法の詳細

\[ \begin{aligned} L&(\mathbf{x},\,z_1, \dots, z_N,\,\lambda_1,\dots,\lambda_N) = J(\mathbf{x},\,z_{1},\dots,z_N)\\&-\lambda_1(z_1 - f_1(\mathbf{x})) - \sum_{n=2}^{N} \lambda_n(z_n - f_n(\mathbf{x}, z_{n-1})) \end{aligned} \]

制御変数 \(\mathbf{x} = (x_1, \dots, x_m)^\mathrm{T}\)
時間ステップ \(n\) での状態変数の一つ \(z_n,\,n = 1, \dots, N\)
観測 \(\mathbf{y} = (y_1, \dots, y_k)^\mathrm{T}\)
コスト函数の例 \(J = \sum_{i=1}^k (y_k - h(z_n(t_n=t_k)))^2\)

鞍点は最適値

\(\partial L/\partial\boldsymbol\lambda\): 制約に用いる順行モデル \(z_n = f_n(\mathbf{x}, z_{n-1})\)
\(\partial L/\partial\mathbf{z}\): 随伴モデル \[\lambda_N = \frac{\partial J}{\partial Z_N},\,\lambda_n = \frac{\partial J}{\partial z_n} + \frac{\partial f_{n+1}}{\partial z_n}\lambda_{n+1}\]
\(\partial L/\partial\mathbf{x}\): \(J\)の制御変数についての勾配

\[ \frac{\partial J}{\partial x_k} = \sum_{i=1}^N \frac{\partial f_i}{\partial x_k}\lambda_i,\,k = 1,\dots, m \]

参考文献

Franks, P. J. S., J. S. Wroblewski, and G. R. Flierl, 1986: Behavior of a simple plankton model with food-level acclimation by herbivores. Marine Biol., 91, 121–129, https://doi.org/10.1007/BF00397577.

Lawson, L. M., Y. H. Spitz, E. E. Hofmann, and R. B. Long, 1995: A data assimilation technique applied to a predator-prey model. Bull. Math. Biol., 57, 593–617, https://doi.org/10.1016/S0092-8240(05)80759-1.

Lotka, A. J., 1920: Analytical note on certain rhythmic relations in organic systems. PNAS, 6, 410–415, https://doi.org/10.1073/pnas.6.7.410.

Nober, F. J., and H. F. Graf, 2005: A new convective cloud field model based on principles of self-organisation. Atmos. Chem. Phys., 5, 2749–2759, https://doi.org/10.5194/acp-5-2749-2005.

Volterra, V., 1926: Fluctuations in the abundance of a species considered mathematically. Nature, 118, 558–560, https://doi.org/10.1038/118558a0.

随伴モデルを作ってみよう

Lotka-Volterraモデル

支配方程式系

変数の定義

離散化

時間発展

位相平面

おさらい

Lagrangeの未定乗数法

制御変数と状態変数、随伴変数

随伴モデル作成のルール

捕食者パラメタの随伴

状態変数の随伴

被食者パラメタの随伴

状態変数の随伴

随伴モデル

数値最適化

課題

付録

質問と回答

未定乗数法の詳細

鞍点は最適値

参考文献

随伴モデルを
作ってみよう