当前位置：首页 > news >正文

基于POD和DMD的压气机叶片瞬态流场分析与神经网络预测

news 2025/7/27 5:55:25

基于POD和DMD的压气机叶片瞬态流场分析与神经网络预测

前些天发现了一个巨牛的人工智能学习网站，通俗易懂，风趣幽默，忍不住分享一下给大家，觉得好请收藏。点击跳转到网站。

1. 引言

1.1 研究背景与意义

压气机作为航空发动机和燃气轮机的核心部件，其内部流场特性直接影响整机性能。压气机叶片流场具有高度复杂的三维非定常特性，包含多种流动现象如激波、边界层分离、二次流等。准确分析和预测这些瞬态流场特性对于提高压气机效率、扩大稳定工作范围具有重要意义。

传统计算流体力学(CFD)方法虽然能够提供详细的流场信息，但计算成本高昂，难以满足工程实时分析和优化的需求。降阶模型(Reduced Order Model, ROM)技术通过提取流场主导特征，可以大幅降低模型复杂度，同时保持关键物理特性。其中，本征正交分解(Proper Orthogonal Decomposition, POD)和动态模态分解(Dynamic Mode Decomposition, DMD)是两种广泛应用于流体力学领域的降阶方法。

近年来，随着深度学习技术的快速发展，神经网络在流场重构和预测方面展现出强大潜力。将传统降阶方法与神经网络相结合，可以充分发挥各自优势，为压气机流场分析提供新的技术途径。

1.2 研究内容与方法

本研究基于压气机叶片99%叶高处的瞬态流场数据(包括坐标、压力和速度)，采用POD和DMD方法对流场进行分解与重构，并对重构结果进行可视化分析。在此基础上，构建神经网络模型实现流场重构预测和稳定性预测。主要研究内容包括：

流场数据预处理与特征分析
POD方法实现与模态分析
DMD方法实现与频率特性分析
基于神经网络的重构预测模型
流场稳定性预测模型
结果可视化与对比分析

1.3 技术路线

本研究的技术路线如图1所示：

原始流场数据 → 数据预处理 → POD/DMD分解 → 模态分析 → 流场重构 → 重构误差分析↓                ↓神经网络模型训练 ← 特征提取 → 稳定性预测

2. 数据预处理与特征分析

2.1 数据加载与初步处理

我们首先加载压气机叶片99%叶高处的瞬态流场数据。假设数据以CSV格式存储，包含时间序列的坐标(x,y)、压力§和速度(u,v,w)信息。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D# 加载数据
def load_data(file_path):data = pd.read_csv(file_path)# 假设数据列名为: time, x, y, p, u, v, wtimes = np.unique(data['time'].values)coords = data[['x', 'y']].drop_duplicates().valuesn_times = len(times)n_points = len(coords)# 初始化数据数组p_data = np.zeros((n_points, n_times))u_data = np.zeros((n_points, n_times))v_data = np.zeros((n_points, n_times))w_data = np.zeros((n_points, n_times))# 填充数据for i, t in enumerate(times):time_data = data[data['time'] == t]p_data[:, i] = time_data['p'].valuesu_data[:, i] = time_data['u'].valuesv_data[:, i] = time_data['v'].valuesw_data[:, i] = time_data['w'].valuesreturn coords, times, p_data, u_data, v_data, w_data# 示例数据加载
file_path = 'compressor_blade_99percent.csv'
coords, times, p_data, u_data, v_data, w_data = load_data(file_path)

2.2 数据可视化

为了直观了解流场特性，我们对某一时刻的流场进行可视化：

def plot_flow_field(coords, field, title, cmap='jet'):plt.figure(figsize=(10, 6))plt.scatter(coords[:, 0], coords[:, 1], c=field, cmap=cmap, s=5)plt.colorbar(label=title)plt.xlabel('x coordinate')plt.ylabel('y coordinate')plt.title(f'{title} distribution at t={times[0]:.4f}s')plt.axis('equal')plt.show()# 绘制压力场
plot_flow_field(coords, p_data[:, 0], 'Pressure')# 绘制速度场
velocity_magnitude = np.sqrt(u_data[:, 0]**2 + v_data[:, 0]**2 + w_data[:, 0]**2)
plot_flow_field(coords, velocity_magnitude, 'Velocity Magnitude')

