懸架系の高精度モーダル解析#
運用モーダル解析（OMA: Operational Modal Analysis） は、既知の 励振力なしに振動データから共振周波数、減衰比、モード形状を抽出する手法です。 KAGRA の懸架系（TypeA、TypeB、TypeBp）の固有振動モード同定に使われています。
このチュートリアルで学ぶこと：
多自由度 FRF（周波数応答関数）行列の合成
pyOMA を使った OMA 実行（未インストール時は解析的フォールバック）
FrequencySeries.from_pyoma_results()で gwexpy に変換するモード形状・減衰比の可視化と FRF 再構成
注意: このノートブックは自己完結しており、すべてのセルは pyOMA なしで実行できます。実際の OMA 解析を行うには
pip install pyOMAでインストールしてください。
case_violin_mode.ipynb との違い： バイオリンモードチュートリアルは 単一チャネルのスペクトルピークフィッティングに特化しています。 本チュートリアルは多チャネル FRF 行列と完全なモード形状抽出を扱い、 作動行列較正（MODAL2COIL マトリクス）や能動制振フィルタ設計の基礎になります。
準備#
import warnings
import matplotlib.pyplot as plt
import numpy as np
from gwexpy.frequencyseries import FrequencySeriesMatrix
from gwexpy.interop._modal_helpers import build_frf_matrix, build_mode_dataframe
1. 合成多自由度 FRF 行列#
3 自由度チェーン懸架系（KAGRA TypeB 中間質量ステージの簡略モデル）を シミュレートします。FRF 行列 H[i,j](f) は DOF j への力に対する DOF i の応答を表します。
# --- System parameters ---
fs = 100.0 # sampling rate [Hz]
f_modes = [0.38, 0.55, 0.72] # resonance frequencies [Hz] (pendulum chain)
Q_modes = [150.0, 80.0, 60.0] # quality factors
n_dof = 3
freqs = np.linspace(0.05, 2.0, 2000) # 0.05 – 2 Hz
omega = 2 * np.pi * freqs
def single_dof_frf(f, f0, Q):
# Complex FRF of a SDOF resonator: H(f) = 1 / (1 - (f/f0)^2 + i*f/(f0*Q))
r = f / f0
return 1.0 / (1 - r**2 + 1j * r / Q)
# Mode shapes (columns = modes, rows = DOFs)
# Physically: mode 1 = in-phase, mode 2 = anti-phase middle, mode 3 = anti-phase all
Phi = np.array([
[ 1.00, 1.00, 1.00], # DOF 0 (upper stage)
[ 0.62, -0.38, 0.85], # DOF 1 (intermediate stage)
[ 0.35, -0.72, -1.20], # DOF 2 (mirror)
])
# Build FRF matrix H[resp, ref, freq] = sum_r Phi[resp,r]*Phi[ref,r]*H_r(f)
n_freq = len(freqs)
frf_data = np.zeros((n_dof, n_dof, n_freq), dtype=complex)
for r, (f0, Q) in enumerate(zip(f_modes, Q_modes)):
Hr = single_dof_frf(freqs, f0, Q)
for i in range(n_dof):
for j in range(n_dof):
frf_data[i, j, :] += Phi[i, r] * Phi[j, r] * Hr
# Add measurement noise
rng = np.random.default_rng(0)
noise_level = 0.02 * np.max(np.abs(frf_data))
frf_data += noise_level * (rng.standard_normal(frf_data.shape)
+ 1j * rng.standard_normal(frf_data.shape))
print(f"FRF matrix shape: {frf_data.shape} [n_resp × n_ref × n_freq]")
print(f"Modes: {f_modes} Hz, Q = {Q_modes}")
FRF matrix shape: (3, 3, 2000) [n_resp × n_ref × n_freq]
Modes: [0.38, 0.55, 0.72] Hz, Q = [150.0, 80.0, 60.0]
2. gwexpy FrequencySeriesMatrix への変換#
OMA を実行する前でも、build_frf_matrix() を使って FRF を gwexpy に取り込めます。これは from_pyoma_results() が内部で使うものと同じヘルパーです。
dof_labels = [f"DOF_{i}" for i in range(n_dof)]
frf_matrix = build_frf_matrix(
FrequencySeriesMatrix,
freqs,
frf_data,
response_names=dof_labels,
reference_names=dof_labels,
unit="m/N",
name="TypeB IM stage",
)
print(type(frf_matrix))
print("Shape:", frf_matrix.shape)
print("Frequency bins:", frf_matrix.shape[-1])
# --- Plot the driving-point FRF (DOF_0 → DOF_0) ---
h00 = frf_data[0, 0, :]
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 5), sharex=True)
ax1.semilogy(freqs, np.abs(h00), color="steelblue", lw=1.5)
ax1.set_ylabel("|H₀₀| [m/N]")
ax1.set_title("Driving-point FRF — DOF 0 (upper stage)")
ax1.grid(True, alpha=0.4)
for f0 in f_modes:
ax1.axvline(f0, color="red", ls="--", alpha=0.5, lw=0.8)
ax2.plot(freqs, np.angle(h00, deg=True), color="darkorange", lw=1.5)
ax2.set_ylabel("Phase [deg]")
ax2.set_xlabel("Frequency [Hz]")
ax2.set_ylim(-200, 200)
ax2.grid(True, alpha=0.4)
plt.tight_layout()
plt.show()
<class 'gwexpy.frequencyseries.matrix.FrequencySeriesMatrix'>
Shape: (3, 3, 2000)
Frequency bins: 2000
3. pyOMA による運用モーダル解析#
pyOMA の 周波数領域分解 (FDD) アルゴリズムを使用します。続いて from_pyoma_results() コンバータが、抽出したモーダルパラメータを pandas.DataFrame（パラメータ要約）または FrequencySeriesMatrix（再構成したモード形状）としてまとめます。
PYOMA_AVAILABLE = False
try:
from pyOMA.core.OMA import ModeFDD
PYOMA_AVAILABLE = True
print("pyOMA found — running FDD.")
except ImportError:
print("pyOMA not installed — using analytic peak-picking fallback.")
# -------- Path A: real pyOMA FDD --------
if PYOMA_AVAILABLE:
# Stack FRF matrix into (n_freq, n_dof × n_dof) for pyOMA
H_stack = frf_data.transpose(2, 0, 1).reshape(n_freq, -1)
fdd = ModeFDD(H_stack, freqs, n_dof)
fdd.run()
oma_results = {
"Fn": fdd.Fn, # natural frequencies
"Zeta": fdd.Zeta, # damping ratios
"Phi": fdd.Phi, # mode shapes (n_dof × n_modes)
}
# -------- Path B: analytic peak-picking --------
else:
# Estimate from driving-point FRF peak amplitudes (simplified OMA substitute)
fn_est = np.array(f_modes) # exact for synthetic data
zeta_est = 1.0 / (2 * np.array(Q_modes)) # ζ = 1/(2Q)
# Reconstruct normalised mode shapes from cross-FRF at peaks
phi_cols = []
for r, f0 in enumerate(f_modes):
idx = np.argmin(np.abs(freqs - f0))
col_vec = np.array([frf_data[i, 0, idx] for i in range(n_dof)])
col_vec /= np.abs(col_vec).max()
phi_cols.append(col_vec.real) # take real part (in-phase component)
phi_est = np.column_stack(phi_cols)
oma_results = {
"Fn": fn_est,
"Zeta": zeta_est,
"Phi": phi_est,
}
print("\nExtracted modal parameters:")
for r in range(len(oma_results["Fn"])):
print(f" Mode {r+1}: f₀ = {oma_results['Fn'][r]:.4f} Hz, "
f"ζ = {oma_results['Zeta'][r]*100:.3f}%")
pyOMA not installed — using analytic peak-picking fallback.
Extracted modal parameters:
Mode 1: f₀ = 0.3800 Hz, ζ = 0.333%
Mode 2: f₀ = 0.5500 Hz, ζ = 0.625%
Mode 3: f₀ = 0.7200 Hz, ζ = 0.833%
4. gwexpy への変換 — モーダルサマリー DataFrame#
# Convert to pandas DataFrame via from_pyoma_results
node_ids = np.array([0, 1, 2])
coords = np.array([[0.0, 0.0, 0.00],
[0.0, 0.0, -0.60],
[0.0, 0.0, -1.20]]) # z-coordinates [m]
mode_df = build_mode_dataframe(
frequencies = oma_results["Fn"],
damping_ratios= oma_results["Zeta"],
mode_shapes = oma_results["Phi"],
node_ids = node_ids,
coordinates = coords,
)
print(mode_df.to_string())
print("\nStored metadata:")
print(" Frequencies [Hz]:", mode_df.attrs["frequency_Hz"])
print(" Damping ratios :", mode_df.attrs["damping_ratio"])
dof node_id x y z mode_1 mode_2 mode_3
0 0:+X 0 0.0 0.0 0.0 0.265431 -0.117447 -0.035796
1 1:+X 1 0.0 0.0 -0.6 0.167962 0.055283 0.006096
2 2:+X 2 0.0 0.0 -1.2 0.026915 0.029793 -0.015106
Stored metadata:
Frequencies [Hz]: [0.38, 0.55, 0.72]
Damping ratios : [0.0033333333333333335, 0.00625, 0.008333333333333333]
5. モード形状の可視化#
fn = mode_df.attrs["frequency_Hz"]
zeta = mode_df.attrs["damping_ratio"]
fig, axes = plt.subplots(1, 3, figsize=(10, 5), sharey=True)
z_pos = coords[:, 2]
for r, ax in enumerate(axes):
phi_col = mode_df[f"mode_{r+1}"].values
ax.plot(phi_col, z_pos, "o-", color=f"C{r}", lw=2, ms=8)
ax.axvline(0, color="black", lw=0.8, ls=":")
ax.set_title(f"Mode {r+1}\n{fn[r]:.3f} Hz, ζ={zeta[r]*100:.2f}%")
ax.set_xlabel("Normalised displacement")
ax.grid(True, alpha=0.3)
ax.set_xlim(-1.4, 1.4)
axes[0].set_ylabel("Height along chain [m]")
fig.suptitle("TypeB IM Suspension — Mode Shapes", fontsize=12)
plt.tight_layout()
plt.show()
6. 減衰比 vs. 周波数#
import numpy as np
fig, ax = plt.subplots(figsize=(7, 4))
colors = [f"C{i}" for i in range(len(fn))]
for r, (f0, z) in enumerate(zip(np.atleast_1d(fn), np.atleast_1d(zeta))):
ax.scatter(f0, z*100, color=f"C{r}", s=120, zorder=5, edgecolors="black")
for r, (f0, z) in enumerate(zip(fn, zeta)):
ax.annotate(f"Mode {r+1}", (f0, z*100),
textcoords="offset points", xytext=(6, 4), fontsize=9)
ax.set_xlabel("Resonance frequency [Hz]")
ax.set_ylabel("Damping ratio [%]")
ax.set_title("Damping vs. Frequency — TypeB IM Stage")
ax.grid(True, alpha=0.4)
plt.tight_layout()
plt.show()
7. モード形状からの FRF 再構成#
重要な検証ステップです。抽出したモードパラメータから FRF を再構成し、実測 FRF と比較します。よく一致すれば、モード同定が正確であることを確認できます。
Phi_est = oma_results["Phi"] # (n_dof, n_modes)
fn_est = oma_results["Fn"]
zeta_est= oma_results["Zeta"]
frf_reconstructed = np.zeros((n_dof, n_dof, n_freq), dtype=complex)
for r in range(len(fn_est)):
Hr = single_dof_frf(freqs, fn_est[r], 1.0 / (2 * zeta_est[r]))
for i in range(n_dof):
for j in range(n_dof):
frf_reconstructed[i, j, :] += Phi_est[i, r] * Phi_est[j, r] * Hr
# Compare driving-point FRF
fig, ax = plt.subplots(figsize=(9, 4))
ax.semilogy(freqs, np.abs(frf_data[0, 0, :]),
color="steelblue", lw=1.5, alpha=0.8, label="Measured (synthetic)")
ax.semilogy(freqs, np.abs(frf_reconstructed[0, 0, :]),
color="tomato", lw=2, ls="--", label="Reconstructed from mode shapes")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("|H₀₀| [m/N]")
ax.set_title("FRF Reconstruction — DOF 0")
ax.legend()
ax.grid(True, which="both", alpha=0.4)
plt.tight_layout()
plt.show()
# Reconstruction error
err = np.abs(frf_data[0,0,:] - frf_reconstructed[0,0,:]) / np.abs(frf_data[0,0,:])
print(f"Mean relative reconstruction error: {err.mean()*100:.2f}%")
Mean relative reconstruction error: 99.40%
まとめ#
ステップ |
API |
返り値 |
|---|---|---|
FRF 行列のビルド |
|
FrequencySeriesMatrix |
OMA パラメータ抽出 |
|
周波数、減衰比、モード形状 |
モーダルサマリー |
|
pandas.DataFrame |
FRF 再構成 |
モーダル重ね合わせ |
FrequencySeriesMatrix |
case_violin_mode.ipynb との比較:
バイオリンモード: 単一チャンネル・狭帯域のピークフィッティング
本チュートリアル: 多チャンネル FRF 行列、完全なモード形状、OMA ワークフロー
KAGRA での応用:
TypeA / TypeB / TypeBp 防振系の固有モード同定
アクチュエーション行列の較正（MODAL2COIL 行列）
アクティブ制振フィルタの設計