🌒

Julia言語の常微分方程式ソルバーを使って円制限三体問題の軌道計算をしてみる〜軌道伝播・境界値問題編〜

2023/01/06に公開

Julia

1. はじめに

本記事では，宇宙軌道力学分野でよく登場する円制限三体問題を題材にして，常微分方程式のソルバーの使い方を紹介する．Julia言語が最も活躍するシチュエーションは「常微分方程式 (Ordinary Differential Equation; ODE) を解く」というような数値計算である．プログラミング言語ごとの常微分方程式ソルバーを比較したPDFを見てもらうと，Julia言語がいかに強力か感じ取れるだろう．

Julia言語を用いて，常微分方程式の数値解法を行う際に，最も広く利用されているライブラリはDifferentialEquations.jlである．その他にもOrdinaryDiffEq.jl，ODE.jl，Sundials.jl等のライブラリがある．ODE.jlはDifferentialEquations.jl等に取って代わられた旧式のライブラリであり，本記事を投稿した時点では，使用することが非推奨とされている．また，OrdinaryDiffEq.jlおよびSundials.jlは，DifferentialEquations.jlの中に組み込まれている（一部の機能を利用するためには，個別にパッケージをインストールする必要がありそう）ため，DifferentialEquations.jlが最も一般的なライブラリとなっている．以下では，DifferentialEquations.jlを利用した解析を行う．DifferentailEquations.jlのライブラリはMITライセンスなので，ライセンスを気にせず，かなり緩く使えそうである．

Juliaでの軌道の描画には，Makie.jlを用いる．描画ライブラリで何を使うかは個人の好みだが，Makie.jlは，Plots.jlのように扱いやすく，PlotlyJS.jlのようにインタラクティブに動かせ，かつ，美しい描画ができるので，個人的なお気に入りである．

本記事で紹介するJulia言語と比較用のPythonスクリプトは，以下のGithubリポジトリのworkディレクトリ内に掲載している．

https://github.com/naoyaozaki/CRTBP.jl/

Juliaプロジェクト化しているため，以下の記事を参考に，Githubからクローン後にCRTBPのプロジェクト環境を読み込めば，必要なパッケージを一括でセットアップすることができる．

2. 円制限三体問題(CRTBP)の運動方程式

本記事では，宇宙軌道力学分野でよく登場する円制限三体問題の運動方程式を扱う．ここでは，円制限三体問題に関する詳しい説明は割愛するが，Juliaの常微分方程式ライブラリにのみ興味のある読者は本章は軽く読み飛ばしても構わない．

お互いに万有引力が作用する３つの天体がどのような運動をするのかを問う問題を三体問題（Three-Body Problem) と呼ぶ．三体問題は，その軌道を与える一般解が求積法では求まらない問題として知られており（＝要は，ある時刻に対する位置・速度を解析的な式で計算できない），宇宙軌道力学の分野でも数値解法が用いられる．三体問題の中でも，第三天体が及ぼす万有引力を無視した問題を制限三体問題 (Restricted Three-Body Problem) と呼び，更に，第一天体と第二天体がその共通重心の周りを円軌道で周回する問題を円制限三体問題 (Circular Restricted Three-Body Problem; CRTBP) と呼ぶ．宇宙軌道力学分野で対象となる「太陽＝地球＝宇宙機の三体問題」や「地球＝月＝宇宙機の三体問題」は，CRTBPで非常によく近似できる．よく近似できるとは言え，あくまで近似であり，実際の宇宙機運用を考えると無視できない誤差が生じる．それにも関わらず，なぜ軌道力学の専門家が，CRTBPを使いたがるかと言うと，CRTBPの運動方程式は，自励系（＝独立変数を陽に含まない常微分方程式）となっているからである．自励系で表されることで，力学系理論等の解析的なテクニックを利用することが可能となる．

CRTBPで用いる座標系

CRTBPでは，第一天体と第二天体の共通重心を中心として，第一天体から見た第二天体の位置ベクトル方向を $x$ 軸，第一天体に対する第二天体の角運動量ベクトル方向を $z$ 軸，右手系をなすように $y$ 軸（＝第二天体の速度ベクトル方向）を取った座標系を用いて，運動方程式を記述する．詳しい導出は省略するが，宇宙機の位置ベクトルを $\boldsymbol{r}=\begin{bmatrix}x & y & z\end{bmatrix}^{\top}$ ，速度ベクトルを $\boldsymbol{v}=\begin{bmatrix}v_x & v_y & v_z\end{bmatrix}^{\top}$ ，運動方程式を規定するパラメータを $\mu=\frac{m_2}{m_1+m_2}$ としたとき，CRTBPにおける宇宙機の運動方程式は，以下のような式で計算できる

