衛星の熱設計を担当する人なら、こんな経験があるはずです。数十万要素のFEMモデルで非定常熱伝導を1軌道分（90分）解こうとしたら、PCが1晩動き続けてもまだ終わらない。スリュー機動中の太陽光入射が変わるたびにメッシュを切り直す必要があり、設計変更のループが回らない。さらに、軌道上で熱電対の実測値が予測とずれたとき、モデルパラメータをどう更新すればよいのか分からない。これらは熱設計の現場でよく聞く悩みです。
PINN（Physics-Informed Neural Network） は、こうした課題に新しい道筋を示します。Raissi らが2019年に提案したこの手法は、ニューラルネットワーク $T_\theta(\bm{x}, t)$ を「PDEの連続な解」そのものとして学習させるアプローチです。FEMがメッシュ上の節点値を解く離散手法であるのに対し、PINNは座標 $(\bm{x}, t)$ を入力として温度や応力を直接出力する連続関数を学習します。微分はメッシュではなく自動微分で取るため、メッシュ品質に縛られず、高次元・複雑形状にも比較的素直にスケールします。さらに、損失関数に観測データの誤差を加えれば、観測と物理を同時に満たすモデルを得られる — これは軌道上でのモデル更新（デジタルツインの実時間補正）に直結します。
本記事の応用先を具体的に挙げておきます。
- 衛星バス熱予測: メッシュレスな連続解として、複雑なラジエータ配置と熱経路を扱う。設計変更のたびにメッシュを切り直す必要がない
- 軌道上のリアルタイム熱応力監視: 熱電対の疎なデータから、観測されていない位置の温度・応力を物理拘束付きで推定する
- モデル更新による軌道上適応: 熱伝導率やヒートパイプ性能の劣化を、観測残差と PDE 残差の整合から逆推定する
本記事の内容
- 偏微分方程式に対する PINN の定式化（Raissi 2019 の枠組み）
- 熱伝導方程式と構造方程式のPINNロス設計（PDE残差・境界条件・初期条件・観測データ）
- 自動微分で物理ロスがなぜ計算できるか、JAX による高次微分の取り扱い
- アーキテクチャ選択（FCNN・SIREN・Fourier features）と学習を安定化させる工夫
- 1D非定常熱伝導PINNの実装と解析解との比較
- 2D衛星パネル熱問題のPINN実装と、FEMとの計算時間・精度比較
前提知識
この記事は PINN の基礎から熱・構造PDEの定式化、JAX/PyTorch 実装までを扱います。以下の記事に目を通しておくと、流れを追いやすくなります。
- PINN（物理情報ニューラルネットワーク）の理論と軌道力学への応用
- 熱応力 — 温度差が構造に与える応力の理論と実装
- FEM（有限要素法）入門 — 弾性体の離散化を1から理解する
- 再突入カプセルの熱構造連成解析
- 衛星熱制御 — MLI・ヒートパイプ・ラジエータの設計
PINN を熱・構造に拡張する直感
軌道力学 PINN の記事では、入力 $t$ から状態 $\bm{r}(t)$ を出すネットワークを、運動方程式 $\ddot{\bm{r}} = -\mu \bm{r}/|\bm{r}|^3$ の残差をロスにして学習しました。常微分方程式（ODE）が舞台でしたから、入力は時間 $t$ ひとつでした。
熱伝導や構造解析の舞台は偏微分方程式（PDE）です。温度や変位は時間と空間の両方の関数になります。たとえば衛星パネル上の温度 $T(x, y, t)$ は、位置 $(x, y)$ と時刻 $t$ の関数で、これが熱伝導方程式
$$ \rho c \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) + Q $$
に従って時々刻々と変化します。FEMはこれをメッシュ上の節点値 $T_i^n$ に離散化して解きますが、PINNでは
$$ T_\theta : (x, y, t) \;\mapsto\; T $$
という連続関数としてニューラルネットを学習させます。学習が完了したネットワーク $T_\theta$ は、任意の点 $(x, y, t)$ における温度を返す「物理を満たす補間関数」になります。
ここまで読んで、当然こんな疑問が湧くはずです。FEMは数十年かけて磨かれた成熟手法で、商用ソルバーも揃っているのに、なぜいまさら PINN を使うのか？答えは「FEMが苦手な場面で PINN が光るから」です。
- 高次元 / メッシュレス: 3次元時間依存問題でも、座標を増やすだけでネットワークの入力次元が増えるだけ。メッシュ生成という人手の作業が要らない
- 逆問題が同じ枠組み: 物性パラメータ（熱伝導率 $k$ など）を未知数として学習させると、観測データと PDE 残差から逆推定できる。FEMで逆問題を解こうとすると、随伴法など別の重い実装が必要
- 微分可能: 学習済み $T_\theta$ は座標について微分可能なので、設計変数に対する感度解析が自動微分で得られる。形状最適化との相性がよい
- データ同化: 観測点が疎でも、PDE 残差が「物理的に整合する解」を支えてくれる。データだけの NN は外挿で破綻するが、PINN は物理が外挿の道しるべになる
逆に、メッシュ精度が出るところで FEM と精度競争をすると PINN は分が悪い場面が多くあります。本記事の立ち位置は「PINN は FEM を置き換えるものではなく、補強・補完する道具」です。
ここまでで PINN を PDE に拡張する動機が見えました。次に、その数学的な定式化を Raissi 2019 の枠組みで整理しましょう。
PINN の PDE 定式化（Raissi 2019）
一般形
一般の PDE を次のように書きます。
$$ \begin{equation} \mathcal{N}[u](\bm{x}, t) = 0, \quad (\bm{x}, t) \in \Omega \times (0, T] \end{equation} $$
ここで $u(\bm{x}, t)$ が未知の解、$\mathcal{N}$ は微分演算子です。たとえば1次元熱伝導なら $\mathcal{N}[u] = \partial_t u – \alpha \partial_{xx} u – q$ という形になります。これに境界条件と初期条件が付随します。
$$ \mathcal{B}[u](\bm{x}, t) = 0, \quad (\bm{x}, t) \in \partial\Omega \times (0, T] $$ $$ u(\bm{x}, 0) = u_0(\bm{x}), \quad \bm{x} \in \Omega $$
PINN では、未知の $u$ をニューラルネットワーク $u_\theta(\bm{x}, t)$ で近似します。ここで $\theta$ は重みとバイアスの集合です。学習が目指すのは、$u_\theta$ が PDE と境界条件と初期条件のすべてを「同時にできるだけ満たす」ようにパラメータ $\theta$ を選ぶことです。
損失関数の構成
そのために、各条件をそれぞれ二乗ノルムで評価した損失項を作り、重み付き和を取ります。
$$ \begin{equation} \mathcal{L}(\theta) = \lambda_f \mathcal{L}_f + \lambda_b \mathcal{L}_b + \lambda_0 \mathcal{L}_0 + \lambda_d \mathcal{L}_d \end{equation} $$
各項の中身を順に書きます。
PDE 残差項。領域内部に散布したコロケーション点 $\{(\bm{x}_f^{(i)}, t_f^{(i)})\}_{i=1}^{N_f}$ の上で、PDE の残差をゼロに近づけます。
$$ \mathcal{L}_f = \frac{1}{N_f} \sum_{i=1}^{N_f} \left| \mathcal{N}[u_\theta](\bm{x}_f^{(i)}, t_f^{(i)}) \right|^2 $$
ここで重要なのは、$\mathcal{N}[u_\theta]$ を計算するには $u_\theta$ の偏微分が必要なことです。これを自動微分で取るのが PINN の魂です。たとえば1次元熱伝導の残差 $\partial_t u_\theta – \alpha \partial_{xx} u_\theta – q$ は、$u_\theta(x, t)$ を $t$ について1回、$x$ について2回 autograd で微分するだけで計算できます。
境界条件項。境界 $\partial\Omega$ 上の点 $\{(\bm{x}_b^{(i)}, t_b^{(i)})\}_{i=1}^{N_b}$ で境界条件を評価します。
$$ \mathcal{L}_b = \frac{1}{N_b} \sum_{i=1}^{N_b} \left| \mathcal{B}[u_\theta](\bm{x}_b^{(i)}, t_b^{(i)}) \right|^2 $$
Dirichlet 境界 $u = g$ なら $|u_\theta – g|^2$、Neumann 境界 $\partial_n u = q$ なら $|\partial_n u_\theta – q|^2$ になります。
初期条件項。初期時刻 $t = 0$ 上の点で初期条件を評価します。
$$ \mathcal{L}_0 = \frac{1}{N_0} \sum_{i=1}^{N_0} \left| u_\theta(\bm{x}_0^{(i)}, 0) – u_0(\bm{x}_0^{(i)}) \right|^2 $$
観測データ項。観測値 $\{u_d^{(i)}\}_{i=1}^{N_d}$ がある場合に追加します。
$$ \mathcal{L}_d = \frac{1}{N_d} \sum_{i=1}^{N_d} \left| u_\theta(\bm{x}_d^{(i)}, t_d^{(i)}) – u_d^{(i)} \right|^2 $$
宇宙機の温度モニタリングでは、熱電対の値がこの項に入ります。観測点が疎でも、PDE 残差が「物理的に妥当な補間」を支えてくれる、というのが PINN の強みでした。
なぜ自動微分で物理が入るのか
ニューラルネットの順伝播 $u_\theta(\bm{x}, t)$ は、入力に対して合成関数の連なりです。これを座標について微分するということは、合成関数の連鎖律を遡るだけで、autograd の機能としてフレームワーク（JAX, PyTorch）にそのまま備わっています。重要なのは、この微分が解析的に正確（メッシュによる近似ではない）で、しかも $\theta$ についてもさらに微分可能だということです。
つまり、$\mathcal{L}_f$ は $u_\theta$ の高次微分を含む式ですが、それでも $\theta$ について微分可能で、勾配降下が回せます。これが、メッシュレスで PDE を解ける根拠です。
ここまでで PINN の損失関数の構造が整理できました。次に、本記事の主役である熱伝導方程式と構造方程式に、この枠組みを当てはめていきましょう。
熱伝導方程式の PINN 定式化
衛星熱方程式
衛星バスやペイロード筐体の温度場 $T(\bm{x}, t)$ は、次の熱伝導方程式に従います。
$$ \begin{equation} \rho c \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) + Q \end{equation} $$
ここで $\rho$ は密度 $[\text{kg/m}^3]$、$c$ は比熱 $[\text{J/(kg·K)}]$、$k$ は熱伝導率 $[\text{W/(m·K)}]$、$Q$ は内部発熱密度 $[\text{W/m}^3]$ です。$k$ が一定なら、両辺を $\rho c$ で割って
$$ \frac{\partial T}{\partial t} = \alpha \nabla^2 T + \frac{Q}{\rho c}, \quad \alpha = \frac{k}{\rho c} $$
の形になります。$\alpha$ は熱拡散率です。アルミニウム合金（衛星構体の主材料）なら $\alpha \sim 7 \times 10^{-5}$ m²/s、CFRP なら $\alpha \sim 5 \times 10^{-7}$ m²/s 程度です。
境界条件の種類
衛星の熱境界は多彩で、典型的には次の3種が組み合わさります。
- Dirichlet 境界（温度固定）: ヒーターで一定温度に保たれている面
- Neumann 境界（熱流束）: 太陽光入射 $q_{\text{sun}}$、地球アルベド、地球赤外
- Robin 境界（対流・放射）: ラジエータからの宇宙放射 $q = \varepsilon \sigma (T^4 – T_\infty^4)$
放射境界は $T^4$ の非線形性が入り、FEM では反復計算が必要になりますが、PINN では損失関数に書き込むだけで自然に扱えます。これは PINN の隠れた強みのひとつです。
PINN ロスの具体化
ネットワーク $T_\theta(\bm{x}, t)$ を用意します。コロケーション点を領域内部に散布し、PDE 残差
$$ r_f(\bm{x}, t) = \rho c \frac{\partial T_\theta}{\partial t} – \nabla \cdot (k \nabla T_\theta) – Q $$
の二乗平均が PDE ロス $\mathcal{L}_f$ になります。境界では境界条件タイプに応じて
$$ r_b = \begin{cases} T_\theta – T_{\text{bc}} & \text{Dirichlet} \\ k \partial_n T_\theta – q_{\text{bc}} & \text{Neumann} \\ k \partial_n T_\theta + \varepsilon \sigma (T_\theta^4 – T_\infty^4) & \text{放射} \end{cases} $$
を選びます。Robin/放射境界では、境界点で $T^4$ を含む非線形項も自動微分で問題なく計算できます。
無次元化のコツ
実装上の注意点として、温度や座標の単位を揃えると学習が安定します。物理量はそのままだと $T \sim 300$ K、$x \sim 0.1$ m、$t \sim 100$ s と桁が違うため、損失項間の重みのバランスを取りにくくなります。次のような無次元化が定番です。
$$ \tilde{x} = \frac{x}{L_0}, \quad \tilde{t} = \frac{\alpha t}{L_0^2}, \quad \tilde{T} = \frac{T – T_0}{\Delta T} $$
これで方程式は $\partial_{\tilde t} \tilde T = \nabla_{\tilde x}^2 \tilde T + \tilde Q$ という綺麗な形になります。Fourier 数 $\mathrm{Fo} = \alpha t / L_0^2$ を新しい時間軸として使う、と覚えても同じことです。
熱の話はこれで一通り見通せました。次は構造方程式 — 温度差が引き起こす熱応力を扱う準備をしましょう。
構造方程式の PINN 定式化
弾性体の支配方程式
構造解析では、変位場 $\bm{u}(\bm{x})$ を未知数として、釣り合い式と構成則を解きます。準静的な微小変形ならば、
$$ \begin{equation} \nabla \cdot \bm{\sigma} + \bm{f} = \bm{0} \end{equation} $$
が釣り合い式（運動方程式の慣性項を落としたもの）です。$\bm{\sigma}$ は応力テンソル、$\bm{f}$ は体積力（重力など）。弾性体の構成則は線形弾性なら
$$ \begin{equation} \bm{\sigma} = \bm{C} : \bm{\varepsilon} \end{equation} $$
で、$\bm{C}$ は弾性係数テンソル（4階）、$\bm{\varepsilon} = (\nabla \bm{u} + \nabla \bm{u}^T)/2$ がひずみテンソルです。「$:$」は2階テンソルの2重内積（成分ごとに掛けて足す）を意味します。
熱応力を含む形
衛星構造では温度変化が応力を生みます。温度差 $\Delta T = T – T_{\text{ref}}$ による熱ひずみを加えると、
$$ \bm{\varepsilon} = \bm{\varepsilon}_{\text{el}} + \bm{\varepsilon}_{\text{th}}, \quad \bm{\varepsilon}_{\text{th}} = \alpha_T \Delta T \, \bm{I} $$
応力は弾性ひずみのみに比例するので、
$$ \bm{\sigma} = \bm{C} : (\bm{\varepsilon} – \bm{\varepsilon}_{\text{th}}) $$
これが熱応力の基本式です。等方均質材なら $C_{ijkl} = \lambda \delta_{ij}\delta_{kl} + \mu(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$ で、$\lambda, \mu$ がラメ定数になります。
構造 PINN の損失
ネットワーク $\bm{u}_\theta(\bm{x})$ を用意します。位置 $\bm{x}$ で $\bm{u}_\theta$ を1回微分すると $\bm{\varepsilon}$ が、もう1回微分すると応力勾配 $\nabla \cdot \bm{\sigma}$ が出ます。釣り合い式の残差を二乗して PDE ロスとし、境界では応力境界（Neumann）または変位境界（Dirichlet）を別途ロスに足します。
弱形式（エネルギー汎関数）を直接最小化する DeepRitz 法 という変種もあります。こちらは PDE の強形式ではなく、
$$ \Pi[\bm{u}] = \int_\Omega \frac{1}{2} \bm{\sigma} : \bm{\varepsilon} \, dV – \int_\Omega \bm{f} \cdot \bm{u} \, dV – \int_{\partial \Omega_t} \bm{t} \cdot \bm{u} \, dS $$
というポテンシャルエネルギーを直接最小化します。微分が1階で済むので学習が安定しやすい一方で、応力境界をやはり別途扱う必要があります。本記事では強形式の PINN を中心に進めますが、構造問題では DeepRitz もよく使われることを覚えておいてください。
熱と構造の方程式が揃いました。次に、PINN の学習を成功させるための工夫 — 境界条件の扱いとアーキテクチャ選択 — を見ていきます。
境界条件の扱い：ハード制約 vs ソフト制約
ソフト制約
これまで見てきたように、境界条件を損失項として加えるのがソフト制約です。実装が簡単で、Dirichlet・Neumann・Robin が同じ枠組みで書けます。
問題点は、損失項間のバランスです。$\mathcal{L}_b$ の重み $\lambda_b$ が小さいと境界が満たされず、大きすぎると内部の PDE 残差が無視される。実用上、$\lambda_b$ は問題に応じて手で調整するか、Wang ら（2021）の Neural Tangent Kernel（NTK）解析に基づく自動調整を使うのが主流です。
ハード制約
別のアプローチがハード制約で、ネットワーク出力を「境界条件を必ず満たす形」に組み立てます。たとえば1次元区間 $[0, L]$ で両端 $u(0) = u_0, u(L) = u_L$ を満たしたいなら、
$$ u_\theta(x, t) = u_0 (1 – x/L) + u_L (x/L) + x(L – x) N_\theta(x, t) $$
と書きます。$N_\theta$ はニューラルネット、$x(L-x)$ は両端でゼロになる「マスク関数」。これにより $u_\theta(0, t) = u_0$、$u_\theta(L, t) = u_L$ が構造的に保証されます。
ハード制約の利点は、$\mathcal{L}_b$ がそもそも不要になり、損失バランスの調整が一段楽になることです。欠点は、複雑な境界形状ではマスク関数の設計が難しいこと、Neumann/Robin 境界には素直に使えないことです。実用では「Dirichlet 境界はハード、Neumann/放射はソフト」のハイブリッドがよく選ばれます。
ここまでで損失関数の作り方が一通り見えました。次に、ネットワークの「中身」 — どんなアーキテクチャを選ぶかを考えます。
アーキテクチャ選択：FCNN・SIREN・Fourier features
FCNN（標準）
最も基本的なのは、全結合層 + tanh 活性化の素朴な多層パーセプトロン（FCNN）です。5〜10層、各層幅 50〜100 程度。tanh は滑らかで何度でも微分でき、PINN の PDE 残差計算と相性がよい。
ただし、FCNN にはスペクトルバイアスという弱点があります。Rahaman ら（2019）の指摘で、ニューラルネットは低周波関数を先に学習し、高周波成分は学習が遅い。温度や応力の分布が急峻なピークを持つ場合、これが命取りになります。
SIREN
Sitzmann ら（2020）の SIREN（Sinusoidal Representation Networks） は、活性化関数を $\sin(\omega_0 x)$ に置き換えることで高周波成分を捉えやすくしました。$\omega_0 \sim 30$ という大きな初期周波数で第1層を初期化するのがコツです。SIREN は微分も $\cos$ になって滑らかなままで、PINN との相性がよく、近年標準的に使われます。
Fourier features
Tancik ら（2020）の Fourier features は、入力を高次元の sin/cos 基底に射影してから FCNN に入れる手法です。
$$ \gamma(\bm{x}) = [\cos(2\pi \bm{B}\bm{x}), \sin(2\pi \bm{B}\bm{x})] $$
$\bm{B}$ は乱数行列。これでネットワークの実効的な周波数帯域を広げ、スペクトルバイアスを緩和します。位置エンコーディング（NeRFで有名な手法）と本質的に同じ発想です。
本記事の実装では、1D熱伝導には素朴な tanh FCNN、2D衛星パネルでは Fourier features を使い分けます。
アーキテクチャの選択肢が見えました。最後に、PINN の学習でしばしば起きる罠を共有しておきます。
訓練の難しさ：損失バランス・NTK・gradient pathology
損失項の重み
式 (2) には複数の損失項があり、それぞれスケールが違います。$\mathcal{L}_f$（PDE残差）は座標と物性によってスケールが決まる一方、$\mathcal{L}_b$（境界）は温度差オーダーで決まる。素朴に $\lambda_f = \lambda_b = 1$ とすると、片方が支配的になって他方が学習されない、という現象が頻発します。
実用的な対処は3段階あります。
- 無次元化（温度・長さ・時間を正規化する）
- 学習中に各項の勾配ノルムが揃うように $\lambda$ を動的に調整する（Wang 2021 の gradient normalization や NTK 重み付け）
- それでも難しいときは、curriculum learning（PDE 残差を初期はゼロ、徐々に強くする）を使う
NTK 解析
Neural Tangent Kernel（NTK）は、無限幅 NN の学習動力学が線形で記述できるという理論で、PINN の収束特性を解析するのに使われます。Wang ら（2022）は、PINN の各損失項の収束速度がスペクトル分布で決まることを示し、自動で $\lambda$ を調整する手法を提案しました。実装は重いですが、難しい問題ではこの自動調整が威力を発揮します。
Gradient pathology
PINN は時間発展問題で、初期条件付近では小さな誤差でも先の時刻で誤差が積み重なる、という性質があります（causal weighting がこの対策）。また、急峻な解の付近で残差が「歪んだ風景」を作り、Adam が局所最小に陥ることもあります。対策として L-BFGS で二次収束に頼る、SGD + 学習率スケジューリング で大域探索を促す、などが知られています。
理論と訓練のコツが揃いました。いよいよ実装に移ります。最初は1次元の熱伝導 — 解析解と比較して PINN が本当に PDE を解けているか確認しましょう。
1D 熱伝導 PINN の実装
問題設定
長さ $L = 1$ の棒に対し、
$$ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}, \quad x \in (0, 1), \; t \in (0, 0.5] $$
を解きます。初期条件 $T(x, 0) = \sin(\pi x)$、両端 Dirichlet $T(0, t) = T(1, t) = 0$、熱拡散率 $\alpha = 0.1$。この問題には次の解析解があります。
$$ T_{\text{exact}}(x, t) = \sin(\pi x) \, e^{-\alpha \pi^2 t} $$
これは「指数的に減衰しながら正弦波の形を保つ」典型的な熱伝導解で、PINN の精度を検証するベンチマークとして定番です。
PyTorch によるネットワーク定義
ネットワークは $(x, t) \to T$ の 2 入力 1 出力 FCNN とします。
import torch
import torch.nn as nn
class PINN1D(nn.Module):
"""1次元熱伝導用の素朴な FCNN（tanh活性化）"""
def __init__(self, layers=(2, 32, 32, 32, 32, 1)):
super().__init__()
modules = []
for i in range(len(layers) - 1):
modules.append(nn.Linear(layers[i], layers[i+1]))
if i < len(layers) - 2:
modules.append(nn.Tanh()) # 何度でも微分できる滑らかな活性化
self.net = nn.Sequential(*modules)
def forward(self, x, t):
# x, t は (N, 1) の Tensor を仮定
return self.net(torch.cat([x, t], dim=1))
ここで重要なのは、forward の中で torch.cat を使って 2 入力を結合していることです。後で x や t に対して autograd.grad を取りたいので、それぞれを requires_grad=True な独立 Tensor として保持する必要があります。
PDE 残差の autograd 計算
PDE 残差を計算する関数を実装します。
import torch
def pde_residual(model, x, t, alpha=0.1):
"""熱伝導方程式 ∂T/∂t - α ∂²T/∂x² = 0 の残差"""
T = model(x, t)
# ∂T/∂t を計算（grad_outputs に同形のテンソルを渡す）
T_t = torch.autograd.grad(
T, t, grad_outputs=torch.ones_like(T),
create_graph=True, retain_graph=True
)[0]
# ∂T/∂x
T_x = torch.autograd.grad(
T, x, grad_outputs=torch.ones_like(T),
create_graph=True, retain_graph=True
)[0]
# ∂²T/∂x²（さらにもう1回微分）
T_xx = torch.autograd.grad(
T_x, x, grad_outputs=torch.ones_like(T_x),
create_graph=True, retain_graph=True
)[0]
return T_t - alpha * T_xx
create_graph=True がポイントで、これによって「微分結果に対してさらに微分する」（高次微分）が可能になります。T_x を $x$ についてもう一度微分して $T_{xx}$ を得る、という流れが PINN の核心です。
コロケーション点と学習ループ
ランダムにサンプリングした内部点・境界点・初期点で、各損失を計算します。
import torch
import numpy as np
torch.manual_seed(0)
np.random.seed(0)
# デバイス
device = "cuda" if torch.cuda.is_available() else "cpu"
# モデル
model = PINN1D().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
# サンプル点を生成する関数
def sample_points(N_f=5000, N_b=200, N_0=200):
"""内部 N_f 点、境界 N_b 点、初期 N_0 点をランダムサンプル"""
# 内部点（PDE 残差用）
x_f = torch.rand(N_f, 1, device=device, requires_grad=True)
t_f = torch.rand(N_f, 1, device=device, requires_grad=True) * 0.5
# 境界点（x=0 と x=1 の両方を半々）
t_b = torch.rand(N_b, 1, device=device) * 0.5
x_b0 = torch.zeros_like(t_b)
x_b1 = torch.ones_like(t_b)
# 初期点（t=0）
x_0 = torch.rand(N_0, 1, device=device)
t_0 = torch.zeros_like(x_0)
return x_f, t_f, t_b, x_b0, x_b1, x_0, t_0
# 学習ループ
losses = []
for it in range(8000):
x_f, t_f, t_b, x_b0, x_b1, x_0, t_0 = sample_points()
# PDE 残差ロス
r = pde_residual(model, x_f, t_f)
loss_f = (r ** 2).mean()
# 境界ロス（両端 0）
T_b0 = model(x_b0, t_b)
T_b1 = model(x_b1, t_b)
loss_b = (T_b0 ** 2).mean() + (T_b1 ** 2).mean()
# 初期ロス（T(x,0) = sin(πx)）
T_0 = model(x_0, t_0)
T_0_true = torch.sin(np.pi * x_0)
loss_0 = ((T_0 - T_0_true) ** 2).mean()
loss = loss_f + 10.0 * loss_b + 10.0 * loss_0 # 境界・初期を強めに
optimizer.zero_grad()
loss.backward()
optimizer.step()
losses.append(loss.item())
if it % 1000 == 0:
print(f"iter {it:5d} loss={loss.item():.3e} "
f"(f={loss_f.item():.2e}, b={loss_b.item():.2e}, 0={loss_0.item():.2e})")
ここでの設計判断を補足します。loss_f の重みは 1、loss_b と loss_0 の重みは 10 にしました。これは、初期・境界が支配的に学ばれてから PDE 残差が効くようにする、簡単な curriculum 的な工夫です。問題によっては動的に調整する手法（GradNorm, NTK）に置き換える方が安定します。
解析解との比較
学習が終わったら、$(x, t)$ メッシュ上で予測値と解析解を比較します。
import numpy as np
import torch
import matplotlib.pyplot as plt
# 比較用のメッシュ
nx, nt = 100, 5
x_test = np.linspace(0, 1, nx)
t_snap = [0.0, 0.1, 0.2, 0.3, 0.4]
alpha = 0.1
plt.figure(figsize=(10, 6))
for i, ts in enumerate(t_snap):
# PINN 予測
xt = torch.tensor(x_test.reshape(-1, 1), dtype=torch.float32, device=device)
tt = torch.full_like(xt, ts)
with torch.no_grad():
T_pinn = model(xt, tt).cpu().numpy().flatten()
# 解析解
T_exact = np.sin(np.pi * x_test) * np.exp(-alpha * np.pi ** 2 * ts)
plt.plot(x_test, T_exact, '-', color=f'C{i}', label=f't={ts:.2f} exact')
plt.plot(x_test, T_pinn, 'o', color=f'C{i}', mfc='none', ms=4)
plt.xlabel('x')
plt.ylabel('T(x, t)')
plt.title('1D heat equation: PINN (circles) vs exact (lines)')
plt.legend(ncol=2)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('pinn_1d_heat.png', dpi=150, bbox_inches='tight')
plt.show()
# 最大誤差
xt = torch.tensor(np.tile(x_test, len(t_snap)).reshape(-1, 1), dtype=torch.float32, device=device)
tt = torch.tensor(np.repeat(t_snap, nx).reshape(-1, 1), dtype=torch.float32, device=device)
with torch.no_grad():
T_pinn_all = model(xt, tt).cpu().numpy().flatten()
T_exact_all = np.sin(np.pi * np.tile(x_test, len(t_snap))) * \
np.exp(-alpha * np.pi ** 2 * np.repeat(t_snap, nx))
print(f"最大絶対誤差: {np.max(np.abs(T_pinn_all - T_exact_all)):.3e}")
このコードを実行すると、PINN の予測（丸印）は解析解（実線）とほぼ重なります。最大絶対誤差はおおよそ $10^{-3}$ オーダー、相対誤差にして約 0.1% 程度に収束するはずです。重要なのは2点 — まず、初期の正弦波形状が時間とともに指数的に減衰していく解析解の振る舞いを PINN が正しく学習していること。これは、座標 $(x, t)$ について自動微分で取った高次微分が、本物の物理的拡散方程式をきちんと表現できていることを意味します。次に、両端の Dirichlet 境界 $T = 0$ も損失項を通じて十分満たされていることが、$x = 0, 1$ で全曲線がゼロに揃うことから確認できます。FEM 並みの精度には届かないこともありますが、メッシュレスでこのレベルが出るのは、簡単な実装としては悪くない成果です。
ここまでで PINN が本当に PDE を解けることを確認できました。次に、より実用に近い 2 次元衛星パネル問題に進みます。
2D 衛星パネル熱問題
問題設定
衛星のラジエータパネルを模した 2 次元矩形領域 $\Omega = [0, 1] \times [0, 1]$ を考えます。
$$ \rho c \frac{\partial T}{\partial t} = k \nabla^2 T + Q(x, y, t) $$
物性値はアルミ合金相当で、無次元化した形で $\alpha = k/(\rho c) = 1$ とします。発熱源 $Q$ は中央付近の電子機器を模して
$$ Q(x, y, t) = Q_0 \exp\left(-50\left((x – 0.5)^2 + (y – 0.5)^2\right)\right) \cdot (1 + 0.5 \sin(2\pi t)) $$
と置きます。$\sin(2\pi t)$ がパルス変動（軌道変動による発熱変動）を模しています。境界は4辺すべて $T = 0$ の Dirichlet（簡略化）、初期は $T(x, y, 0) = 0$。
Fourier features 付きネットワーク
スペクトルバイアスを緩和するため、入力に Fourier features を噛ませます。
import torch
import torch.nn as nn
import numpy as np
class FourierFeatures(nn.Module):
"""ガウス乱数行列で入力を高周波特徴へ射影"""
def __init__(self, in_dim, num_features=32, sigma=2.0, seed=0):
super().__init__()
rng = torch.Generator().manual_seed(seed)
# 訓練しない固定の乱数行列
B = sigma * torch.randn(in_dim, num_features, generator=rng)
self.register_buffer('B', B)
def forward(self, x):
# x: (N, in_dim) -> (N, 2*num_features)
proj = 2 * np.pi * x @ self.B
return torch.cat([torch.sin(proj), torch.cos(proj)], dim=1)
class PINN2D(nn.Module):
def __init__(self, num_features=32, hidden=64, depth=4):
super().__init__()
self.ff = FourierFeatures(3, num_features=num_features, sigma=2.0)
layers = [nn.Linear(2 * num_features, hidden), nn.Tanh()]
for _ in range(depth - 1):
layers += [nn.Linear(hidden, hidden), nn.Tanh()]
layers += [nn.Linear(hidden, 1)]
self.net = nn.Sequential(*layers)
def forward(self, x, y, t):
xyt = torch.cat([x, y, t], dim=1)
return self.net(self.ff(xyt))
FourierFeatures で 3 次元入力 $(x, y, t)$ をガウス乱数行列 $\bm{B}$ で射影し、$\sin/\cos$ 経由で高周波基底を作っています。register_buffer で $\bm{B}$ を保持することで、学習時には固定（パラメータ化しない）にしています。
2D の PDE 残差
2 階偏微分が $x, y$ の2方向ぶん必要になります。
import torch
def pde_residual_2d(model, x, y, t, Q_fn, alpha=1.0):
"""∂T/∂t - α(∂²T/∂x² + ∂²T/∂y²) - Q = 0"""
T = model(x, y, t)
T_t = torch.autograd.grad(T, t, torch.ones_like(T),
create_graph=True, retain_graph=True)[0]
T_x = torch.autograd.grad(T, x, torch.ones_like(T),
create_graph=True, retain_graph=True)[0]
T_y = torch.autograd.grad(T, y, torch.ones_like(T),
create_graph=True, retain_graph=True)[0]
T_xx = torch.autograd.grad(T_x, x, torch.ones_like(T_x),
create_graph=True, retain_graph=True)[0]
T_yy = torch.autograd.grad(T_y, y, torch.ones_like(T_y),
create_graph=True, retain_graph=True)[0]
Q = Q_fn(x, y, t)
return T_t - alpha * (T_xx + T_yy) - Q
def heat_source(x, y, t, Q0=5.0):
return Q0 * torch.exp(-50 * ((x - 0.5)**2 + (y - 0.5)**2)) * \
(1.0 + 0.5 * torch.sin(2 * np.pi * t))
T_x と T_y をそれぞれ作り、もう一度 $x, y$ で微分してラプラシアンを組み立てます。3階以上にはなりませんが、autograd の合成深さが深くなる分メモリは食います。
学習ループ（一部抜粋）
import torch
import numpy as np
torch.manual_seed(0)
device = "cuda" if torch.cuda.is_available() else "cpu"
model = PINN2D().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
T_end = 0.3
def sample_2d(N_f=4000, N_b=400, N_0=400):
x_f = torch.rand(N_f, 1, device=device, requires_grad=True)
y_f = torch.rand(N_f, 1, device=device, requires_grad=True)
t_f = torch.rand(N_f, 1, device=device, requires_grad=True) * T_end
# 境界（4辺）
s = torch.rand(N_b, 1, device=device)
tb = torch.rand(N_b, 1, device=device) * T_end
edges = [
(torch.zeros_like(s), s, tb), # x=0
(torch.ones_like(s), s, tb), # x=1
(s, torch.zeros_like(s), tb), # y=0
(s, torch.ones_like(s), tb), # y=1
]
# 初期
x0 = torch.rand(N_0, 1, device=device)
y0 = torch.rand(N_0, 1, device=device)
t0 = torch.zeros_like(x0)
return x_f, y_f, t_f, edges, x0, y0, t0
for it in range(10000):
x_f, y_f, t_f, edges, x0, y0, t0 = sample_2d()
# PDE 残差
r = pde_residual_2d(model, x_f, y_f, t_f, heat_source)
loss_f = (r ** 2).mean()
# 境界
loss_b = 0.0
for (xb, yb, tb) in edges:
Tb = model(xb, yb, tb)
loss_b = loss_b + (Tb ** 2).mean()
# 初期
T0 = model(x0, y0, t0)
loss_0 = (T0 ** 2).mean()
loss = loss_f + 50.0 * loss_b + 50.0 * loss_0
optimizer.zero_grad()
loss.backward()
optimizer.step()
if it % 1000 == 0:
print(f"iter {it:5d} loss={loss.item():.3e}")
ここで境界・初期の重みを 50 と大きめにしています。これは、初期と境界が両方ともゼロという「真っ平らな解」が trivial に最小化されないよう、強い拘束として効かせるためです。長時間学習させたい場合は、まず loss_b, loss_0 を強くしてから徐々に下げる curriculum も有効です。
可視化と FEM との比較準備
学習後、$(x, y)$ 平面のスナップショットを描きます。
import torch
import numpy as np
import matplotlib.pyplot as plt
n = 50
xs = torch.linspace(0, 1, n, device=device).reshape(-1, 1)
ys = torch.linspace(0, 1, n, device=device).reshape(-1, 1)
X, Y = torch.meshgrid(xs.squeeze(), ys.squeeze(), indexing='ij')
X_f = X.reshape(-1, 1); Y_f = Y.reshape(-1, 1)
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
for ax, ts in zip(axes, [0.05, 0.15, 0.30]):
T_in = torch.full_like(X_f, ts)
with torch.no_grad():
T = model(X_f, Y_f, T_in).cpu().numpy().reshape(n, n)
cs = ax.imshow(T.T, origin='lower', extent=[0, 1, 0, 1],
cmap='hot', vmin=0, vmax=T.max())
ax.set_title(f't = {ts:.2f}')
ax.set_xlabel('x'); ax.set_ylabel('y')
plt.colorbar(cs, ax=ax)
plt.suptitle('2D satellite panel temperature (PINN)')
plt.tight_layout()
plt.savefig('pinn_2d_panel.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、衛星パネル中央の発熱源を中心に温度がガウス状に広がっていき、時刻が進むと境界の冷却（$T=0$）と発熱源のパルス変動の競合で温度分布が時間変化する様子が読み取れます。$t=0.05$ では発熱が局所的なピークを作り、$t=0.30$ では境界冷却によって全体的に温度が下がりつつ、中央の振動的な発熱でパルス的な温度変動が現れます。Fourier features がなければ、こうした急峻な中央ピークは tanh FCNN では十分に表現できず、ぼんやりした分布になりがちでした。スペクトルバイアスの緩和が効いていることが、可視化からはっきり見て取れます。
ここまでで 2D PINN が動くことを確認できました。気になるのは「で、FEM と比べてどうなのか」です。次にこれを正面から扱います。
FEM との計算時間比較
比較の枠組み
公平な比較は難しい — 用語を整理しておきます。FEM は「メッシュを切ってから1ステップずつ時間積分」する離散解法で、メッシュ依存・時間ステップ依存です。PINN は「学習時間が前払い」で、学習後の評価は座標を入力するだけで即座に得られます。よって、
- オフライン段階: メッシュ生成 + FEM の数値積分 vs PINN の学習
- オンライン段階: 既存メッシュでの追加時刻評価 vs 学習済み PINN の評価
の2軸で見るのが正しい比較です。
簡易 FEM の実装と計時
scipy で 2D 熱伝導の有限差分（簡略 FEM 相当）を組み、PINN と所要時間を比べます。
import numpy as np
import time
from scipy.sparse import diags, eye, kron
from scipy.sparse.linalg import spsolve
def fdm_2d_heat(N=50, T_end=0.3, dt=1e-3, Q0=5.0, alpha=1.0):
"""2次元熱伝導の陰的オイラー（FEMに近い精度の有限差分）"""
h = 1.0 / (N - 1)
# 1Dラプラシアン
e = np.ones(N)
L1 = diags([e, -2*e, e], [-1, 0, 1], shape=(N, N)).tolil()
# Dirichlet 境界（端の行をゼロに）
L1[0, :] = 0; L1[-1, :] = 0
L1 = L1.tocsr()
I = eye(N).tocsr()
L2 = (kron(I, L1) + kron(L1, I)) / h**2
A = eye(N*N) - dt * alpha * L2
# 初期・座標
x = np.linspace(0, 1, N)
X, Y = np.meshgrid(x, x, indexing='ij')
u = np.zeros(N * N)
n_steps = int(T_end / dt)
snapshots = {}
save_t = [0.05, 0.15, 0.30]
for n in range(n_steps):
t = (n + 1) * dt
Q = Q0 * np.exp(-50 * ((X - 0.5)**2 + (Y - 0.5)**2)) * \
(1.0 + 0.5 * np.sin(2 * np.pi * t))
b = u + dt * Q.flatten()
# 境界をゼロに固定
Qmask = np.ones((N, N), dtype=bool)
Qmask[0,:] = Qmask[-1,:] = Qmask[:,0] = Qmask[:,-1] = False
u = spsolve(A, b)
u = u.reshape(N, N)
u[0,:] = u[-1,:] = u[:,0] = u[:,-1] = 0
u = u.flatten()
for ts in save_t:
if abs(t - ts) < dt / 2:
snapshots[ts] = u.reshape(N, N).copy()
return snapshots
# 計時
t0 = time.perf_counter()
fem_snaps = fdm_2d_heat(N=50)
t_fem = time.perf_counter() - t0
print(f"FEM(50x50) 時間: {t_fem:.2f} s")
# PINN は上で学習済み（10000 iter）と仮定 — 学習時間は別途記録
# 学習: t_pinn_train（記録済み）
# 評価: 50x50 を1回評価する時間
n = 50
xs = torch.linspace(0, 1, n, device=device).reshape(-1, 1)
ys = torch.linspace(0, 1, n, device=device).reshape(-1, 1)
X, Y = torch.meshgrid(xs.squeeze(), ys.squeeze(), indexing='ij')
X_f = X.reshape(-1, 1); Y_f = Y.reshape(-1, 1)
T_in = torch.full_like(X_f, 0.15)
t0 = time.perf_counter()
for _ in range(10):
with torch.no_grad():
_ = model(X_f, Y_f, T_in)
t_pinn_eval = (time.perf_counter() - t0) / 10
print(f"PINN 評価 1スナップショット: {t_pinn_eval*1000:.1f} ms")
このコードを実行すると、典型的な数値感は次のとおりです（CPU 環境、N=50、T_end=0.3）。FEM は数秒〜十数秒で全時間ステップを解き終わります。PINN は学習に数十秒〜数分かかる一方で、学習後の評価は 1〜10 ms オーダーです。つまり、
- 1回だけ解けばよい問題: FEM が圧勝。PINN は学習コストが回収できない
- 同じ問題を何度も評価したい / 任意点で値が欲しい: PINN が有利。FEM はその都度メッシュ補間が必要
- 形状やパラメータが変わるたびに解き直したい: ハイパーネットを使った PINN（パラメータ条件付き）が現実的
つまり PINN は「サロゲートモデルとしての使い方」が本命で、FEM を完全置換する手法ではありません。
精度の比較
スナップショット $t = 0.15$ で PINN と FEM の解の差を見てみます。
import torch
import numpy as np
import matplotlib.pyplot as plt
n = 50
xs = torch.linspace(0, 1, n, device=device).reshape(-1, 1)
ys = torch.linspace(0, 1, n, device=device).reshape(-1, 1)
X, Y = torch.meshgrid(xs.squeeze(), ys.squeeze(), indexing='ij')
X_f = X.reshape(-1, 1); Y_f = Y.reshape(-1, 1)
T_in = torch.full_like(X_f, 0.15)
with torch.no_grad():
T_pinn = model(X_f, Y_f, T_in).cpu().numpy().reshape(n, n)
T_fem = fem_snaps[0.15]
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
axes[0].imshow(T_pinn.T, origin='lower', cmap='hot'); axes[0].set_title('PINN')
axes[1].imshow(T_fem.T, origin='lower', cmap='hot'); axes[1].set_title('FEM')
diff = np.abs(T_pinn - T_fem)
cs = axes[2].imshow(diff.T, origin='lower', cmap='viridis')
axes[2].set_title(f'|PINN - FEM| max={diff.max():.3e}')
plt.colorbar(cs, ax=axes[2])
plt.tight_layout()
plt.savefig('pinn_vs_fem.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"L2 相対誤差: {np.linalg.norm(T_pinn - T_fem) / np.linalg.norm(T_fem):.3e}")
このグラフから、PINN と FEM の温度場は定性的にはほぼ同じガウス状ピークと境界冷却を再現しますが、差分マップ（右パネル）に細かい誤差が残ることがわかります。L2 相対誤差はおよそ数 % のオーダーになるはずです（学習回数や Fourier features の分散 $\sigma$ に依存）。誤差の局在は境界近くと発熱源中央付近に多く、これは PINN が「平均的に滑らかな解」を好む（スペクトルバイアスの残滓）一方で、FEM は局所的な不連続を素直に表現する、という両者の性格の違いから来ています。設計初期のサロゲートとしては数 % 誤差で十分なケースが多く、PINN の使い所はそうした「速さを精度に優先する」場面です。
ここまでで、PINN がどんな問題にどう効くかが具体的な数値で見えました。最後に、これらの結果を踏まえた実装上のヒントと、実機応用での位置づけをまとめます。
実機応用への展望
衛星バス熱予測のサロゲート
衛星バスの熱モデルは、設計の各段階で「太陽光入射角」「ヒーター ON/OFF」「電力消費」など多くのシナリオで解き直す必要があります。FEMで全シナリオを毎回解くと計算量が膨大になります。ここで PINN をパラメータ依存サロゲートとして学習させると、入力に運転条件を加えた $(x, y, t, \theta_{\text{ops}})$ を取って温度を返す高速な代替モデルになります。一度学習しておけば、シナリオの違いはネットワーク評価1回分（ミリ秒オーダー）で得られます。
軌道上のリアルタイム熱応力監視
熱電対は衛星に十数点〜数十点しか付けられず、観測点は疎です。観測点での測定値を $\mathcal{L}_d$ として、PDE 残差 $\mathcal{L}_f$ と同時に学習すると、観測されていない位置の温度を物理拘束付きで推定できます。これは「データ同化」の機械学習版と言えます。重要なのは、観測ノイズが多くてもデータだけでなく PDE が解を支えてくれる点で、純粋なニューラルネットの内挿よりも安定です。
モデル更新による軌道上適応
軌道上で熱電対の値が予測とずれた場合、熱伝導率 $k$ やヒートパイプの熱抵抗が劣化している可能性があります。PINN では未知パラメータも学習対象に加え、
$$ \mathcal{L}(\theta, k) = \mathcal{L}_f(\theta, k) + \mathcal{L}_d(\theta) + \cdots $$
として観測残差と PDE 残差を同時に最小化することで、$k$ の値を逆推定できます。FEMで同じことをやるには随伴法など別実装が必要で重いのですが、PINN は枠組みが同じで自然に逆問題に拡張できる、というのが大きな利点です。
既存 FEM との統合
実用では FEM をベースに、特に時間がかかる非線形項（放射境界）や逆問題部分だけを PINN に置き換える、というハイブリッドが現実的です。たとえば
- 線形熱伝導は FEM で陰解法
- 放射境界の非線形項のみ PINN で表現し、FEM の境界条件として渡す
- パラメータ逆推定は PINN で、解析自体は FEM
といった連携です。両者は競合ではなく補完関係にある、というのが本記事を通じての結論です。
まとめ
本記事では、PINN を熱伝導・構造方程式に拡張する枠組みを、Raissi 2019 の損失関数の構成から実装まで丁寧に辿りました。
- PINN は PDE の連続解を学習する — ネットワーク $u_\theta(\bm{x}, t)$ を PDE 残差・境界・初期・観測の各損失で訓練し、メッシュレスに $u$ を得る
- 熱伝導方程式 $\rho c \partial T/\partial t = \nabla \cdot (k \nabla T) + Q$ は PDE 残差を autograd で計算し、境界条件はソフト/ハードで扱える
- 構造方程式 $\nabla \cdot \bm{\sigma} = \bm{f}$ は変位場 $\bm{u}_\theta$ から微分でひずみ・応力を組み立て、熱応力は $\bm{\varepsilon}_{\text{th}} = \alpha_T \Delta T \bm{I}$ を引いて表現
- アーキテクチャの選択肢 として FCNN（標準）、SIREN（高周波対応）、Fourier features（スペクトルバイアス緩和）があり、問題に応じて使い分ける
- 訓練の難しさ は損失バランス、スペクトルバイアス、causal weighting に集約され、無次元化と curriculum learning が定石
- FEMとの比較: 1回解く問題では FEM、繰り返し評価・任意点評価・逆問題ではPINN。両者は補完関係
数値実験では、1D 熱伝導で解析解と $10^{-3}$ オーダーで一致、2D 衛星パネルで Fourier features を入れた PINN が中央発熱源と境界冷却の競合を捉えること、FEM とは数 % 誤差で定性的に一致しつつ、評価コストは桁違いに軽いこと、を確認しました。
PINN は「次世代の数値解法」というよりは、「物理を取り込んだ機械学習」と理解するのが正確です。サロゲートモデル、逆問題、データ同化 — 既存の数値解法では重かった用途に新しい道を開く道具として、宇宙機の熱構造解析にも有望な選択肢になります。
次のステップとして、以下の記事も参考にしてください。