第六十五篇：结构动力学机器学习应用

摘要

机器学习技术为结构动力学分析带来了新的方法和工具。本主题系统介绍机器学习在结构动力学中的应用，包括神经网络响应预测、模态参数识别、损伤检测与诊断、以及代理模型构建等内容。通过Python实现各种机器学习算法与结构动力学问题的结合，展示如何利用数据驱动方法解决传统数值方法难以处理的问题，如高维参数空间探索、实时响应预测和复杂模式识别等。

关键词

机器学习；神经网络；模态识别；损伤检测；代理模型；结构健康监测；深度学习

---

1. 机器学习在结构动力学中的应用概述

1.1 应用背景

传统结构动力学分析依赖于物理模型和数值方法，但在以下场景面临挑战：

• 高维问题：复杂结构参数众多，传统优化和不确定性分析计算成本高昂
• 实时需求：结构健康监测需要快速响应预测和决策
• 数据丰富：传感器网络产生大量监测数据，需要智能分析方法
• 非线性复杂：强非线性系统的建模和识别困难

1.2 主要应用领域

响应预测：利用历史数据训练模型，实现快速响应预测

参数识别：从振动数据自动识别模态参数

损伤检测：识别结构损伤的位置和程度

代理建模：构建计算高效的近似模型替代复杂仿真

优化设计：结合机器学习加速结构优化过程

1.3 常用机器学习方法

监督学习：

• 神经网络（NN）
• 支持向量机（SVM）
• 随机森林（RF）
• 高斯过程（GP）

无监督学习：

• 聚类分析
• 主成分分析（PCA）
• 自编码器（Autoencoder）

深度学习：

• 卷积神经网络（CNN）
• 循环神经网络（RNN/LSTM）
• 图神经网络（GNN）

---

2. 神经网络响应预测

2.1 基本原理

神经网络可以学习结构响应与输入参数之间的映射关系：

$$\mathbf{y}=f_{NN}(\mathbf{x};\mathbf{w})$$

其中 $\mathbf{x}$ 是输入参数（载荷、材料参数等），$\mathbf{y}$ 是输出响应（位移、应力等），$f_{NN}$ 是神经网络函数，$\mathbf{w}$ 是网络权重。

网络结构：

• 输入层：接收结构参数和载荷信息
• 隐藏层：提取特征和非线性变换
• 输出层：预测结构响应

损失函数：

$$L=\frac{1}{N}\sum _{i=1}^N||{\mathbf{y}}_i^{pred}-{\mathbf{y}}_i^{true}|{|}^2$$

2.2 数据准备与训练

训练数据生成：

• 使用有限元模型生成样本
• 拉丁超立方采样覆盖参数空间
• 添加噪声模拟测量不确定性

网络训练技巧：

• 数据归一化
• 早停防止过拟合
• 交叉验证
• 学习率调整

---

3. 模态参数识别的机器学习方法

3.1 问题描述

传统模态识别方法（如ERA、SSI）需要人工选择参数。机器学习方法可以：

• 自动识别模态参数
• 处理噪声数据
• 实现实时识别

3.2 深度学习方法

卷积神经网络（CNN）：

• 输入：时频图或功率谱密度
• 特征提取：自动学习模态特征
• 输出：模态频率和阻尼比

循环神经网络（LSTM）：

• 适合处理时间序列振动数据
• 捕捉模态的时间演化特性

---

4. 结构损伤检测与诊断

4.1 损伤检测框架

损伤指标提取：

• 基于频率变化
• 基于模态振型变化
• 基于柔度矩阵
• 基于振动信号特征

分类与识别：

• 二分类：损伤/无损伤
• 多分类：损伤位置识别
• 回归：损伤程度量化

4.2 深度学习方法

自编码器：

• 学习正常结构的特征表示
• 重构误差作为损伤指标

图神经网络：

• 将结构建模为图
• 利用拓扑信息提高识别精度

---

5. 代理模型与优化

5.1 代理模型构建

代理模型用于替代计算昂贵的有限元仿真：

