高周波エンジニアのためのAI・機械学習入門（GPU編6）PythonとKeras3.0を使って3次のLCバンドパスフィルタ（BPF）のSパラメータから物理法則を取り入れたPINN(Physics Informed Neural Networks)で素子の値（L、C）推定する。使ったのはLC回路網をSパラメータに直して入力と比較。

今回はPINN。

Physics-Informed Neural Networks (PINNs)

といっても今回は回路なので、微分方程式でなくてLCの値を回路網理論でSパラメータに直したものを物理法則とする。

Google Gemini 3 Proと相談しながらすぐできた。

コードはこちら。

import os

os.environ["KERAS_BACKEND"] = "jax"  # JAXバックエンド

import keras

from keras import ops

import numpy as np

import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split

from sklearn.preprocessing import StandardScaler

plt.rcParams['font.family'] = 'Noto Sans CJK JP'

# ==========================

# データ読み込み

# ==========================

data_label = np.load("data_label.npz")

data = data_label["data"]

label = data_label["label"]

x_train, x_test, y_train, y_test = train_test_split(data, label, test_size=0.3, random_state=0)

# ==========================

# ラベルスケーリング

# ==========================

scaler_y = StandardScaler()

y_train_f = scaler_y.fit_transform(y_train)

y_test_f  = scaler_y.transform(y_test)

# PINN内で逆変換するために平均と標準偏差をTensorとして保持

# float32型にキャストしておく

y_mean = keras.ops.convert_to_tensor(scaler_y.mean_, dtype="float32")

y_scale = keras.ops.convert_to_tensor(scaler_y.scale_, dtype="float32")

# ==========================

# 物理法則 (Circuit Theory)

# ==========================

def calculate_s_parameters_from_LC(freq_norm, lc_params):

    """

    freq_norm: 正規化周波数 (batch, 200)

    lc_params: 物理単位のLC値 (batch, 6) -> [L1, C1, L2, C2, L3, C3]

    """

    f_stop = 20.0e9 # 20 GHz

    z0 = 50.0 + 0j

    # 周波数を物理単位(Hz)に戻して角周波数計算

    freq_hz = freq_norm * f_stop

    omega = 2.0 * np.pi * freq_hz

    omega = ops.cast(omega, "complex64")

    epsilon = 1e-9 # 0除算防止

    # (batch, 1) に拡張してブロードキャスト

    L1 = ops.expand_dims(lc_params[:, 0], axis=1) * 1e-9

    C1 = ops.expand_dims(lc_params[:, 1], axis=1) * 1e-12

    L2 = ops.expand_dims(lc_params[:, 2], axis=1) * 1e-9

    C2 = ops.expand_dims(lc_params[:, 3], axis=1) * 1e-12

    L3 = ops.expand_dims(lc_params[:, 4], axis=1) * 1e-9

    C3 = ops.expand_dims(lc_params[:, 5], axis=1) * 1e-12

    # Series (L1, C1)

    Zs1 = 1j * omega * L1 + 1.0 / (1j * omega * C1 + epsilon)

    # Shunt (L2, C2)

    Yp2 = 1j * omega * C2 + 1.0 / (1j * omega * L2 + epsilon)

    # Series (L3, C3)

    Zs3 = 1j * omega * L3 + 1.0 / (1j * omega * C3 + epsilon)

    # ABCD行列 (A, B, C, D) の要素計算

    A = 1.0 + Zs1 * Yp2

    B = (1.0 + Zs1 * Yp2) * Zs3 + Zs1

    C = Yp2

    D = Yp2 * Zs3 + 1.0

    # Sパラメータ変換

    delta = A + B/z0 + C*z0 + D

    S11 = (A + B/z0 - C*z0 - D) / delta

    S21 = 2.0 / delta

    return S11, S21

# ==========================

# PINN モデル定義 (compute_lossのみオーバーライド)

# ==========================