2.3 数据标准化

为了提高后续分析的数值稳定性，我们对数据进行标准化处理：

from sklearn.preprocessing import StandardScalerdef standardize_data(data):scaler = StandardScaler()data_std = scaler.fit_transform(data)return data_std, scaler# 标准化压力数据
p_data_std, p_scaler = standardize_data(p_data)# 标准化速度数据
u_data_std, u_scaler = standardize_data(u_data)
v_data_std, v_scaler = standardize_data(v_data)
w_data_std, w_scaler = standardize_data(w_data)

3. POD方法实现与分析

3.1 POD理论基础

本征正交分解(POD)是一种将高维动态系统降维到低维空间的技术，通过寻找最优基函数来捕捉流场中的主要能量结构。POD的核心是求解数据协方差矩阵的特征值问题，得到POD模态及其对应的能量贡献。

数学上，给定流场快照矩阵 $\in \mathbb{R}^{m \times n}$ ，其中 $m$ 为空间点数， $n$ 为时间步数。POD求解以下特征问题：

$XXTϕi=λiϕiXX^T \phi_i = \lambda_i \phi_i$

其中 $ϕi\phi_i$ 为POD模态， $λi\lambda_i$ 为对应的特征值，表征模态的能量贡献。

3.2 POD实现

我们使用奇异值分解(SVD)来实现POD分析：

def pod_analysis(data, n_modes=None):# 数据矩阵: m points × n snapshotsU, s, Vh = np.linalg.svd(data, full_matrices=False)# 计算能量贡献energy = s**2 / np.sum(s**2)cum_energy = np.cumsum(energy)if n_modes is None:# 自动确定模态数，覆盖95%能量n_modes = np.argmax(cum_energy >= 0.95) + 1# 选择前n_modes个模态modes = U[:, :n_modes]coefficients = np.diag(s[:n_modes]) @ Vh[:n_modes, :]return modes, coefficients, energy, cum_energy# 对压力场进行POD分析
p_modes, p_coeffs, p_energy, p_cum_energy = pod_analysis(p_data_std)# 对速度场进行POD分析
u_modes, u_coeffs, u_energy, u_cum_energy = pod_analysis(u_data_std)
v_modes, v_coeffs, v_energy, v_cum_energy = pod_analysis(v_data_std)
w_modes, w_coeffs, w_energy, w_cum_energy = pod_analysis(w_data_std)

3.3 POD模态可视化

可视化前几个POD模态及其能量贡献：

def plot_pod_modes(coords, modes, energy, n_modes=4):fig, axes = plt.subplots(2, 2, figsize=(12, 10))axes = axes.ravel()for i in range(n_modes):ax = axes[i]sc = ax.scatter(coords[:, 0], coords[:, 1], c=modes[:, i], cmap='coolwarm', s=5)fig.colorbar(sc, ax=ax)ax.set_title(f'POD Mode {i+1}, Energy: {energy[i]*100:.1f}%')ax.axis('equal')plt.tight_layout()plt.show()# 绘制压力POD模态
plot_pod_modes(coords, p_modes, p_energy)# 绘制速度POD模态
plot_pod_modes(coords, u_modes, u_energy)def plot_energy_spectrum(energy, cum_energy):plt.figure(figsize=(10, 5))plt.subplot(1, 2, 1)plt.bar(range(1, len(energy)+1), energy[:len(energy)])plt.xlabel('Mode number')plt.ylabel('Energy contribution')plt.title('POD Energy Spectrum')plt.subplot(1, 2, 2)plt.plot(range(1, len(cum_energy)+1), cum_energy[:len(cum_energy)])plt.axhline(y=0.95, color='r', linestyle='--')plt.xlabel('Mode number')plt.ylabel('Cumulative energy')plt.title('Cumulative Energy')plt.tight_layout()plt.show()# 绘制压力能量谱
plot_energy_spectrum(p_energy, p_cum_energy)