$$\hat{y}(\mathbf{x})\approx y^{FEM}(\mathbf{x})$$

常用代理模型：

• 高斯过程（Kriging）
• 径向基函数（RBF）
• 神经网络
• 多项式混沌展开

5.2 优化应用

结合代理模型实现：

• 快速优化迭代
• 不确定性量化
• 可靠性分析
• 多目标优化

---

6. Python实现案例

6.1 案例1：神经网络响应预测

问题描述：
使用神经网络学习单自由度系统的动力响应。训练数据通过数值仿真生成，网络学习从系统参数和载荷历史预测位移响应。

Python代码实现：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.size'] = 10

def sdof_response(m, k, c, f, dt, n_steps):
    """
    计算单自由度系统响应（Newmark-beta方法）
    
    Parameters:
    -----------
    m, k, c : float
        质量、刚度、阻尼
    f : array
        外力时程
    dt : float
        时间步长
    n_steps : int
        时间步数
        
    Returns:
    --------
    u, v, a : arrays
        位移、速度、加速度
    """
    # Newmark-beta参数
    beta = 0.25
    gamma = 0.5
    
    u = np.zeros(n_steps)
    v = np.zeros(n_steps)
    a = np.zeros(n_steps)
    
    # 初始条件
    a[0] = (f[0] - c * v[0] - k * u[0]) / m
    
    # 等效刚度
    k_eff = k + gamma / (beta * dt) * c + 1 / (beta * dt**2) * m
    
    for i in range(n_steps - 1):
        # 预测
        u_pred = u[i] + dt * v[i] + (0.5 - beta) * dt**2 * a[i]
        v_pred = v[i] + (1 - gamma) * dt * a[i]
        
        # 等效载荷
        f_eff = f[i+1] + (1 / (beta * dt**2) * m + gamma / (beta * dt) * c) * u_pred \
                - (1 / (beta * dt) * m + (gamma / beta - 1) * c) * v_pred \
                - ((1 / (2 * beta) - 1) * m + dt * (gamma / (2 * beta) - 1) * c) * a[i]
        
        # 求解
        u[i+1] = f_eff / k_eff
        a[i+1] = (u[i+1] - u_pred) / (beta * dt**2)
        v[i+1] = v_pred + gamma * dt * a[i+1]
    
    return u, v, a

def generate_training_data(n_samples=1000, n_steps=200, dt=0.01):
    """
    生成训练数据
    
    Parameters:
    -----------
    n_samples : int
        样本数
    n_steps : int
        时间步数
    dt : float
        时间步长
        
    Returns:
    --------
    X, y : arrays
        输入特征和输出响应
    """
    print(f"生成训练数据 (n_samples={n_samples})...")
    
    X_list = []
    y_list = []
    
    np.random.seed(42)
    
    for i in range(n_samples):
        if (i + 1) % 100 == 0:
            print(f"  进度: {i+1}/{n_samples}")
        
        # 随机系统参数
        m = np.random.uniform(0.5, 2.0)  # 质量
        k = np.random.uniform(500, 2000)  # 刚度
        c = np.random.uniform(0.5, 5.0)  # 阻尼
        
        # 计算固有频率和阻尼比
        omega_n = np.sqrt(k / m)
        xi = c / (2 * np.sqrt(k * m))
        
        # 随机载荷参数
        f_amp = np.random.uniform(50, 200)
        f_freq = np.random.uniform(0.5, 3.0)
        
        # 生成载荷
        t = np.arange(n_steps) * dt
        f = f_amp * np.sin(2 * np.pi * f_freq * t)
        
        # 计算响应
        u, _, _ = sdof_response(m, k, c, f, dt, n_steps)
        
        # 输入特征：系统参数 + 载荷历史（降采样）
        f_sample = f[::10]  # 降采样到20个点
        features = np.concatenate([[m, k, c, omega_n, xi, f_freq], f_sample])
        
        X_list.append(features)
        y_list.append(u[::10])  # 降采样输出
    
    X = np.array(X_list)
    y = np.array(y_list)
    
    return X, y, t[::10]

