はしくれエンジニアもどきのメモ

情報系技術・哲学・デザインなどの勉強メモ・備忘録です。

微分方程式モデルでPursuitCurve問題を解く:hawk-pigeon問題

微分方程式モデルでPursuitCurve問題を解く:hawk-pigeon問題

微分方程式モデルでPursuitCurve問題:hawk-pigeon問題(鳩を追いかける鷹の追跡曲線)を解くメモ.

hawk-pigeonモデルとここでは呼んでいるが他の専門書では

  • 商船を追いかける海賊船
  • 飛んでいる飛行機を撃ち落とすための迎撃ミサイル

などの追跡曲線で紹介されている.

今回は以下をまとめる.

  • その微分方程式がどんな形になるか
  • どう解くか
  • capture time(捕捉時間)はどうなるか

図1. hawk-pigeon problem

参考:

環境

  • Windows10 ver2004(OSビルド 19041.330)
  • Docker for win v19.03.12
    • jupyter/scipy-notebook

Pursuit Curveとは

Pursuit Curveとは日本語で追跡曲線のことである. ある追跡者(chaser)が環境の外部要因を受けながら対象(target)を追いかけるときの軌跡を求める問題である.

解析としては,初期位置や初期速度により追跡者が対象に到達可能化どうか,到達可能な場合そのときの捕捉時間(capture time)などを調べる.

対象物への追跡曲線の作り方・考え方としては,追跡曲線の接線が常に対象物と交わるにすることで その軌跡はいずれ接線と平行するようになり,十分な速度があれば対象物に到達できる.

接線は常に対象と交わる

問題設定

2次元平面$(x,y)$を考え,

  • 鳩Q(target)が等速$v>0$でこの座標$(b=0,vt)$上を飛んでいるとする.
  • 鷹P(chaser)は時間$t=0$で初期位置$(a>0, 0)$で鳩Qを見つけそれを等速$w>0$で追跡しようとする.

このときの鷹Pの追跡曲線を求める.

微分方程式を立てる

鷹の軌道$P(x_{h}, y_{h})$の接線の傾きは

$$ \frac{dy(x_{h})}{dx}= \frac{y - y_{h}}{x - x_{h}} $$

鳩Q$(b=0, vt)$と交わるので

$$ \frac{dy(x_h)}{dx} =\frac{y-y_h}{x-x_h}\\ = \frac{vt-y_h}{b-x_h}\\ = \frac{vt-y_h}{0-x_h}\\ = \frac{y_h-vt}{x_h}\\ $$

時間$t$の式に直す.


\begin{eqnarray*}
\frac{vt-y}{b-x} &=& \frac{dy}{dx}\\
vt-y &=& (b-x)\frac{dy}{dx}\\
vt &=& y + (b-x)\frac{dy}{dx}\\
\end{eqnarray*}

\begin{eqnarray}
t = \frac{y}{v} - \frac{x-b}{v}\frac{dy}{dx}, (1)
\end{eqnarray}

合成速度$w$は

$$ w = \sqrt{\dot{x}^{2} + \dot{y}^{2}}\ = \sqrt{(\frac{dx}{dt})^{2} + (\frac{dy}{dt})^{2}} $$

時間tで積分する.

$$ wt = \int_{0}^{t} \sqrt{(\frac{dx}{dt'})^{2} + (\frac{dy}{dt'})^{2}}dt' $$

時間[0, t]でx軸で[a,a-x]まで移動したとすると

$$ wt = \int_{a}^{a-x} \sqrt{(dx')^{2} + (dy)^{2}}\\ = \int_{a}^{a-x} \sqrt{1 + (\frac{dy}{dx'})^{2}}dx' $$

$a-x=x'$として変数変換する.