\left\{ \begin{align*} \dot{\boldsymbol{r}} &= \boldsymbol{v}\\ \dot{\boldsymbol{v}} &= -\frac{1-\mu}{\|\boldsymbol{r}_1\|^3}\boldsymbol{r}_1 -\frac{\mu}{\|\boldsymbol{r}_2\|^3}\boldsymbol{r}_2 + \boldsymbol{a}_{\textrm{coli}} + \boldsymbol{a}_{\textrm{cf}} \end{align*} \right.

ただし，

\begin{align*} \boldsymbol{r}_1 &= \begin{bmatrix}x+\mu & y & z\end{bmatrix}^{\top}\\ \boldsymbol{r}_2 &= \begin{bmatrix}x-(1-\mu) & y & z\end{bmatrix}^{\top}\\ \boldsymbol{a}_{\textrm{coli}} &= \begin{bmatrix}2v_y & -2v_x & 0\end{bmatrix}^{\top}\\ \boldsymbol{a}_{\textrm{cf}} &= \begin{bmatrix}x & y & 0\end{bmatrix}^{\top} \end{align*}

である． $m_1$ と $m_2$ は，それぞれ第一天体と第二天体の質量を表すが，太陽＝地球系のCRTBPでは $\mu=3.036\times 10^{-6}$ ，地球＝月系のCRTBPでは $\mu=1.215\times 10^{-2}$ である．なお， $\boldsymbol{a}_{\textrm{coli}}$ はコリオリ力， $\boldsymbol{a}_{\textrm{cf}}$ は遠心力に起因する加速度である．この運動方程式は，第一天体と第二天体の距離が1，第一天体と第二天体の公転軌道周期が $2\pi$ となるように無次元化されている．

3. 円制限三体問題の軌道伝播

まず円制限三体問題の運動方程式をJuliaの関数として記述しよう．Julia言語の配列は1オリジン（1-based numbering） な定義になっている．MATLABやFORTRANユーザーは1オリジンの書き方に慣れていると思うが，PythonやC言語ユーザーは0オリジン(zero-based numbering) の書き方に慣れているため，注意が必要である．また，Julia言語では， $\mu$ 等のギリシャ文字が使えたり，数字×変数のときの掛け算を省略できたりする．関数を定義する際にはParameterizedFunctions.jlを使って，分かりやすく記述する方法もあるようだが，計算速度の観点では不利に働くため，ここではParameterizedFunctionsを利用しない方法を紹介する．

using LinearAlgebra

# Equation of Motion of CRTBP
function eom_crtbp!(dx,x,μ,t) #関数内で引数`dx`を変更するので，慣例にならって`!`を関数名の後ろに着けている
    # Preparation
    r1 = x[1:3] + [μ,0,0]
    r2 = x[1:3] - [1-μ,0,0]
    a_cori = [2x[5],-2x[4],0]
    a_cf = [x[1],x[2],0]
    # dx/dt
    dx[1:3] = x[4:6]
    dx[4:6] = -(1-μ)*r1/(norm(r1)^3) - μ*r2/(norm(r2)^3) + a_cori + a_cf
end

早速，DifferentialEquations.jlで，月軌道ゲートウェイが建設される予定の地球・月系のNear-Rectilinear Halo Orbit(NRHO)の初期条件を与えて，軌道伝播をしてみよう．ここでは，Pythonコードと比較するために，DP8()という手法を用いて，常微分方程式を解いている．公式ドキュメント曰く，同程度の精度を目指す場合はVern7()の方が性能が良いらしい．

using DifferentialEquations
using GLMakie
Makie.inline!(true)

# 問題定義 (地球・月系のNear-Rectilinear Halo Orbitを与える初期条件)
x0 = [0.987384153663276, 0.0, 0.008372273063008, 0.0, 1.67419265037912, 0.0]
tspan = (0.0, 0.5*pi)
μ = 0.01215058426994 # 地球・月系のμ
prob = ODEProblem(eom_crtbp!, x0, tspan, μ)

# 計算
@time sol = solve(prob, DP8(), reltol=1e-10, abstol=1e-10)

# 描画
fig = Figure()
ax = Axis3(fig[1, 1]; aspect=:data, perspectiveness=0.0, title="CRTBP Example 1 (NRHO)", xlabel=L"x", ylabel=L"y", zlabel=L"z")
scatterlines!(ax, 1 - μ, 0, 0; markersize=20, color=:blue)
lines!(ax, sol[1, :], sol[2, :], sol[3, :], linewidth=1, color=:red)
fig

Near-Rectilinear Halo Orbitの軌道伝播結果

本プログラムの実行時間は5.030ms(ミリ秒)であった（初回は時間が掛かるので注意）．同様の計算をPythonのscipy.integrate.solve_ivpの'DOP853'を用いて計算してみた結果，その実行時間は32msとなった．Julia言語の方がPythonよりも，6倍ほど計算時間が早いことが分かる．上記の例では，@timeで時刻計測を行なっているが，CPU時間で評価する@CPUtimeでも同程度の差がある．

4. `DifferentialEquations.jl`の拡張機能〜境界値問題の解法〜

DifferentialEquations.jlには，便利な拡張機能がある．その１つの（常微分方程式の）境界値問題(Boundary Value Problem)のソルバーを用いれば，ある境界条件を満たすような初期条件を計算することが可能である．宇宙軌道力学の分野でも，ランベール問題(Lambert's Problem)をはじめとして，二点境界値問題(Two-Point Boundary Value Problem; TPBVP) がよく登場する．この拡張機能を利用するためには，

(CRTBP) pkg> add BoundaryValueDiffEq

によって，BoundaryValueDiffEq.jlパッケージを導入する必要がある．本記事を執筆した段階では，BoundaryValueDiffEqは軌道伝播時間（tspan）を可変にした境界値問題を解くことができない．個人的にはBoundaryValueDiffEqを使わずに自分でシューティング法を実装するアプローチが好みだが，そのアプローチはここでは紹介しない（必要があれば別の記事にまとめようと思う）．

「軌道伝播時間を可変にした境界値問題を解くことができない」と言ったが，ちょっとしたテクニックを使えば，BoundaryValueDiffEqでも軌道伝播時間を可変にした問題を解くことが可能となる．例えば， $\frac{dt}{d\tau}=T$ (定数)となるように時間微分をスケーリングすることを考える（ $T$ は定数であるため， $T$ の時間微分はゼロである）．この $T$ を含めて，状態量を $[\boldsymbol{r}, \boldsymbol{v}]^{\top}$ から $[\boldsymbol{r}, \boldsymbol{v}, T]^{\top}$ に拡張した運動方程式を扱うことで，軌道伝播時間を可変にした境界値問題を解くことが可能となる．

\left\{ \begin{align*} \frac{d\boldsymbol{r}}{d\tau} &= T \boldsymbol{v}\\ \frac{d\boldsymbol{v}}{d\tau} &= T\left(-\frac{1-\mu}{\|\boldsymbol{r}_1\|^3}\boldsymbol{r}_1 -\frac{\mu}{\|\boldsymbol{r}_2\|^3}\boldsymbol{r}_2 + \boldsymbol{a}_{\textrm{coli}} + \boldsymbol{a}_{\textrm{cf}}\right)\\ \frac{dT}{d\tau} &= 0 \end{align*} \right.

さて，この境界値問題のソルバーを用いて，ハロー軌道 (Halo Orbit) というCRTBPの周期軌道を生成してみよう $^1$ ．CRTBP上の軌道は， $x$ - $z$ 平面に対する対称性を有する．したがって， $x$ - $z$ 平面を通過する周期軌道であるハロー軌道は，以下のような境界条件（＝ $x$ - $z$ 平面に対する垂直交差条件）を満たす必要がある．

\left\{ \begin{align*} y(t_0) &= 0\\ v_x(t_0) & = 0\\ v_z(t_0) & = 0\\ y(t_f) &= 0\\ v_x(t_f) & = 0\\ v_z(t_f) & = 0\\ \end{align*} \right.

なお， $t_f$ がハロー軌道半周後の時刻となるように設定し，半周分の軌道に関する境界値問題を解くことが通例である（残りの半周分はの境界値問題は，CRTBPの対称性から満たされることが保証されている）．状態量が7自由度あり，境界条件は6自由度の拘束となっているため，もう１自由度拘束しなければ，境界値問題が解けない．様々な拘束方法が考えられるが，ここでは， $x$ - $z$ 平面上の初期位置の $z$ 方向成分を固定した境界条件

z(t_0) = \bar{z}_0

を追加することを考える．他にも $x$ 方向成分を固定したり，CRTBPの保存量であるヤコビ定数(Jacobi Constant) を固定したり，遷移時間を固定するような境界条件を置くこともある．

では，BoundaryValueDiffEqライブラリのシューティング法（Shooting()）を用いて，ハロー軌道を生成してみよう．BoundaryValueDiffEqライブラリのBVProblemという構造体で問題を定義して，solve()関数で問題を解く流れで問題が解かれる．

bvp = BVProblem(f,bc!,u0,tspan,p=NullParameters();kwargs...)

ここで，f(du,u,p,t)はODEの関数(今回の場合はCRTBPの運動方程式eom_crtbp_w_time!(dx,x,p,t))，bc!(residual, u, p, t)は境界条件を規定する関数，u0は初期予想値，tspanはODEの計算時間，pはf(du,u,p,t)やbc!(residual, u, p, t)に渡すパラメータである．bc!(residual, u, p, t)は，u, p, tを入力としてredisualの配列がゼロとなるように境界条件である．

あとはsol = solve(bvp, Shooting(Vern7()))と実行すれば，境界値問題を解いてくれる．solve()関数の第一引数は境界値問題の構造体bvpで，第二引数は境界値問題のソルバーである．3種類くらいしかないが，ここではODEソルバーとしてVern7()を使ったシューティング法Shooting(Vern7()))を使う．

実際にハロー軌道を生成するコードを以下に示す．

using BoundaryValueDiffEq
using DifferentialEquations
using LinearAlgebra
using GLMakie
GLMakie.activate!()
Makie.inline!(true)

# CRTBPの運動方程式の定義（第7変数に時刻のスケーリングを追加することで可変時刻問題に対応できるように変更）
function eom_crtbp_w_time!(dx,x,p,t)
    # Preparation
    μ = p[1]
    r1 = x[1:3] + [μ,0,0]
    r2 = x[1:3] - [1-μ,0,0]
    a_cori = [2x[5],-2x[4],0]
    a_cf = [x[1],x[2],0]
    # dx/dt
    dx[1:3] = x[7]*(x[4:6])
    dx[4:6] = x[7]*(-(1-μ)*r1/(norm(r1)^3) - μ*r2/(norm(r2)^3) + a_cori + a_cf)
    dx[7] = 0 # 第7変数に時刻のスケーリング値を追加
end

# z0を固定した上での垂直交差(Perpendicular Crossing)の境界条件
function bc!(residual, x, p, t) # u[1] is the beginning of the time span, and u[end] is the ending
    # 初期時刻における垂直交差の境界条件
    residual[1] = x[1][2]
    residual[2] = x[1][4]
    residual[3] = x[1][6]
    # 終端時刻における垂直交差の境界条件
    residual[4] = x[end][2]
    residual[5] = x[end][4]
    residual[6] = x[end][6]
    # ハロー軌道を１つに特定するためにz0固定という条件を追加（x0固定，ヤコビ定数固定，遷移時間固定等に変えても構わない）
    residual[7] = x[1][3] - p[2] 
end

# 初期条件
x0 = [0.987384153663276, 0.0, 0.008372273063008, 0.0, 1.67419265037912, 0.0, 1.0] # 初期予想
tspan = (0.0, 0.25 * pi)
μ = 0.01215058426994 # 地球・月系のμ

# BoundaryValueDiffEqのライブラリを用いて境界値問題を解く
z0 = 0.07
bvp = BVProblem(eom_crtbp_w_time!, bc!, x0, tspan, [μ,z0])
sol = solve(bvp, Shooting(Vern7()), abstol=1e-10, reltol=1e-10)

# 描画を目的とした軌道逆伝播
tspan_neg = (0.0, -0.25 * pi)
prob = ODEProblem(eom_crtbp_w_time!, sol[:,1], tspan_neg, μ)
sol_neg = solve(prob, DP8(), reltol=1e-10, abstol=1e-10)

# 描画
fig = Figure(; resolution=(1024, 1024))
ax = Axis3(fig[1, 1]; aspect=:data, perspectiveness=0.0, title="Example of Halo Orbit", xlabel=L"x", ylabel=L"y", zlabel=L"z")
scatterlines!(ax, 1 - μ, 0, 0; markersize=20, color=:blue)
lines!(ax, sol[1, :], sol[2, :], sol[3, :], linewidth=1, color=:red)
lines!(ax, sol_neg[1, :], sol_neg[2, :], sol_neg[3, :], linewidth=1, color=:blue)
fig

境界値問題を用いたハロー軌道の計算結果

5. まとめ

本記事では，宇宙軌道力学分野で用いられる円制限三体問題を題材として，常微分方程式の数値解析ライブラリDifferentialEquations.jlの使い方を紹介をした．また，DifferentialEquations.jlの拡張機能の１つである境界値問題ライブラリBoundaryValueDiffEq.jlの使い方に関しても紹介し，ハロー軌道の生成を行った．DifferentialEquations.jlには，その他にも便利な拡張機能として，軌道伝播の途中で"ある条件"を満たした時に"ある処理"を実行するイベント処理(Event Handler) 機能があり，以下の記事で紹介しているため，参照いただきたい．

参考

Kathleen Connor Howell, "Three-dimensional, Periodic, ‘Halo’ Orbits." Celestial Mechanics, Vol.32, pp.53–71, 1984. https://doi.org/10.1007/BF01358403