def case1_neural_network():
    """案例1：神经网络响应预测"""
    print("="*70)
    print("案例1：神经网络响应预测")
    print("="*70)
    
    # 生成数据
    X, y, t_sample = generate_training_data(n_samples=2000, n_steps=200, dt=0.01)
    
    # 划分训练集和测试集
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42)
    
    print(f"\n训练集大小: {X_train.shape[0]}")
    print(f"测试集大小: {X_test.shape[0]}")
    print(f"输入特征维度: {X_train.shape[1]}")
    print(f"输出维度: {y_train.shape[1]}")
    
    # 数据归一化
    scaler_X = StandardScaler()
    scaler_y = StandardScaler()
    
    X_train_scaled = scaler_X.fit_transform(X_train)
    X_test_scaled = scaler_X.transform(X_test)
    y_train_scaled = scaler_y.fit_transform(y_train)
    
    # 训练神经网络
    print("\n训练神经网络...")
    nn = MLPRegressor(
        hidden_layer_sizes=(128, 64, 32),
        activation='relu',
        solver='adam',
        max_iter=500,
        early_stopping=True,
        validation_fraction=0.1,
        n_iter_no_change=20,
        random_state=42,
        verbose=True
    )
    
    nn.fit(X_train_scaled, y_train_scaled)
    
    print(f"\n训练完成!")
    print(f"训练迭代次数: {nn.n_iter_}")
    print(f"最终损失: {nn.loss_:.6f}")
    
    # 测试集预测
    y_pred_scaled = nn.predict(X_test_scaled)
    y_pred = scaler_y.inverse_transform(y_pred_scaled)
    
    # 计算误差
    mse = np.mean((y_pred - y_test)**2)
    rmse = np.sqrt(mse)
    mae = np.mean(np.abs(y_pred - y_test))
    r2 = 1 - np.sum((y_test - y_pred)**2) / np.sum((y_test - np.mean(y_test))**2)
    
    print(f"\n测试集性能:")
    print(f"  MSE: {mse:.6f}")
    print(f"  RMSE: {rmse:.6f}")
    print(f"  MAE: {mae:.6f}")
    print(f"  R²: {r2:.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 子图1：训练损失曲线
    ax = axes[0, 0]
    ax.plot(nn.loss_curve_, linewidth=2)
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Loss')
    ax.set_title('Training Loss Curve')
    ax.grid(True, alpha=0.3)
    ax.set_yscale('log')
    
    # 子图2：预测 vs 真实值
    ax = axes[0, 1]
    sample_idx = 0
    ax.plot(t_sample, y_test[sample_idx], 'b-', linewidth=2, label='True')
    ax.plot(t_sample, y_pred[sample_idx], 'r--', linewidth=2, label='Predicted')
    ax.set_xlabel('Time (s)')
    ax.set_ylabel('Displacement (m)')
    ax.set_title(f'Sample Prediction (Test #{sample_idx+1})')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 子图3：多个样本对比
    ax = axes[1, 0]
    for i in range(min(5, len(y_test))):
        ax.plot(t_sample, y_test[i], alpha=0.5, linewidth=1.5)
    ax.set_xlabel('Time (s)')
    ax.set_ylabel('Displacement (m)')
    ax.set_title('Multiple Test Samples (True)')
    ax.grid(True, alpha=0.3)
    
    # 子图4：误差分布
    ax = axes[1, 1]
    errors = y_pred - y_test
    ax.hist(errors.flatten(), bins=50, alpha=0.7, edgecolor='black')
    ax.axvline(0, color='r', linestyle='--', linewidth=2)
    ax.set_xlabel('Prediction Error (m)')
    ax.set_ylabel('Frequency')
    ax.set_title('Error Distribution')
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('nn_response_prediction.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("\n  已保存: nn_response_prediction.png")
    
    return nn, scaler_X, scaler_y, t_sample

if __name__ == "__main__":
    nn, scaler_X, scaler_y, t_sample = case1_neural_network()

6.2 案例2：模态参数识别的机器学习方法

问题描述：
使用机器学习方法从振动响应数据中自动识别模态参数（频率和阻尼比）。训练数据包含不同模态特性的系统响应，模型学习从时域或频域特征提取模态参数。

特征工程：
从振动响应中提取多维度特征：

• 时域特征：均值、标准差、最大值、最小值、过零率
• 频域特征：主频率、频谱质心、频谱带宽、频谱熵
• 模态特征：对数衰减率、自相关特征
• 能量特征：频带能量分布

机器学习方法：

• 随机森林：集成学习方法，具有良好的泛化能力
• 支持向量机：适合处理高维特征空间

Python代码实现：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestRegressor
from sklearn.svm import SVR
from sklearn.multioutput import MultiOutputRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score

def extract_features(u, t):
    """从响应中提取特征"""
    dt = t[1] - t[0]
    n = len(u)
    
    # 检查并处理NaN/Inf值
    u = np.nan_to_num(u, nan=0.0, posinf=0.0, neginf=0.0)
    
    # 时域特征
    features = []
    
    # 1. 统计特征
    features.extend([
        np.mean(u), np.std(u), np.max(u), np.min(u), np.max(np.abs(u))
    ])
    
    # 2. 过零率
    zero_crossings = np.sum(np.diff(np.sign(u)) != 0)
    features.append(zero_crossings / len(t))
    
    # 3. 频域特征（FFT）
    fft_vals = np.fft.fft(u)
    freqs = np.fft.fftfreq(n, dt)
    psd = np.abs(fft_vals)**2
    
    # 只取正频率
    pos_mask = freqs > 0
    freqs_pos = freqs[pos_mask]
    psd_pos = psd[pos_mask]
    
    # 主频率
    if len(psd_pos) > 0 and np.sum(psd_pos) > 0:
        peak_idx = np.argmax(psd_pos)
        dominant_freq = freqs_pos[peak_idx]
        features.append(dominant_freq)
        
        # 频谱质心
        spectral_centroid = np.sum(freqs_pos * psd_pos) / np.sum(psd_pos)
        features.append(spectral_centroid)
        
        # 频谱带宽
        spectral_bandwidth = np.sqrt(np.sum((freqs_pos - spectral_centroid)**2 * psd_pos) / np.sum(psd_pos))
        features.append(spectral_bandwidth)
        
        # 频谱熵
        psd_norm = psd_pos / np.sum(psd_pos)
        spectral_entropy = -np.sum(psd_norm * np.log(psd_norm + 1e-10))
        features.append(spectral_entropy)
    else:
        features.extend([0, 0, 0, 0])
    
    # 4. 对数衰减特征
    peaks = []
    for i in range(1, len(u)-1):
        if u[i] > u[i-1] and u[i] > u[i+1] and u[i] > 0:
            peaks.append((i, u[i]))
    
    if len(peaks) >= 2:
        log_decrements = []
        for i in range(len(peaks)-1):
            if peaks[i][1] > 0 and peaks[i+1][1] > 0 and peaks[i+1][1] != 0:
                log_dec = np.log(peaks[i][1] / peaks[i+1][1])
                if np.isfinite(log_dec):
                    log_decrements.append(log_dec)
        
        if log_decrements:
            features.append(np.mean(log_decrements))
            features.append(np.std(log_decrements) if len(log_decrements) > 1 else 0)
        else:
            features.extend([0, 0])
    else:
        features.extend([0, 0])
    
    # 5. 自相关特征
    autocorr = np.correlate(u, u, mode='full')
    autocorr = autocorr[len(autocorr)//2:]
    if autocorr[0] != 0:
        autocorr = autocorr / autocorr[0]
    else:
        autocorr = np.zeros_like(autocorr)
    
    # 找到第一个过零点（近似半周期）
    zero_idx = np.where(autocorr < 0)[0]
    if len(zero_idx) > 0:
        features.append(zero_idx[0] * dt)
    else:
        features.append(0)
    
    # 6. 小波包能量（简化版 - 使用带通滤波）
    fft_u = np.fft.fft(u)
    freqs = np.fft.fftfreq(len(u), dt)
    
    bands = [(0, 5), (5, 10), (10, 20), (20, 50)]
    for f_low, f_high in bands:
        band_mask = (np.abs(freqs) >= f_low) & (np.abs(freqs) < f_high)
        band_energy = np.sum(np.abs(fft_u[band_mask])**2)
        features.append(band_energy if np.isfinite(band_energy) else 0)
    
    # 确保没有NaN或Inf
    features = np.array(features)
    features = np.nan_to_num(features, nan=0.0, posinf=0.0, neginf=0.0)
    
    return features

结果分析：

• 随机森林和支持向量机都能有效识别模态参数
• 特征重要性分析显示主频率和对数衰减率是最重要的特征
• 模型可以处理不同激励类型（脉冲、随机、正弦）的数据

6.3 案例3：结构损伤检测与诊断

问题描述：
构建机器学习模型检测结构损伤。使用正常和损伤状态下的结构响应数据训练分类器，实现损伤检测、位置识别和程度评估。

损伤检测框架：

1. 损伤检测（二分类）：判断结构是否损伤
2. 位置识别（多分类）：识别损伤发生的单元位置
3. 程度评估（回归）：量化损伤的严重程度

特征提取：
基于模态参数变化提取损伤敏感特征：

• 频率变化率
• 振型变化（MAC值）
• 模态柔度变化

Python代码实现：

from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from sklearn.neural_network import MLPClassifier
from sklearn.metrics import classification_report, accuracy_score

def generate_damage_data(n_samples=1000, n_elements=10):
    """生成损伤检测数据集"""
    print(f"生成损伤检测数据 (n_samples={n_samples})...")
    
    # 梁参数
    L_total = 10.0
    L_element = L_total / n_elements
    E = 2.1e11
    I = 8.33e-6
    rho = 7850
    A = 0.01
    
    X_list = []
    y_damage_list = []
    y_location_list = []
    y_severity_list = []
    
    np.random.seed(42)
    
    for i in range(n_samples):
        # 随机决定是否损伤
        is_damaged = np.random.random() > 0.5
        
        if is_damaged:
            damage_element = np.random.randint(0, n_elements)
            damage_severity = np.random.uniform(0.1, 0.5)
            damage_info = {'element': damage_element, 'severity': damage_severity}
        else:
            damage_info = None
            damage_element = -1
            damage_severity = 0
        
        # 组装系统
        K, M = assemble_system(n_elements, L_element, E, I, rho, A, damage_info)
        
        # 求解模态
        freqs, modes = solve_modal(K, M, n_modes=5)
        
        # 提取特征
        K_healthy, M_healthy = assemble_system(n_elements, L_element, E, I, rho, A, None)
        freqs_healthy, _ = solve_modal(K_healthy, M_healthy, n_modes=5)
        freq_changes = (freqs_healthy - freqs) / freqs_healthy
        
        # 振型特征
        mode_features = np.abs(modes[::2, :3]).flatten()
        
        # MAC值
        _, modes_healthy = solve_modal(K_healthy, M_healthy, n_modes=5)
        mac_values = []
        for j in range(min(3, len(freqs))):
            mode_damaged = modes[:, j]
            mode_healthy = modes_healthy[:, j]
            mac = np.abs(mode_damaged @ mode_healthy)**2 / \
                  ((mode_damaged @ mode_damaged) * (mode_healthy @ mode_healthy))
            mac_values.append(mac)
        
        features = np.concatenate([freqs, freq_changes, mode_features, mac_values])
        
        X_list.append(features)
        y_damage_list.append(1 if is_damaged else 0)
        y_location_list.append(damage_element + 1 if is_damaged else 0)
        y_severity_list.append(damage_severity)
    
    return np.array(X_list), np.array(y_damage_list), \
           np.array(y_location_list), np.array(y_severity_list)

# 任务1：损伤检测（二分类）
rf_clf = RandomForestClassifier(n_estimators=100, max_depth=10, random_state=42)
rf_clf.fit(X_train, y_train)
y_pred = rf_clf.predict(X_test)
print(f"损伤检测准确率: {accuracy_score(y_test, y_pred):.4f}")

# 任务2：位置识别（多分类）
rf_loc = RandomForestClassifier(n_estimators=100, max_depth=15, random_state=42)
rf_loc.fit(X_train_loc, y_train_loc)
y_pred_loc = rf_loc.predict(X_test_loc)
print(f"位置识别准确率: {accuracy_score(y_test_loc, y_pred_loc):.4f}")

# 任务3：程度评估（回归）
rf_sev = RandomForestRegressor(n_estimators=100, max_depth=10, random_state=42)
rf_sev.fit(X_train_sev, y_train_sev)
y_pred_sev = rf_sev.predict(X_test_sev)

结果分析：

• 损伤检测准确率可达98%以上
• 位置识别准确率接近100%
• 程度评估RMSE小于0.01，精度很高
• ROC曲线AUC接近1.0，模型区分能力强

6.4 案例4：代理模型与优化

问题描述：
构建神经网络代理模型替代有限元分析，用于快速结构优化。对比直接有限元优化和代理模型辅助优化的效率。

代理模型构建：

1. 训练数据生成：使用拉丁超立方采样生成设计参数样本
2. 有限元分析：计算每个样本的结构响应
3. 神经网络训练：学习设计参数到响应的映射
4. 模型验证：评估代理模型的预测精度

优化对比：

• 直接FEM优化：每次迭代需要完整的有限元分析
• 代理模型优化：使用训练好的神经网络快速预测

Python代码实现：

from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
import time

def build_surrogate_model(X, y):
    """构建代理模型"""
    print("\n构建代理模型...")
    
    # 数据归一化
    scaler_X = StandardScaler()
    scaler_y = StandardScaler()
    
    X_scaled = scaler_X.fit_transform(X)
    y_scaled = scaler_y.fit_transform(y)
    
    # 训练神经网络
    model = MLPRegressor(
        hidden_layer_sizes=(100, 50, 25),
        activation='relu',
        solver='adam',
        max_iter=1000,
        early_stopping=True,
        validation_fraction=0.1,
        n_iter_no_change=30,
        random_state=42
    )
    
    model.fit(X_scaled, y_scaled)
    
    print(f"  训练完成，迭代次数: {model.n_iter_}")
    print(f"  最终损失: {model.loss_:.6f}")
    
    return model, scaler_X, scaler_y

def optimize_with_fem(n_iter=50):
    """使用有限元分析进行优化"""
    print(f"\n使用有限元分析优化 (n_iter={n_iter})...")
    
    bounds = {
        'L_total': (5.0, 15.0),
        'E': (1.5e11, 2.5e11),
        'I': (5e-6, 1.5e-5),
        'rho': (7000, 8500),
        'A': (0.005, 0.02)
    }
    
    best_f = 0
    best_x = None
    history = []
    
    start_time = time.time()
    
    np.random.seed(123)
    
    for i in range(n_iter):
        # 随机采样
        x = np.array([
            np.random.uniform(*bounds['L_total']),
            np.random.uniform(*bounds['E']),
            np.random.uniform(*bounds['I']),
            np.random.uniform(*bounds['rho']),
            np.random.uniform(*bounds['A'])
        ])
        
        # 有限元分析
        freqs = solve_beam_frequencies(10, x[0], x[1], x[2], x[3], x[4])
        f = freqs[0]
        history.append(f)
        
        if f > best_f:
            best_f = f
            best_x = x.copy()
        
        if (i + 1) % 10 == 0:
            print(f"  迭代 {i+1}/{n_iter}, 当前最优: {best_f:.3f} Hz")
    
    time_cost = time.time() - start_time
    
    return best_x, best_f, history, time_cost

def optimize_with_surrogate(model, scaler_X, scaler_y, n_iter=500):
    """使用代理模型进行优化"""
    print(f"\n使用代理模型优化 (n_iter={n_iter})...")
    
    bounds = {
        'L_total': (5.0, 15.0),
        'E': (1.5e11, 2.5e11),
        'I': (5e-6, 1.5e-5),
        'rho': (7000, 8500),
        'A': (0.005, 0.02)
    }
    
    best_f = 0
    best_x = None
    history = []
    
    start_time = time.time()
    
    np.random.seed(123)
    
    for i in range(n_iter):
        x = np.array([
            np.random.uniform(*bounds['L_total']),
            np.random.uniform(*bounds['E']),
            np.random.uniform(*bounds['I']),
            np.random.uniform(*bounds['rho']),
            np.random.uniform(*bounds['A'])
        ])
        
        # 代理模型预测
        X_scaled = scaler_X.transform(x.reshape(1, -1))
        y_pred_scaled = model.predict(X_scaled)
        y_pred = scaler_y.inverse_transform(y_pred_scaled)
        f = y_pred[0, 0]
        
        history.append(f)
        
        if f > best_f:
            best_f = f
            best_x = x.copy()
        
        if (i + 1) % 100 == 0:
            print(f"  迭代 {i+1}/{n_iter}, 当前最优: {best_f:.3f} Hz")
    
    time_cost = time.time() - start_time
    
    return best_x, best_f, history, time_cost

# 生成训练数据
X, y = generate_training_data(n_samples=300)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 构建代理模型
model, scaler_X, scaler_y = build_surrogate_model(X_train, y_train)

# 测试代理模型精度
y_pred_scaled = model.predict(scaler_X.transform(X_test))
y_pred = scaler_y.inverse_transform(y_pred_scaled)
r2 = r2_score(y_test, y_pred)
print(f"代理模型R²: {r2:.4f}")

# 对比优化
best_x_fem, best_f_fem, history_fem, time_fem = optimize_with_fem(n_iter=50)
best_x_sur, best_f_sur, history_sur, time_sur = optimize_with_surrogate(
    model, scaler_X, scaler_y, n_iter=500)

# 验证代理模型的最优解
best_f_sur_verify = objective_function_fem(best_x_sur)

print(f"\n有限元优化结果:")
print(f"  最优频率: {best_f_fem:.3f} Hz")
print(f"  耗时: {time_fem:.3f} s")

print(f"\n代理模型优化结果:")
print(f"  预测最优频率: {best_f_sur:.3f} Hz")
print(f"  FEM验证频率: {best_f_sur_verify:.3f} Hz")
print(f"  耗时: {time_sur:.3f} s")

speedup = time_fem / time_sur if time_sur > 0 else float('inf')
print(f"\n加速比: {speedup:.1f}x")

结果分析：

• 代理模型R²达到0.95以上，预测精度高
• 代理模型优化可以进行更多次迭代（500次 vs 50次）
• 代理模型找到的最优解经过FEM验证，结果可靠
• 对于复杂问题，代理模型可显著降低计算成本

---

7. 小结

本主题系统介绍了机器学习在结构动力学中的应用，包括：

7.1 方法对比

应用场景	推荐方法	优点	注意事项
响应预测	神经网络	高精度、非线性拟合强	需要大量训练数据
模态识别	随机森林/SVM	特征重要性可解释	特征工程关键
损伤检测	集成学习	鲁棒性好、处理不平衡数据	需要多工况数据
代理模型	神经网络/GP	预测快、适合优化	需要覆盖设计空间

7.2 关键要点

1. 数据质量：机器学习模型的性能高度依赖于训练数据的质量和数量
2. 特征工程：合适的特征提取是获得良好性能的关键
3. 模型选择：根据问题特点选择合适的机器学习算法
4. 验证评估：使用独立的测试集评估模型泛化能力
5. 物理约束：机器学习预测应满足物理约束和工程经验

7.3 工程应用建议

• 数据收集：建立系统的数据采集和标注流程
• 模型更新：定期使用新数据更新模型
• 不确定性量化：考虑模型预测的不确定性
• 混合方法：结合物理模型和数据驱动方法
• 实时应用：对于实时监测，需要优化模型推理速度

机器学习为结构动力学分析提供了新的工具和方法，与传统数值方法相结合，可以更好地解决复杂工程问题，实现智能化的结构健康监测和优化设计。

---