3.4 流场重构

利用POD模态和系数重构原始流场：

def reconstruct_field(modes, coefficients, n_modes=None):if n_modes is None:n_modes = modes.shape[1]return modes[:, :n_modes] @ coefficients[:n_modes, :]# 使用不同数量的模态进行重构
n_modes_list = [1, 5, 10, 20]
for n in n_modes_list:p_recon = reconstruct_field(p_modes, p_coeffs, n)# 反标准化p_recon = p_scaler.inverse_transform(p_recon.T).T# 计算重构误差error = np.linalg.norm(p_data - p_recon) / np.linalg.norm(p_data)print(f'Reconstruction with {n} modes, relative error: {error*100:.2f}%')# 可视化重构结果plot_flow_field(coords, p_recon[:, 0], f'Reconstructed Pressure (Modes={n})')

4. DMD方法实现与分析

4.1 DMD理论基础

动态模态分解(DMD)是一种从时间序列数据中提取动态特征的方法，可以识别流场中的主导频率和增长/衰减模式。DMD假设数据可以由线性动力学系统生成，并求解该系统的特征值和特征向量。

DMD的核心是求解以下近似关系：

$\approx AX$

其中 $X$ 和 $X^{'}$ 分别是时间序列的当前和下一时刻数据矩阵， $A$ 为线性动力学算子。DMD模态是 $A$ 的特征向量，对应的特征值表征模态的动态特性。

4.2 DMD实现

我们使用精确DMD算法实现：

def dmd_analysis(data, dt, rank=None):# 数据矩阵: m points × n snapshotsX = data[:, :-1]Y = data[:, 1:]# SVD降维U, s, Vh = np.linalg.svd(X, full_matrices=False)if rank is not None:U = U[:, :rank]s = s[:rank]Vh = Vh[:rank, :]# 构建降维算子A_tilde = U.T @ Y @ Vh.T @ np.diag(1/s)# 求解特征值和特征向量eigvals, eigvecs = np.linalg.eig(A_tilde)# 计算DMD模态modes = Y @ Vh.T @ np.diag(1/s) @ eigvecs# 计算频率和增长率omega = np.log(eigvals) / dtfrequencies = np.imag(omega) / (2*np.pi)growth_rates = np.real(omega)return modes, eigvals, frequencies, growth_rates# 计算时间步长
dt = times[1] - times[0]# 对压力场进行DMD分析
p_modes_dmd, p_eigvals, p_freqs, p_growth = dmd_analysis(p_data_std, dt, rank=20)# 对速度场进行DMD分析
u_modes_dmd, u_eigvals, u_freqs, u_growth = dmd_analysis(u_data_std, dt, rank=20)

4.3 DMD模态可视化

可视化DMD模态及其频率特性：

def plot_dmd_modes(coords, modes, freqs, growth_rates, n_modes=4):fig, axes = plt.subplots(2, 2, figsize=(12, 10))axes = axes.ravel()# 按能量排序(使用模态幅值作为代理)mode_energy = np.abs(modes).sum(axis=0)sorted_idx = np.argsort(mode_energy)[::-1]for i in range(n_modes):idx = sorted_idx[i]ax = axes[i]mode = modes[:, idx]# 可视化模态实部sc = ax.scatter(coords[:, 0], coords[:, 1], c=np.real(mode), cmap='coolwarm', s=5)fig.colorbar(sc, ax=ax)ax.set_title(f'DMD Mode {i+1}\nFreq: {freqs[idx]:.2f}Hz, Growth: {growth_rates[idx]:.2f}')ax.axis('equal')plt.tight_layout()plt.show()# 绘制压力DMD模态
plot_dmd_modes(coords, p_modes_dmd, p_freqs, p_growth)def plot_dmd_spectrum(frequencies, growth_rates):plt.figure(figsize=(8, 6))plt.scatter(frequencies, growth_rates, alpha=0.6)plt.axhline(y=0, color='k', linestyle='--')plt.xlabel('Frequency (Hz)')plt.ylabel('Growth Rate')plt.title('DMD Spectrum')plt.grid(True)plt.show()# 绘制DMD频谱
plot_dmd_spectrum(p_freqs, p_growth)

4.4 基于DMD的流场预测

利用DMD模态进行流场预测：

def dmd_predict(modes, eigvals, initial_condition, n_steps):# 计算初始条件在DMD模态上的投影系数b = np.linalg.pinv(modes) @ initial_condition# 初始化预测矩阵prediction = np.zeros((modes.shape[0], n_steps), dtype=complex)# 逐步预测for i in range(n_steps):prediction[:, i] = modes @ (b * eigvals**i)return prediction# 使用第一个快照作为初始条件
initial_condition = p_data_std[:, 0]# 预测后续50个时间步
n_predict = 50
p_pred_dmd = dmd_predict(p_modes_dmd, p_eigvals, initial_condition, n_predict)# 可视化预测结果
def plot_prediction(coords, true_data, pred_data, time_indices):n_plots = len(time_indices)fig, axes = plt.subplots(2, n_plots, figsize=(5*n_plots, 8))for i, t_idx in enumerate(time_indices):# 真实数据ax = axes[0, i]sc = ax.scatter(coords[:, 0], coords[:, 1], c=true_data[:, t_idx], cmap='jet', s=5)fig.colorbar(sc, ax=ax)ax.set_title(f'True t={times[t_idx]:.3f}s')ax.axis('equal')# 预测数据ax = axes[1, i]sc = ax.scatter(coords[:, 0], coords[:, 1], c=np.real(pred_data[:, t_idx]), cmap='jet', s=5)fig.colorbar(sc, ax=ax)ax.set_title(f'Predicted t={times[t_idx]:.3f}s')ax.axis('equal')plt.tight_layout()plt.show()# 选择几个时间点进行可视化
plot_indices = [0, 10, 20, 30]
plot_prediction(coords, p_data_std, p_pred_dmd, plot_indices)

5. 基于神经网络的流场重构与预测

5.1 神经网络模型设计

我们构建一个自动编码器(Autoencoder)结构进行流场特征提取和重构：

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.model_selection import train_test_split# 定义自动编码器模型
class FlowFieldAutoencoder(nn.Module):def __init__(self, input_dim, latent_dim):super(FlowFieldAutoencoder, self).__init__()# 编码器self.encoder = nn.Sequential(nn.Linear(input_dim, 256),nn.ReLU(),nn.Linear(256, 128),nn.ReLU(),nn.Linear(128, latent_dim))# 解码器self.decoder = nn.Sequential(nn.Linear(latent_dim, 128),nn.ReLU(),nn.Linear(128, 256),nn.ReLU(),nn.Linear(256, input_dim))def forward(self, x):encoded = self.encoder(x)decoded = self.decoder(encoded)return decoded# 定义LSTM预测模型
class FlowPredictor(nn.Module):def __init__(self, input_dim, hidden_dim, output_dim, n_layers=2):super(FlowPredictor, self).__init__()self.lstm = nn.LSTM(input_dim, hidden_dim, n_layers, batch_first=True)self.fc = nn.Linear(hidden_dim, output_dim)def forward(self, x):lstm_out, _ = self.lstm(x)output = self.fc(lstm_out[:, -1, :])return output

5.2 数据准备与模型训练

准备训练数据并训练自动编码器：

# 准备自动编码器数据
X_ae = p_data_std.T  # 时间步 × 空间点
X_train_ae, X_test_ae = train_test_split(X_ae, test_size=0.2, random_state=42)# 转换为PyTorch张量
train_data_ae = torch.FloatTensor(X_train_ae)
test_data_ae = torch.FloatTensor(X_test_ae)
train_loader_ae = DataLoader(TensorDataset(train_data_ae, train_data_ae), batch_size=32, shuffle=True)# 初始化模型
input_dim = X_ae.shape[1]
latent_dim = 10  # 与POD模态数相当
autoencoder = FlowFieldAutoencoder(input_dim, latent_dim)# 训练参数
criterion = nn.MSELoss()
optimizer = optim.Adam(autoencoder.parameters(), lr=0.001)
n_epochs = 100# 训练循环
train_losses = []
test_losses = []for epoch in range(n_epochs):# 训练阶段autoencoder.train()train_loss = 0for batch in train_loader_ae:inputs, targets = batchoptimizer.zero_grad()outputs = autoencoder(inputs)loss = criterion(outputs, targets)loss.backward()optimizer.step()train_loss += loss.item()train_loss /= len(train_loader_ae)train_losses.append(train_loss)# 测试阶段autoencoder.eval()with torch.no_grad():test_outputs = autoencoder(test_data_ae)test_loss = criterion(test_outputs, test_data_ae).item()test_losses.append(test_loss)if (epoch+1) % 10 == 0:print(f'Epoch {epoch+1}/{n_epochs}, Train Loss: {train_loss:.4f}, Test Loss: {test_loss:.4f}')# 绘制训练曲线
plt.figure(figsize=(10, 5))
plt.plot(train_losses, label='Train Loss')
plt.plot(test_losses, label='Test Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Autoencoder Training')
plt.legend()
plt.show()

5.3 流场重构与可视化

使用训练好的自动编码器进行流场重构：

# 使用自动编码器重构流场
autoencoder.eval()
with torch.no_grad():p_recon_nn = autoencoder(torch.FloatTensor(p_data_std.T)).numpy().T# 反标准化
p_recon_nn = p_scaler.inverse_transform(p_recon_nn.T).T# 计算重构误差
error_nn = np.linalg.norm(p_data - p_recon_nn) / np.linalg.norm(p_data)
print(f'Neural Network Reconstruction Error: {error_nn*100:.2f}%')# 可视化重构结果
plot_flow_field(coords, p_recon_nn[:, 0], 'NN Reconstructed Pressure')

5.4 流场预测模型

构建LSTM模型进行流场时间序列预测：

# 准备预测数据
def create_sequences(data, seq_length):sequences = []for i in range(len(data) - seq_length):seq = data[i:i+seq_length]target = data[i+seq_length]sequences.append((seq, target))return sequencesseq_length = 5
sequences = create_sequences(p_data_std, seq_length)
X_seq = np.array([s[0] for s in sequences])
y_seq = np.array([s[1] for s in sequences])# 划分训练测试集
X_train_seq, X_test_seq, y_train_seq, y_test_seq = train_test_split(X_seq, y_seq, test_size=0.2, random_state=42)# 转换为PyTorch张量
train_data_seq = torch.FloatTensor(X_train_seq)
train_target_seq = torch.FloatTensor(y_train_seq)
test_data_seq = torch.FloatTensor(X_test_seq)
test_target_seq = torch.FloatTensor(y_test_seq)train_loader_seq = DataLoader(TensorDataset(train_data_seq, train_target_seq), batch_size=32, shuffle=True)# 初始化LSTM模型
input_dim = p_data_std.shape[0]
hidden_dim = 64
output_dim = input_dim
lstm_predictor = FlowPredictor(input_dim, hidden_dim, output_dim)# 训练参数
criterion = nn.MSELoss()
optimizer = optim.Adam(lstm_predictor.parameters(), lr=0.001)
n_epochs = 50# 训练循环
train_losses_seq = []
test_losses_seq = []for epoch in range(n_epochs):# 训练阶段lstm_predictor.train()train_loss = 0for batch in train_loader_seq:inputs, targets = batchoptimizer.zero_grad()outputs = lstm_predictor(inputs)loss = criterion(outputs, targets)loss.backward()optimizer.step()train_loss += loss.item()train_loss /= len(train_loader_seq)train_losses_seq.append(train_loss)# 测试阶段lstm_predictor.eval()with torch.no_grad():test_outputs = lstm_predictor(test_data_seq)test_loss = criterion(test_outputs, test_target_seq).item()test_losses_seq.append(test_loss)if (epoch+1) % 10 == 0:print(f'Epoch {epoch+1}/{n_epochs}, Train Loss: {train_loss:.4f}, Test Loss: {test_loss:.4f}')# 绘制训练曲线
plt.figure(figsize=(10, 5))
plt.plot(train_losses_seq, label='Train Loss')
plt.plot(test_losses_seq, label='Test Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('LSTM Predictor Training')
plt.legend()
plt.show()

5.5 流场预测与可视化

使用训练好的LSTM模型进行多步预测：

def multi_step_prediction(model, initial_seq, n_steps):predictions = []current_seq = initial_seq.clone()with torch.no_grad():for _ in range(n_steps):# 预测下一步next_step = model(current_seq.unsqueeze(0))predictions.append(next_step.squeeze(0).numpy())# 更新序列current_seq = torch.cat([current_seq[1:], next_step], dim=0)return np.array(predictions).T# 使用测试集第一个序列作为初始条件
initial_sequence = test_data_seq[0]# 预测20步
n_predict_steps = 20
p_pred_lstm = multi_step_prediction(lstm_predictor, initial_sequence, n_predict_steps)# 反标准化
p_pred_lstm = p_scaler.inverse_transform(p_pred_lstm.T).T# 可视化预测结果
plot_indices = [0, 5, 10, 15]
true_data = p_data[:, seq_length:seq_length+n_predict_steps]
plot_prediction(coords, true_data, p_pred_lstm, plot_indices)

6. 稳定性分析与预测

6.1 基于POD系数的稳定性分析

POD时间系数可以反映流场的动态特性，可用于稳定性分析：

def analyze_stability_pod(coefficients, window_size=10):n_modes = coefficients.shape[0]stability_metrics = []for i in range(n_modes):# 计算系数的移动标准差mode_coeff = coefficients[i, :]std_dev = np.array([np.std(mode_coeff[j:j+window_size]) for j in range(len(mode_coeff)-window_size)])# 稳定性指标: 标准差的变化率stability = np.gradient(std_dev)stability_metrics.append(stability)return stability_metrics# 分析压力场POD系数稳定性
p_stability = analyze_stability_pod(p_coeffs)# 可视化稳定性指标
plt.figure(figsize=(12, 6))
for i, stability in enumerate(p_stability[:4]):plt.plot(times[:-10], stability, label=f'Mode {i+1}')
plt.xlabel('Time')
plt.ylabel('Stability Metric')
plt.title('POD Mode Stability Analysis')
plt.legend()
plt.grid(True)
plt.show()

6.2 基于神经网络的稳定性预测

构建稳定性预测模型：

class StabilityPredictor(nn.Module):def __init__(self, input_dim, hidden_dim, n_layers=2):super(StabilityPredictor, self).__init__()self.lstm = nn.LSTM(input_dim, hidden_dim, n_layers, batch_first=True)self.fc = nn.Linear(hidden_dim, 1)  # 预测稳定性指标def forward(self, x):lstm_out, _ = self.lstm(x)output = self.fc(lstm_out[:, -1, :])return output# 准备稳定性标签 (使用POD第一模态稳定性指标作为示例)
stability_labels = p_stability[0][1:]  # 去掉第一个时间点
stability_features = p_data_std[:, :-1].T  # 使用前一时间步的流场预测下一时间步的稳定性# 创建序列数据
seq_length = 5
stability_sequences = create_sequences(np.column_stack((stability_features, stability_labels)), seq_length)
X_stab = np.array([s[0][:, :-1] for s in stability_sequences])  # 特征
y_stab = np.array([s[1][-1] for s in stability_sequences])     # 标签# 划分训练测试集
X_train_stab, X_test_stab, y_train_stab, y_test_stab = train_test_split(X_stab, y_stab, test_size=0.2, random_state=42)# 转换为PyTorch张量
train_data_stab = torch.FloatTensor(X_train_stab)
train_target_stab = torch.FloatTensor(y_train_stab).unsqueeze(1)
test_data_stab = torch.FloatTensor(X_test_stab)
test_target_stab = torch.FloatTensor(y_test_stab).unsqueeze(1)train_loader_stab = DataLoader(TensorDataset(train_data_stab, train_target_stab), batch_size=32, shuffle=True)# 初始化模型
input_dim = p_data_std.shape[0]
hidden_dim = 32
stability_predictor = StabilityPredictor(input_dim, hidden_dim)# 训练参数
criterion = nn.MSELoss()
optimizer = optim.Adam(stability_predictor.parameters(), lr=0.001)
n_epochs = 30# 训练循环
train_losses_stab = []
test_losses_stab = []for epoch in range(n_epochs):# 训练阶段stability_predictor.train()train_loss = 0for batch in train_loader_stab:inputs, targets = batchoptimizer.zero_grad()outputs = stability_predictor(inputs)loss = criterion(outputs, targets)loss.backward()optimizer.step()train_loss += loss.item()train_loss /= len(train_loader_stab)train_losses_stab.append(train_loss)# 测试阶段stability_predictor.eval()with torch.no_grad():test_outputs = stability_predictor(test_data_stab)test_loss = criterion(test_outputs, test_target_stab).item()test_losses_stab.append(test_loss)if (epoch+1) % 5 == 0:print(f'Epoch {epoch+1}/{n_epochs}, Train Loss: {train_loss:.4f}, Test Loss: {test_loss:.4f}')# 绘制训练曲线
plt.figure(figsize=(10, 5))
plt.plot(train_losses_stab, label='Train Loss')
plt.plot(test_losses_stab, label='Test Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Stability Predictor Training')
plt.legend()
plt.show()# 可视化预测结果
stability_predictor.eval()
with torch.no_grad():y_pred_stab = stability_predictor(test_data_stab).numpy().flatten()plt.figure(figsize=(10, 5))
plt.scatter(y_test_stab, y_pred_stab, alpha=0.6)
plt.plot([y_test_stab.min(), y_test_stab.max()], [y_test_stab.min(), y_test_stab.max()], 'r--')
plt.xlabel('True Stability Metric')
plt.ylabel('Predicted Stability Metric')
plt.title('Stability Prediction Results')
plt.grid(True)
plt.show()

7. 结果分析与讨论

7.1 方法对比分析

我们对POD、DMD和神经网络三种方法的重构和预测性能进行对比：

# 计算各方法的重构误差
p_recon_pod = reconstruct_field(p_modes, p_coeffs, n_modes=10)
p_recon_pod = p_scaler.inverse_transform(p_recon_pod.T).T
pod_error = np.linalg.norm(p_data - p_recon_pod) / np.linalg.norm(p_data)p_recon_dmd = np.real(p_pred_dmd[:, :len(times)])
p_recon_dmd = p_scaler.inverse_transform(p_recon_dmd.T).T
dmd_error = np.linalg.norm(p_data - p_recon_dmd[:, :p_data.shape[1]]) / np.linalg.norm(p_data)nn_error = np.linalg.norm(p_data - p_recon_nn) / np.linalg.norm(p_data)# 绘制误差对比
methods = ['POD (10 modes)', 'DMD', 'Neural Network']
errors = [pod_error, dmd_error, nn_error]plt.figure(figsize=(8, 5))
plt.bar(methods, errors)
plt.ylabel('Relative Reconstruction Error')
plt.title('Method Comparison')
plt.xticks(rotation=15)
plt.grid(True, axis='y')
plt.show()

7.2 计算效率分析

比较各方法的计算时间：

import time# POD计算时间
start = time.time()
_, _, _, _ = pod_analysis(p_data_std, n_modes=10)
pod_time = time.time() - start# DMD计算时间
start = time.time()
_, _, _, _ = dmd_analysis(p_data_std, dt, rank=20)
dmd_time = time.time() - start# 神经网络推断时间
autoencoder.eval()
sample_input = torch.FloatTensor(p_data_std[:, 0:1].T)
start = time.time()
with torch.no_grad():_ = autoencoder(sample_input)
nn_time = (time.time() - start) * p_data.shape[1]  # 估计全部时间步时间times = [pod_time, dmd_time, nn_time]plt.figure(figsize=(8, 5))
plt.bar(methods, times)
plt.ylabel('Computation Time (s)')
plt.title('Computational Efficiency Comparison')
plt.xticks(rotation=15)
plt.grid(True, axis='y')
plt.show()