class CircuitPINN(keras.Model):

    def __init__(self, inputs, outputs, y_mean, y_scale, physics_weight=0.1, **kwargs):

        super().__init__(inputs=inputs, outputs=outputs, **kwargs)

        self.y_mean = y_mean

        self.y_scale = y_scale

        self.physics_weight = physics_weight

        # Loss監視用メトリクス

        self.data_loss_tracker = keras.metrics.Mean(name="data_loss")

        self.phys_loss_tracker = keras.metrics.Mean(name="phys_loss")

    # Keras 3 ではここをオーバーライドするのが推奨

    def compute_loss(self, x=None, y=None, y_pred=None, sample_weight=None):

        # 1. Data Loss (教師データとの誤差)

        # y: 教師データ(正規化済み), y_pred: 予測値(正規化済み)

        data_loss = ops.mean(ops.square(y - y_pred))

        # 2. Physics Loss (物理法則との整合性)

        # 正規化された予測値を物理単位(nH, pF)に戻す

        y_pred_physical = y_pred * self.y_scale + self.y_mean

        # L, Cは正の値であるべきなので絶対値をとる（またはReLU）

        y_pred_physical = ops.abs(y_pred_physical)

        # 入力 x から周波数と正解Sパラメータを抽出

        # x shape: (batch, 200, 5) -> [freq, S11r, S11i, S21r, S21i]

        freq_norm = x[:, :, 0]

        s11_true_real = x[:, :, 1]

        s11_true_imag = x[:, :, 2]

        s21_true_real = x[:, :, 3]

        s21_true_imag = x[:, :, 4]

        # 予測されたLC値からSパラメータを計算

        s11_pred, s21_pred = calculate_s_parameters_from_LC(freq_norm, y_pred_physical)

        # Sパラメータの誤差を計算 (実部・虚部)

        phys_loss = (

            ops.mean(ops.square(s11_true_real - ops.real(s11_pred))) +

            ops.mean(ops.square(s11_true_imag - ops.imag(s11_pred))) +

            ops.mean(ops.square(s21_true_real - ops.real(s21_pred))) +

            ops.mean(ops.square(s21_true_imag - ops.imag(s21_pred)))

        )

        total_loss = data_loss + self.physics_weight * phys_loss

        # メトリクスの更新

        self.data_loss_tracker.update_state(data_loss)

        self.phys_loss_tracker.update_state(phys_loss)

        return total_loss

    @property

    def metrics(self):

        # fit() のプログレスバーに表示されるメトリクス一覧

        return [self.data_loss_tracker, self.phys_loss_tracker]

# ==========================

# モデル構築

# ==========================

keras.utils.set_random_seed(1)

inputs = keras.Input(shape=(200, 5))

x = keras.layers.Flatten()(inputs)

x = keras.layers.Dense(256, activation="relu")(x)

x = keras.layers.BatchNormalization()(x)

x = keras.layers.Dense(256, activation="relu")(x)

x = keras.layers.BatchNormalization()(x)

#x = keras.layers.Dropout(0.1)(x)

outputs = keras.layers.Dense(6)(x)

# モデル初期化

model = CircuitPINN(

    inputs=inputs,

    outputs=outputs,

    y_mean=y_mean,

    y_scale=y_scale,

    physics_weight=0.1  # 重みは適宜調整してください

)

# コンパイル (loss引数はcompute_lossで定義しているので不要)

model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.001))

batch_size = 32

epochs = 300

# 学習

history = model.fit(

    x_train,

    y_train_f,

    batch_size=batch_size,

    epochs=epochs,

    validation_split=0.15,

)

# ==========================

# 評価・可視化

# ==========================

# 推論

y_pred_f = model.predict(x_test)

y_pred = scaler_y.inverse_transform(y_pred_f)

# R2スコア

metric = keras.metrics.R2Score()

metric.update_state(y_test, y_pred)

print("\nR2 Score:", metric.result())

# 誤差率計算

eps = 1e-12

valid_mask = np.abs(y_test) > eps

pct_error = np.zeros_like(y_test)

pct_error[valid_mask] = np.abs(

    (y_test[valid_mask] - y_pred[valid_mask]) / y_test[valid_mask] * 100

)

print("\n% Error per element (mean):", pct_error.mean(axis=0))

# プロット

row, column = 2, 3

legend = ["L1", "C1", "L2", "C2", "L3", "C3"]

fig, ax = plt.subplots(2, 3, figsize=(15,9))

for i in range(row):

    for j in range(column):

        count = column * i + j

        maxvalue = max(y_pred[:, count].max(), y_test[:,count].max())

        ax[i,j].scatter(y_pred[:, count], y_test[:,count], c="r", s=5)

        ax[i,j].plot([0,maxvalue], [0,maxvalue], "--", c="black")

        ax[i,j].set_xlabel("推定した値")

        ax[i,j].set_ylabel("実際の値")

        ax[i,j].set_xlim(0, maxvalue)

        ax[i,j].set_ylim(0, maxvalue)

        ax[i,j].grid()

        ax[i,j].legend([legend[count] + f" 平均誤差{pct_error.mean(axis=0)[count]:.2f}%"])

fig.tight_layout()

plt.show()

結果はこちら。うーん、単純なDNNより悪くなってる…こんな単純な問題だと逆効果か。