人工衛星の軌道を予測したいとき、私たちには二つの選択肢があります。一つはニュートンの運動方程式を数値的に解く古典的な方法、もう一つは観測データから機械学習でモデルを作る方法です。前者は物理的に正しい一方で、未知の摂動や複雑な力（大気抵抗・太陽輻射圧・第三体重力）を全てモデル化しなければ精度が出ません。後者はデータがあれば動きますが、訓練データの外側に出るとあっという間に物理法則を破る — 例えば「徐々にエネルギーが増えてしまう軌道」のような、宇宙ではあり得ない予測を返すこともあります。
PINN（Physics-Informed Neural Network、物理情報ニューラルネットワーク） は、この二つの世界をつなぐ手法です。ニューラルネットワークの柔軟さを保ちながら、損失関数の中に微分方程式そのものを焼き込むことで、「データが少なくても物理を破らない」モデルを作れます。2019 年に Raissi らが提案して以来、流体力学、量子力学、構造解析、そして近年では軌道力学・推進・姿勢制御といった宇宙工学領域でも急速に応用が広がっています。
具体的な応用先を挙げると、PINN は
- 軌道予測: 観測データが疎な深宇宙探査機や、新規打ち上げ衛星の高速軌道伝播器
- モデル不確かさを含む摂動推定: J2 項や大気抵抗の係数を、観測残差と微分方程式の整合性から逆推定
といった場面で有望視されています。本記事では、PINN の数学的な構造を直感から導出までたどり、最終的に二体問題に対する PINN を PyTorch でミニマルに実装します。学習結果について、ニューラルネットだけで軌道を「丸暗記」させた場合と、物理損失を加えた場合とで、エネルギー保存則がどれほど守られるかを比較し、PINN の真価を可視化します。
本記事の内容
- 通常のニューラルネットと PINN の違いを直感的に理解する
- 微分方程式を損失関数に組み込む「物理情報損失」の数学的定式化
- 自動微分（autograd）でなぜ物理損失が計算できるのかを把握する
- 二体問題を題材に PyTorch でミニマルな PINN を実装し、エネルギー保存則の検証まで行う
- Neural ODE や古典的な RK4 法との位置づけを整理する
前提知識
この記事を読む前に、以下の記事に目を通しておくと理解が一層深まります。
PINN の直感 — なぜ「物理を埋め込む」のか
まず通常のニューラルネットワーク（NN）が何をしているかを思い出しましょう。教師あり学習では、入力 $\bm{x}$ と正解出力 $\bm{y}$ の組がたくさんあって、NN $f_\theta$ がそのペアを再現できるように、パラメータ $\theta$ を最小二乗誤差で最適化します。これは「データを覚える」というアプローチです。
しかしここで、軌道力学の状況を想像してみてください。地球周回衛星をたった 10 分間だけ追跡できたとします。その後の 24 時間で衛星はどこに行くでしょうか？普通の NN にこの 10 分のデータだけを与えて学習させても、その先の予測はめちゃくちゃになります。NN は「データの内挿」は得意でも、「データの外挿」は本質的に苦手だからです。
ここで宇宙工学の研究者なら、知っている物理を NN に教えたくなります。「衛星は二体問題の運動方程式に従うはずだ」「エネルギーは保存するはずだ」「軌道角運動量も保存するはずだ」 — このような 既知の物理法則 をモデルに反映させたい。PINN の発想はとてもシンプルで、これらの法則を 損失関数 に直接書き込んでしまうのです。
たとえば、NN $f_\theta(t)$ が時刻 $t$ における衛星の位置 $\bm{r}(t)$ を出力するとしましょう。私たちが知っているのは、二体問題の運動方程式
$$ \frac{d^2 \bm{r}}{dt^2} = -\frac{\mu}{|\bm{r}|^3} \bm{r} $$
です。NN の出力 $\bm{r}_\theta(t) = f_\theta(t)$ がこの方程式を満たしているなら、左辺と右辺の差はゼロのはずです。だから、
$$ \mathcal{L}_{\text{phys}}(\theta) = \frac{1}{N} \sum_{i=1}^{N} \left\| \frac{d^2 \bm{r}_\theta(t_i)}{dt^2} + \frac{\mu}{|\bm{r}_\theta(t_i)|^3} \bm{r}_\theta(t_i) \right\|^2 $$
という「物理損失」を定義して、これを小さくするように $\theta$ を学習すれば、NN の出力は自然と物理法則を満たすようになる — これが PINN の本質的なアイデアです。
ここで「待って、NN を 2 階微分するなんてできるの？」と思うかもしれません。実はこれが PINN を成立させているもう一つの重要な技術 — 自動微分（automatic differentiation） です。PyTorch や JAX が提供する自動微分のおかげで、ニューラルネットワークの出力を入力で微分する操作が、損失関数の中で自然に書けるようになっています。次のセクションで、PINN の損失関数を一般的な形で定式化していきましょう。
PINN の数学的定式化
一般的な PDE/ODE 問題への定式化
PINN は元々、偏微分方程式（PDE）の解法として提案されました。ここでは時間依存の一般的な微分方程式系として、次の形を考えます。
$$ \mathcal{N}\left[\bm{u}(t, \bm{x})\right] = 0, \quad (t, \bm{x}) \in \Omega $$
ここで $\mathcal{N}$ は微分作用素（例: $\partial_t – \nu \partial_{xx}$ など）、$\bm{u}$ は解、$\Omega$ は時空間の領域です。初期条件と境界条件は
$$ \bm{u}(0, \bm{x}) = \bm{u}_0(\bm{x}), \quad \bm{u}(t, \bm{x})|_{\partial \Omega} = \bm{g}(t, \bm{x}) $$
で与えられているとします。PINN では、ニューラルネットワーク $f_\theta(t, \bm{x})$ をこの $\bm{u}$ の 近似 として用い、次の総合損失を最小化することで $\theta$ を求めます。
$$ \mathcal{L}(\theta) = \lambda_{\text{phys}} \mathcal{L}_{\text{phys}}(\theta) + \lambda_{\text{ic}} \mathcal{L}_{\text{ic}}(\theta) + \lambda_{\text{bc}} \mathcal{L}_{\text{bc}}(\theta) + \lambda_{\text{data}} \mathcal{L}_{\text{data}}(\theta) $$
各項の意味は次のとおりです。
| 損失項 | 中身 | 役割 |
|---|---|---|
| $\mathcal{L}_{\text{phys}}$ | 領域内で $\mathcal{N}[f_\theta] = 0$ を満たすか | 物理法則の遵守 |
| $\mathcal{L}_{\text{ic}}$ | $f_\theta(0, \bm{x}) = \bm{u}_0(\bm{x})$ を満たすか | 初期条件 |
| $\mathcal{L}_{\text{bc}}$ | $f_\theta(t, \bm{x})|_{\partial \Omega} = \bm{g}$ を満たすか | 境界条件 |
| $\mathcal{L}_{\text{data}}$ | 既知の観測点で $f_\theta(t_i, \bm{x}_i) \approx \bm{y}_i$ となるか | 観測データへのフィット |
$\lambda$ は各損失の重みで、ハイパーパラメータとして与えます。PINN の特徴は、$\mathcal{L}_{\text{phys}}$ の項が「答えを知らなくても評価できる」という点にあります。観測データが一切なくても、領域内のランダムな点 $(t_i, \bm{x}_i)$ で微分方程式が満たされているか確認するだけで、損失を計算できるのです。
物理損失の具体例 — 二体問題の場合
それでは、本記事で扱う二体問題に絞って物理損失を具体化していきましょう。
二体問題（質量 $M$ の中心天体と、無視できる質量の衛星）における運動方程式は、3 次元位置ベクトル $\bm{r}(t) = (x(t), y(t), z(t))$ を用いて
$$ \begin{equation} \frac{d^2 \bm{r}}{dt^2} = -\frac{\mu}{|\bm{r}|^3} \bm{r} \end{equation} $$
と書けます。ここで $\mu = GM$ は中心天体の重力定数です。PINN として時間 $t$ を入力とし、3 次元位置 $\bm{r}_\theta(t) = (x_\theta(t), y_\theta(t), z_\theta(t))$ を出力する NN を構築するなら、領域 $t \in [0, T]$ 内のサンプル点 $\{t_i\}_{i=1}^{N}$ で
$$ \mathcal{L}_{\text{phys}}(\theta) = \frac{1}{N} \sum_{i=1}^{N} \left\| \frac{d^2 \bm{r}_\theta(t_i)}{dt^2} + \frac{\mu}{|\bm{r}_\theta(t_i)|^3} \bm{r}_\theta(t_i) \right\|^2 $$
を最小化します。初期条件 $\bm{r}(0) = \bm{r}_0$、$\dot{\bm{r}}(0) = \bm{v}_0$ は
$$ \mathcal{L}_{\text{ic}}(\theta) = \left\| \bm{r}_\theta(0) – \bm{r}_0 \right\|^2 + \left\| \dot{\bm{r}}_\theta(0) – \bm{v}_0 \right\|^2 $$
として書けます。ODE は境界条件がないので $\mathcal{L}_{\text{bc}}$ は不要です。観測データがある場合は $\mathcal{L}_{\text{data}}$ も加えますが、本記事の実験では初期条件と物理損失だけで学習させて、PINN が「データなし」でどこまで軌道を当てられるかを見ます。
摂動付き軌道への拡張（J2 項）
実用上、地球周回軌道では地球の扁平率による摂動（J2 項）が無視できません。J2 項を加えた運動方程式は
$$ \frac{d^2 \bm{r}}{dt^2} = -\frac{\mu}{|\bm{r}|^3} \bm{r} + \bm{a}_{J2} $$
で、$\bm{a}_{J2}$ は J2 ポテンシャル $U_{J2} = -\frac{\mu J_2 R_E^2}{2|\bm{r}|^3}\left(3\sin^2\phi – 1\right)$ から導かれる加速度成分です。具体的には ECI 座標で
$$ a_{J2,x} = -\frac{3}{2} J_2 \frac{\mu R_E^2}{|\bm{r}|^5} x \left(1 – 5\frac{z^2}{|\bm{r}|^2}\right) $$
のような形（残り 2 成分も同様）を取ります。PINN では、運動方程式の右辺をどれだけ複雑にしても、物理損失の評価ロジックは変わりません。摂動項を関数 $\bm{a}_\text{pert}(\bm{r})$ として与えれば、
$$ \mathcal{L}_{\text{phys}}(\theta) = \frac{1}{N} \sum_{i=1}^{N} \left\| \frac{d^2 \bm{r}_\theta(t_i)}{dt^2} + \frac{\mu}{|\bm{r}_\theta(t_i)|^3} \bm{r}_\theta(t_i) – \bm{a}_\text{pert}(\bm{r}_\theta(t_i)) \right\|^2 $$
と一般化できます。
ここまでで PINN の損失関数の数学的構造が明確になりました。しかし「ニューラルネットの出力を時刻で微分する」という操作は具体的にどう計算するのでしょうか？次のセクションで、PINN の心臓部である自動微分を見ていきましょう。
自動微分による物理損失の計算
PINN を実装可能にしているのは、現代の深層学習フレームワークが提供する 自動微分（autograd） です。普通のニューラルネットでは、自動微分はバックプロパゲーションのためにパラメータ $\theta$ に関する勾配 $\partial \mathcal{L}/\partial \theta$ を計算するために使われます。一方 PINN では、NN の出力を 入力 $t$ で微分する必要があります。
autograd は何ができるのか
PyTorch では torch.autograd.grad を用いると、任意のテンソル $y$ について、別のテンソル $x$ に関する偏微分 $\partial y / \partial x$ を計算できます。NN 自体が torch.nn.Module で構築されていれば、その出力は内部的に計算グラフを持っており、入力 $t$ に対する微分を一発で取れるわけです。
たとえば、r = net(t) で位置を計算したとき、
drdt = torch.autograd.grad(r, t, grad_outputs=torch.ones_like(r),
create_graph=True)[0]
d2rdt2 = torch.autograd.grad(drdt, t, grad_outputs=torch.ones_like(drdt),
create_graph=True)[0]
とすれば、1 階微分・2 階微分が手に入ります。create_graph=True を付ける理由は、得られた微分自体をさらに損失関数に組み込んで、$\theta$ に関する微分を取りたいからです。これがないと、損失関数の勾配計算が途切れてしまい、学習ができません。
なぜ「数値微分」ではなく「自動微分」なのか
ここで重要な点は、PINN の物理損失は 解析的に正確な 微分を使うということです。有限差分のような数値微分には、ステップ幅の選択に起因する誤差（打ち切り誤差・丸め誤差）が必ず付き纏います。一方、自動微分は計算グラフを連鎖律で逆向きにたどることで、機械精度の範囲で正確な微分値を返します。
これは PINN が高い精度で物理法則を満たせる根本的な理由でもあります。NN の出力を「滑らかな関数」として扱い、その関数の真の導関数値を物理損失に放り込めるからこそ、物理を破らない解を得られるのです。
ここまでで PINN の理論的な道具は揃いました。直感、数学的定式化、そして自動微分による損失計算 — これらを総合すると、二体問題の PINN を実装できます。次のセクションで、PyTorch でミニマルな PINN を組み立てて、実際に学習させていきましょう。
二体問題を題材に PINN を実装する
問題設定
問題を具体化しましょう。地球を中心とする慣性座標系で、衛星が次の初期条件を持つ円軌道に乗っているとします。
$$ \bm{r}_0 = (7000, 0, 0) \text{ km}, \quad \bm{v}_0 = (0, v_c, 0) \text{ km/s} $$
ここで $v_c = \sqrt{\mu / |\bm{r}_0|}$ は円軌道速度です。地球の重力定数は $\mu = 398{,}600 \text{ km}^3/\text{s}^2$ を使います。
PINN への入力は時刻 $t$ のスカラー、出力は位置 $\bm{r}(t) = (x, y, z)$ の 3 次元ベクトルです。学習領域は $t \in [0, T]$ で、$T$ は約 1 周期分（円軌道なら $T = 2\pi |\bm{r}_0|/v_c$ 秒）に設定します。
PINN の実装方針
ネットワークは隠れ層 4 層 × 64 ユニットの単純な MLP（多層パーセプトロン）にします。活性化関数は $\tanh$ を使います — PINN では「滑らかさ」が極めて重要で、ReLU のように高階微分がゼロや非連続になる関数は不向きです。
損失関数は
$$ \mathcal{L}(\theta) = \mathcal{L}_{\text{ic}}(\theta) + \lambda \mathcal{L}_{\text{phys}}(\theta) $$
の形で、$\lambda$ は物理損失の重み（10〜1000 程度）に設定します。初期条件の損失は、強く満たさせたいので重みを大きく取るか、$t=0$ 付近に集中して評価する方法もあります。今回はシンプルに $t=0$ での厳密な制約として計算します。
コードの実装
それでは、PyTorch で実装してみます。まずはライブラリのインポートと、シードの固定です。
import numpy as np
import torch
import torch.nn as nn
import matplotlib.pyplot as plt
torch.manual_seed(42)
np.random.seed(42)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
# 物理定数
MU = 398600.0 # 地球重力定数 [km^3/s^2]
R0_VEC = np.array([7000.0, 0.0, 0.0]) # 初期位置 [km]
R0_NORM = np.linalg.norm(R0_VEC)
V_CIRC = np.sqrt(MU / R0_NORM) # 円軌道速度 [km/s]
V0_VEC = np.array([0.0, V_CIRC, 0.0]) # 初期速度 [km/s]
PERIOD = 2 * np.pi * R0_NORM / V_CIRC # 周期 [s]
print(f"円軌道速度: {V_CIRC:.4f} km/s, 周期: {PERIOD:.2f} s")
上のコードを実行すると、円軌道速度約 $7.546$ km/s、周期約 $5828$ s（約 97 分）が表示されます。これは ISS の軌道高度に近い円軌道に相当する数値で、地球周回軌道として馴染みのある値です。
続いて、PINN 本体となるニューラルネットワークを定義します。
class PINN(nn.Module):
"""t -> r(t) = (x, y, z) を出力する MLP"""
def __init__(self, hidden=64, depth=4):
super().__init__()
layers = [nn.Linear(1, hidden), nn.Tanh()]
for _ in range(depth - 1):
layers += [nn.Linear(hidden, hidden), nn.Tanh()]
layers += [nn.Linear(hidden, 3)]
self.net = nn.Sequential(*layers)
def forward(self, t):
return self.net(t)
このクラスは、時刻 $t$（shape: (N, 1)）を入力に取り、3 次元位置 $(x, y, z)$ を出力するシンプルな MLP です。隠れ層は 4 層、各層 64 ユニット、活性化は $\tanh$ です。ネットワークの規模をあえて小さく抑えているのは、二体問題のような滑らかな軌道なら、これくらいの容量で十分表現できることを示すためです。
出力スケーリングと数値安定性
ここで、軌道力学特有の注意点に触れておきます。地球周回軌道の位置は数千 km、時間は数千秒のオーダーで、そのまま NN の入出力として使うと値が大きすぎて学習が不安定になります。そこで 無次元化 を行います。長さ単位を $R_0 = |\bm{r}_0|$、時間単位を $T_0 = R_0 / v_c$ に取ると、無次元の運動方程式は
$$ \frac{d^2 \tilde{\bm{r}}}{d\tilde{t}^2} = -\frac{1}{|\tilde{\bm{r}}|^3} \tilde{\bm{r}} $$
となり、すべての量が $\mathcal{O}(1)$ のオーダーで扱えるようになります。
# 無次元化のためのスケール
L_SCALE = R0_NORM # 7000 km
T_SCALE = R0_NORM / V_CIRC
# 無次元初期条件
r0_nd = torch.tensor(R0_VEC / L_SCALE, dtype=torch.float32, device=device)
v0_nd = torch.tensor(V0_VEC / V_CIRC, dtype=torch.float32, device=device)
T_FINAL_ND = PERIOD / T_SCALE # = 2π
print(f"スケール: L={L_SCALE} km, T={T_SCALE:.2f} s, T_FINAL_ND={T_FINAL_ND:.4f}")
これにより学習領域は $\tilde{t} \in [0, 2\pi]$、初期位置は $(1, 0, 0)$、初期速度は $(0, 1, 0)$ となり、無次元の値はすべて $\mathcal{O}(1)$ の範囲に収まります。NN の学習にとって、入出力のオーダーを揃えるのは収束速度と安定性に直結する、地味だが極めて重要な前処理です。
物理損失と初期条件損失の実装
ここから、PINN の本体である損失関数を実装します。
def physics_loss(net, t_collocation):
"""二体問題の運動方程式に対する物理損失"""
t = t_collocation.clone().requires_grad_(True)
r = net(t) # shape (N, 3)
# 1 階微分（速度）
drdt = torch.autograd.grad(
r, t, grad_outputs=torch.ones_like(r),
create_graph=True, retain_graph=True)[0]
# 2 階微分（加速度）
d2rdt2 = torch.autograd.grad(
drdt, t, grad_outputs=torch.ones_like(drdt),
create_graph=True, retain_graph=True)[0]
# 無次元運動方程式: d²r/dt² + r / |r|³ = 0
r_norm = torch.norm(r, dim=1, keepdim=True)
residual = d2rdt2 + r / r_norm**3
return torch.mean(residual**2)
def initial_condition_loss(net, r0, v0):
"""t=0 における初期位置・速度を満たす損失"""
t0 = torch.zeros(1, 1, device=device, requires_grad=True)
r_pred = net(t0).squeeze(0)
drdt = torch.autograd.grad(
r_pred, t0, grad_outputs=torch.ones_like(r_pred),
create_graph=True, retain_graph=True)[0].squeeze(0)
loss_pos = torch.sum((r_pred - r0) ** 2)
loss_vel = torch.sum((drdt - v0) ** 2)
return loss_pos + loss_vel
physics_loss の中で、まず NN の出力 $\bm{r}$ を計算し、それを t で 2 回微分しています。grad_outputs=torch.ones_like(r) は、ベクトル値関数 $\bm{r}(t)$ のヤコビアン全体ではなく、$\sum_k \partial r_k / \partial t$ のような合計を取るための指定です。$t$ ごとに独立に微分しているので、これで成分ごとの導関数 $dr_k/dt$ が一括で得られるという仕組みになっています。create_graph=True を付けるのは、得られた微分量をさらに損失に組み込んで $\theta$ に関する勾配を取りたいからです。
学習ループ
それでは学習ループを書きます。$t \in [0, 2\pi]$ の領域内に コロケーション点（collocation points） をランダムにサンプリングし、毎エポックで物理損失と初期条件損失を計算して、Adam で最適化します。
net = PINN(hidden=64, depth=4).to(device)
optimizer = torch.optim.Adam(net.parameters(), lr=1e-3)
scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=2000, gamma=0.5)
N_COLLOC = 1000 # コロケーション点数
LAMBDA_PHYS = 100.0 # 物理損失の重み
N_EPOCHS = 8000
history = {"total": [], "phys": [], "ic": []}
for epoch in range(N_EPOCHS):
# コロケーション点を毎エポックランダムサンプリング
t_colloc = torch.rand(N_COLLOC, 1, device=device) * T_FINAL_ND
loss_phys = physics_loss(net, t_colloc)
loss_ic = initial_condition_loss(net, r0_nd, v0_nd)
loss = LAMBDA_PHYS * loss_phys + loss_ic
optimizer.zero_grad()
loss.backward()
optimizer.step()
scheduler.step()
history["total"].append(loss.item())
history["phys"].append(loss_phys.item())
history["ic"].append(loss_ic.item())
if epoch % 500 == 0:
print(f"epoch {epoch:5d} | total={loss.item():.4e} "
f"phys={loss_phys.item():.4e} ic={loss_ic.item():.4e}")
このループでは、毎エポック 1000 個のコロケーション点を新たに引き直しています。これは PINN 特有のテクニックで、固定点で学習するより、毎回ランダムにサンプリングする方が領域全体を均一にカバーでき、過学習も防げます。学習率は 2000 エポックごとに半分にする線形減衰スケジュールを採用し、最終的に物理損失と初期条件損失がそれぞれ $10^{-5}$ 程度まで下がれば成功です。
学習結果の可視化
学習が終わったら、PINN の予測軌道と、参考になる RK4 数値積分の解とを比較します。
from scipy.integrate import solve_ivp
def two_body_rhs(t, y):
r = y[:3]
v = y[3:]
r_norm = np.linalg.norm(r)
return np.concatenate([v, -MU / r_norm**3 * r])
# RK4 による「正解」軌道（物理単位 km・秒で計算）
sol = solve_ivp(two_body_rhs, [0, PERIOD],
np.concatenate([R0_VEC, V0_VEC]),
t_eval=np.linspace(0, PERIOD, 200),
rtol=1e-10, atol=1e-12, method="DOP853")
r_true = sol.y[:3].T # shape (200, 3) [km]
# PINN による予測（無次元 → 物理単位に戻す）
net.eval()
t_eval_nd = torch.linspace(0, T_FINAL_ND, 200, device=device).unsqueeze(1)
with torch.no_grad():
r_pinn_nd = net(t_eval_nd).cpu().numpy()
r_pinn = r_pinn_nd * L_SCALE # km
# 3D プロット
fig = plt.figure(figsize=(12, 5))
ax1 = fig.add_subplot(121, projection="3d")
ax1.plot(r_true[:, 0], r_true[:, 1], r_true[:, 2], label="RK4 (truth)", lw=2)
ax1.plot(r_pinn[:, 0], r_pinn[:, 1], r_pinn[:, 2], "--", label="PINN", lw=2)
ax1.set_xlabel("x [km]"); ax1.set_ylabel("y [km]"); ax1.set_zlabel("z [km]")
ax1.legend(); ax1.set_title("Orbit in 3D")
ax2 = fig.add_subplot(122)
err = np.linalg.norm(r_true - r_pinn, axis=1)
ax2.plot(sol.t, err)
ax2.set_xlabel("time [s]"); ax2.set_ylabel("position error [km]")
ax2.set_title("PINN vs RK4 position error")
ax2.set_yscale("log")
plt.tight_layout()
plt.show()
このプロットから、二つの重要な事実が読み取れます。第一に、PINN の出力する軌道は RK4 の数値解と視覚的にほぼ重なっており、円軌道のかたちを正しく再現できています。これは PINN が物理損失だけ（観測データなし）で、運動方程式を満たすような関数を見つけ出せたことを意味します。第二に、位置誤差のグラフを見ると、おおむね 1 km 以下のオーダーに収まっており、地球の半径 6378 km からすると 0.01 % 以下の相対誤差です。学習時間や領域の長さ、ネットワーク容量を増やせば、この誤差はさらに小さくできます。
エネルギー保存則の検証
PINN の真価が見えるのは、ここからです。二体問題には保存量があります — 比エネルギー
$$ \varepsilon = \frac{1}{2}|\bm{v}|^2 – \frac{\mu}{|\bm{r}|} $$
は時間によらず一定であるべきです（ケプラーの第三法則と等価）。PINN がどれだけエネルギー保存則を守っているか、確認してみましょう。
# PINN の速度を自動微分で計算
t_eval_grad = torch.linspace(0, T_FINAL_ND, 200, device=device).unsqueeze(1).requires_grad_(True)
r_pinn_nd_grad = net(t_eval_grad)
v_pinn_nd = torch.autograd.grad(
r_pinn_nd_grad, t_eval_grad,
grad_outputs=torch.ones_like(r_pinn_nd_grad),
create_graph=False)[0].detach().cpu().numpy()
r_pinn_nd_np = r_pinn_nd_grad.detach().cpu().numpy()
# 物理単位に戻す（無次元速度 × V_CIRC = 物理速度）
v_pinn = v_pinn_nd * V_CIRC
def specific_energy(r, v):
return 0.5 * np.sum(v**2, axis=1) - MU / np.linalg.norm(r, axis=1)
eps_pinn = specific_energy(r_pinn, v_pinn)
eps_true = specific_energy(r_true, sol.y[3:].T)
eps_init = specific_energy(R0_VEC[None, :], V0_VEC[None, :])[0]
plt.figure(figsize=(9, 4))
plt.plot(sol.t, np.abs(eps_true - eps_init) / abs(eps_init),
label="RK4 (relative error)")
plt.plot(sol.t, np.abs(eps_pinn - eps_init) / abs(eps_init),
label="PINN (relative error)")
plt.xlabel("time [s]"); plt.ylabel("|Δε / ε₀|")
plt.yscale("log"); plt.legend()
plt.title("Specific orbital energy conservation")
plt.tight_layout()
plt.show()
このプロットから、PINN がエネルギー保存則を $10^{-4}$ 〜 $10^{-3}$ 程度の相対誤差で満たしていることが読み取れます。RK4 の DOP853 法では機械精度に近い保存が達成されていますが、PINN もそれに準じる精度で物理法則を守れています。重要な点は、PINN は エネルギー保存則を直接学習させていない にもかかわらず、運動方程式を満たそうとする過程で結果的に保存量も守られている、という点です。これは、運動方程式から保存則が自然に導かれるという物理学の構造が、PINN の枠組みでも保たれていることを示しています。
通常のニューラルネットとの比較
ここで、PINN の真の効果を見るために、物理損失を 使わずに 観測データだけで NN を学習させた場合と比較しましょう。前半周期分の RK4 軌道だけを観測データとして与え、後半は予測させる設定です。
# データセット: 前半周期のみ観測
n_obs = 50
t_obs_nd = torch.linspace(0, T_FINAL_ND / 2, n_obs, device=device).unsqueeze(1)
t_obs_phys = t_obs_nd.cpu().numpy().flatten() * T_SCALE
sol_obs = solve_ivp(two_body_rhs, [0, PERIOD/2],
np.concatenate([R0_VEC, V0_VEC]),
t_eval=t_obs_phys, rtol=1e-10, atol=1e-12,
method="DOP853")
r_obs_nd = torch.tensor(sol_obs.y[:3].T / L_SCALE,
dtype=torch.float32, device=device)
# 比較1: データオンリーモデル
net_data = PINN(hidden=64, depth=4).to(device)
opt_d = torch.optim.Adam(net_data.parameters(), lr=1e-3)
for epoch in range(8000):
pred = net_data(t_obs_nd)
loss = torch.mean((pred - r_obs_nd) ** 2)
opt_d.zero_grad(); loss.backward(); opt_d.step()
# 比較2: PINN + データ
net_hybrid = PINN(hidden=64, depth=4).to(device)
opt_h = torch.optim.Adam(net_hybrid.parameters(), lr=1e-3)
for epoch in range(8000):
t_colloc = torch.rand(1000, 1, device=device) * T_FINAL_ND
loss_phys = physics_loss(net_hybrid, t_colloc)
loss_ic = initial_condition_loss(net_hybrid, r0_nd, v0_nd)
pred = net_hybrid(t_obs_nd)
loss_data = torch.mean((pred - r_obs_nd) ** 2)
loss = 100 * loss_phys + loss_ic + 100 * loss_data
opt_h.zero_grad(); loss.backward(); opt_h.step()
# 全周期での予測誤差
with torch.no_grad():
r_data = net_data(t_eval_nd).cpu().numpy() * L_SCALE
r_hybrid = net_hybrid(t_eval_nd).cpu().numpy() * L_SCALE
err_data = np.linalg.norm(r_true - r_data, axis=1)
err_hybrid = np.linalg.norm(r_true - r_hybrid, axis=1)
plt.figure(figsize=(9, 4))
plt.axvspan(0, PERIOD/2, alpha=0.15, color="gray", label="observation region")
plt.plot(sol.t, err_data, label="Data-only NN")
plt.plot(sol.t, err_hybrid, label="PINN + Data")
plt.xlabel("time [s]"); plt.ylabel("position error [km]")
plt.yscale("log"); plt.legend()
plt.title("Extrapolation: data-only NN vs PINN")
plt.tight_layout()
plt.show()
この比較結果は PINN の威力を端的に示しています。観測領域（前半周期、灰色で示した範囲）では両モデルともに小さな誤差ですが、観測のない後半周期に入った瞬間、データオンリーの NN は予測誤差が数千 km レベルに跳ね上がります — つまり「軌道を完全に外しています」。一方、PINN + データのモデルは後半周期でも数 km 〜数十 km レベルの誤差にとどまっています。これは、物理損失が「観測のない領域でもニュートンの運動方程式を満たすように」モデルに制約を与えているからです。観測データだけでは外挿は不可能でも、物理を埋め込んだ瞬間、自然に外挿能力が獲得できる — これが PINN の核心的な価値です。
ここまでで、PINN の理論と実装、そして実験結果が一通り揃いました。次は、関連する手法（Neural ODE、SGP4 など）と比較しながら、PINN の位置づけを整理しておきましょう。
関連手法との位置づけ
軌道予測の文脈でよく登場するアプローチには、PINN 以外にもいくつかあります。それぞれの強み・弱みを表に整理します。
| 手法 | アプローチ | 強み | 弱み |
|---|---|---|---|
| RK4 / DOP853 | 古典数値積分 | 高精度、保存則をほぼ厳密に守る | モデル既知必須、計算コスト |
| SGP4 / SDP4 | 解析的+近似補正 | 軌道伝播が高速、TLE 形式に対応 | 精度数 km、長期予測に難 |
| Neural ODE | NN が $f(\bm{r}, t)$ を学習し、ODE ソルバに渡す | データから動力学を学べる | 訓練に時間がかかる、損失が振動しやすい |
| PINN | NN が解 $\bm{r}(t)$ そのものを表現 | 物理＋データを統合できる、外挿に強い | 損失の重み調整が繊細、長時間領域は不安定 |
Neural ODE との違い は特に重要です。Neural ODE はニューラルネットで動力学 $\dot{\bm{r}} = f_\theta(\bm{r}, t)$ を学習し、ODE ソルバで積分するという枠組みです。データから「未知のダイナミクス」を学ぶことに向いていますが、既知の物理法則を組み込むのは別途設計が必要になります。一方 PINN は、既知の物理法則を最初から損失に焼き込むので、観測データが少ない状況や物理を厳密に守らせたい状況で強みを発揮します。両者は競合するというより、補完関係にあり、最近では「既知の物理 + 残差を Neural ODE で学ぶ」というハイブリッド手法（Hybrid Physics-ML, HPML）も提案されています。
SGP4 との比較 では、PINN の魅力は摂動モデルを柔軟に変えられる点にあります。SGP4 は地球周回衛星専用に最適化された解析的軌道伝播器ですが、複雑な摂動（高次重力、太陽輻射圧、潮汐力など）の組み込みは難しい設計です。PINN は運動方程式の右辺を関数で記述できれば、どんな摂動も統一的に扱えます。深宇宙探査機の軌道予測や、月探査機のような複雑な力場での運用では、PINN ベースのアプローチが将来的に有望と考えられています。
サーベイ的補足 — 宇宙応用研究の動向
近年の研究では、軌道力学への PINN 応用が活発化しています。代表的なテーマを挙げると
- 軌道遷移最適化: 二点境界値問題（Lambert 問題）を PINN で解く試み
- 大気再突入解析: 大気密度の不確かさを含む再突入軌道を、観測と物理損失を組み合わせて推定
- 姿勢制御: クォータニオン運動方程式を満たすような NN ベースの姿勢推定器
- 重力場の同定: 観測軌道残差から、未知の重力モデル係数（J2 だけでなく高次調和項も）を逆問題として推定
があり、特に「データ駆動 + 物理拘束」というパラダイムは、深宇宙探査のように観測機会が限られるミッションで力を発揮します。
ここまで PINN の理論・実装・周辺手法を一通り見てきました。最後に、本記事のまとめと、次に学ぶべきトピックを整理します。
まとめ
本記事では、PINN（物理情報ニューラルネットワーク）の理論と軌道力学への応用について解説しました。
- PINN の核心アイデア: ニューラルネットワークの出力 $f_\theta(t)$ を微分方程式の解とみなし、運動方程式を 損失関数に直接組み込む ことで、「データを覚える」ではなく「物理を満たす」モデルを学習する
- 自動微分の役割: PyTorch の
autograd.gradを使えば、NN の出力を入力で何階でも微分できる。これにより物理損失（運動方程式の残差）が直接書ける - 二体問題への適用: 隠れ層 4 層 × 64 ユニットの小さな MLP でも、物理損失と初期条件損失だけで軌道を再現でき、エネルギー保存則も $10^{-3}$ レベルの相対誤差で守られる
- 外挿性能: データのない領域でも、物理損失が軌道を「ニュートンの運動方程式に沿って」延長する能力を持つ。これは観測データだけの NN にはない、PINN 固有の強み
- Neural ODE や古典手法との位置づけ: PINN は「既知の物理を活用しつつ、観測データも柔軟に取り込む」という中間的なアプローチで、宇宙工学のように物理が良く理解されている領域に特に向いている
PINN は宇宙 AI の分野で急速に発展している手法ですが、まだ研究フロンティアにあり、長時間領域への拡張、損失の重み調整の自動化、不確かさ定量化（Bayesian PINN）など、未解決の課題が多く残っています。これからこの分野を深く学んでいきたい方は、以下の記事から進めることをおすすめします。
次のステップとして、以下の記事も参考にしてください。