$$ dx' = -dx\\ wt = \int_{0}^{x'} \sqrt{1 + (\frac{dy}{dx})^{2}}(-dx)\\ = - \int_{0}^{x'} \sqrt{1 + (\frac{dy}{dx})^{2}}dx $$

時間tの式に直すと


\begin{eqnarray}
t= -\frac{1}{w} \int_0^x \sqrt{1 + (\frac{dy}{dx'})^2}dx', (2)
\end{eqnarray}

(1)式と(2)式で等式を作ると,

$$ -\frac{1}{w} \int_{0}^{x} \sqrt{1 + (\frac{dy}{dx'})^{2}}dx' = \frac{y}{v} - \frac{x-b}{v}\frac{dy}{dx} $$

非線形微分方程式になる.

置換して解ける形に直す. $p(x) = \frac{dy}{dx}$として置換する.

$$ -\frac{1}{w} \int_{0}^{x} \sqrt{1 + p(x')^{2}}dx' = \frac{y}{v} - \frac{x-b}{v}p(x)\ $$

両辺を$x$で微分する.

$$ -\frac{1}{w} (\sqrt{1 + p(x)^{2}} - \sqrt{1 + p(a)^{2}}) = \frac{1}{v}\frac{dy}{dx} - \frac{x-b}{v}\frac{dp}{dx} - \frac{1}{v}p(x) $$

$dp/dx$の式に直す.


\begin{eqnarray*}
\frac{x-b}{v}\frac{dp}{dx} &=& \frac{1}{w} (\sqrt{1 + p(x)^2})\\
(x-b)\frac{dp}{dx} &=& \frac{v}{w} (\sqrt{1 + p(x)^2})\\
(x-b)\frac{dp}{dx} &=& u (\sqrt{1 + p(x)^2}), u=\frac{v}{w}\\
\end{eqnarray*}

微分方程式を解く

変数分離可能な形になったので解くことができる.

変数分離に直す.

$$ \frac{dp}{\sqrt{1 + p(x)^{2}}} = \frac{u}{x-b} dx\ $$

$x>b$として積分する.

$$ \ln{[p + \sqrt{1 + p(x)^{2}}]} + C = u\ln{(x-b)}\ $$

初期条件$x=a, y=0, p(a)=dy/dx = 0, b=0 $として定数を求める.


\begin{eqnarray*}
\ln{[p + \sqrt{1 + p(x)^2}]} + C &=& u\ln{(x-b)}\\\\
\ln{[0 + \sqrt{1 + 0^2}]} + C &=& u\ln{(x-b)}\\\\
C &=& u\ln{(a-b)}, (a>b) 
\end{eqnarray*}

定数を代入して


\begin{eqnarray*}
\ln{[p + \sqrt{1 + p(x)^2}]} + u\ln{(a-b)} &=& u\ln{(x-b)}\\\\
\ln{[p + \sqrt{1 + p(x)^2}]} - \ln{(x-b)^{u}}+\ln{(a-b)^u} &=& 0
\end{eqnarray*}

変形する


\begin{eqnarray*}
0 = \ln{[p + \sqrt{1 + p(x)^2}]} - \ln{(x-b)^u}+\ln{(a-b)^u}\\
0 = \ln{[p + \sqrt{1 + p(x)^2}]} + \ln{\left(\frac{a-b}{x-b}\right)^u}
\end{eqnarray*}

$\ln{1}=0$より

$$ [p + \sqrt{1 + p(x)^{2}}] \cdot \left(\frac{a-b}{x-b}\right)^{u} = 1 $$

左辺をpの式に直す.

$$ [p + \sqrt{1 + p(x)^{2}}] = \left(\frac{a-b}{x-b}\right)^{-u}\ $$

$q = \left(\frac{a-b}{x-b}\right)^{-u}$とおいて


\begin{eqnarray*}
p + \sqrt{1+p^2} &=& q \\
\sqrt{1+p^2} &=& q - p \\
1+p^2 &=& (q-p)^2\\
1+p^2 &=& q^2 - 2pq + p^2\\
2pq &=& q^2 -1 \\
\end{eqnarray*}

$$ p = \frac{1}{2}\left( q - \frac{1}{q} \right) $$

pとqの変数をもとに戻すと

$$ \frac{dy}{dx} = \frac{1}{2}\left( \left(\frac{a-b}{x-b}\right)^{-u} - \left(\frac{a-b}{x-b}\right)^{u} \right) $$

積分すると


\begin{eqnarray*}
y(x)
&=& \frac{1}{2}\left( \int \left(\frac{a-b}{x-b}\right)^{-u}dx - \int\left(\frac{a-b}{x-b}\right)^{u} dx \right) +C\\
&=& \frac{1}{2}\left( \left( \frac{(x-b)^{u+1}}{(u+1)(a-b)^{u}} \right) - \left(\frac{(a-b)^{u}(x-b)^{-u+1}}{(-u+1)}\right) \right) +C
\end{eqnarray*}

初期条件$x=a, y(a)=0$より,定数を求める.


\begin{eqnarray*}
0 &=& \frac{1}{2}\left( \left( \frac{(a-b)^{u+1}}{(u+1)(a-b)^{u}} \right) - \left(\frac{(a-b)^{u}(a-b)^{-u+1}}{(-u+1)}\right) \right)+C\\
C &=& \frac{1}{2}\left( \left( \frac{(a-b)}{(-u+1)} \right) - \left(\frac{(a-b)}{(u+1)}\right) \right)\\
&=& \frac{1}{2}\left( \frac{(a-b)(1+u) - (a-b)(1-u)}{1-u^2} \right)\\
&=& \frac{1}{2}\left( \frac{2u(a-b)}{1-u^2} \right)\\
&=& \left( \frac{u(a-b)}{1-u^2} \right)\\
\end{eqnarray*}

まとめると,


\begin{eqnarray}
y(x)
&=& \frac{1}{2}\left( \left( \frac{(x-b)^{u+1}}{(u+1)(a-b)^{u}} \right) - \left(\frac{(a-b)^{u}(x-b)^{-u+1}}{(-u+1)}\right) \right) +C\\
C&=& \left( \frac{u(a-b)}{1-u^2} \right)\\
\end{eqnarray}

鳩Qの初期x座標が$b=0$とすると


\begin{eqnarray}
y(x)
&=& \frac{1}{2}\left( \left( \frac{x^{u+1}}{(u+1)a^{u}} \right) - \left(\frac{a^{u}x^{-u+1}}{(-u+1)}\right) \right) +C\\
C&=& \left( \frac{ua}{1-u^2} \right)\\
\end{eqnarray}

可視化と分析

鳩と鷹の速度$v, w$を変えて導出した式を使って可視化する.


\begin{eqnarray}
y(x)
&=& \frac{1}{2}\left( \left( \frac{x^{u+1}}{(u+1)a^{u}} \right) - \left(\frac{a^{u}x^{-u+1}・・}{(-u+1)}\right) \right) +C\\
C&=& \left( \frac{ua}{1-u^2} \right)\\
\end{eqnarray}
import numpy as np
import matplotlib.pyplot as plt

v = 1
w = 1.5
u = v/w
a = 1
C = ( u * a / (1-u**2) )
# print(C)

x = np.linspace(0, a, 1000)
y = .5 * ((x**(u+1)/((u+1)*a**(u)) ) - ((x**(-u+1)*(a)**u)/((-u+1)))) + C
cap_t = a/(w*(1-u**2))
print(f"cap_t = {cap_t}")

fig = plt.figure(figsize=(5,5))
plt.plot(x, y, label=f"$v/w={u:.3}$")

v = 1
w = 0.999
u = v/w
C = ( u * a / (1-u**2) )
y = .5 * ((x**(u+1)/((u+1)*a**(u)) ) - ((x**(-u+1)*(a)**u)/((-u+1)))) + C
plt.plot(x, y, label=f"$v/w={u:.3}$")
cap_t = a/(w*(1-u**2))
print(f"cap_t = {cap_t}")

v = 1.5
w = 1
u = v/w
C = ( u * a / (1-u**2) )
y = .5 * ((x**(u+1)/((u+1)*a**(u)) ) - ((x**(-u+1)*(a)**u)/((-u+1)))) + C
plt.plot(x, y, label=f"$v/w={u:.3}$")
cap_t = a/(w*(1-u**2))
print(f"cap_t = {cap_t}")

plt.xlim(-.1, 2)
plt.ylim(-.1, 2)
plt.grid(True)
plt.legend()
plt.show()

図2. 可視化

速度比3パターンを比較する.

  • $u < w$: 鳩の横軸$b=0$に近づくのが早く鳩に到達する.
  • $u=w$: yの式から計算不可能,つまり鳩に追いつかない.
  • $u>w$: 計算は可能だが鳩の横軸$b=0$に近づくのが遅く,横軸にたどり着いたとしても垂直直線y軸から速度差により追いつけない.

capture time(捕捉時間)を求める

capture timeを求める.

$x=b, y=vt$


\begin{eqnarray*}
y(x=b) &=& 0 +C\\\\
vt &=& 0 + C\\\\
vt &=& C\\\\
t &=& \frac{C}{v}\\\\
&=& \frac{1}{v}\frac{u(a-b)}{1-u^2}\\\\
&=& \frac{(a-b)}{w(1-u^2)}
\end{eqnarray*}

$b=0$とすると


\begin{eqnarray}
t &=& \frac{(a-b)}{w(1-u^2)}
\end{eqnarray}
  • $1-u^{2} > 0 \to 1 > u^{2} \to w^{2} > v^{2} \to w > v$のとき, 計算可能.つまり鳩に追いつく.

  • $1-u^{2} \leq 0 \to 1 \leq u^{2} \to w^{2} \leq v^{2} \to w \leq v$のときcapture timeはマイナスになり計算不可,つまり鷲は鳩に追いつけない.