首页 > 编程语言 >SciPy优化C参数的梁结构频率响应曲线拟合

SciPy优化C参数的梁结构频率响应曲线拟合

来源：互联网 2026-06-30 08:00:07

在结构动力学分析中，当梁端存在弹性约束——例如螺纹连接或夹持界面——时，其等效刚度往往需要通过实验标定。这是一个常见且棘手的工程问题。本文介绍的是一套完整的技术路线：基于欧拉-伯努利梁理论与边界条件行列式模型，借助 scipy.optimize.differential_evolution 全局优化

在结构动力学分析中，当梁端存在弹性约束——例如螺纹连接或夹持界面——时，其等效刚度往往需要通过实验标定。这是一个常见且棘手的工程问题。本文介绍的是一套完整的技术路线：基于欧拉-伯努利梁理论与边界条件行列式模型，借助 scipy.optimize.differential_evolution 全局优化工具箱，对非线性刚度参数 c1–c4 进行系统寻优，最终使理论预测的一阶固有频率与实测频响数据精准匹配。

在振动建模中，一个核心场景是：梁端存在弹性约束，其平动刚度 \(k_t\) 和转动刚度 \(k_r\) 均非常数，而是由预紧力 \(p\) 调制的非线性函数。在该案例中，刚度模型由以下四参数函数描述：

\[k_n(p) = \frac{c_1 \tanh(c_2 p)}{1 + c_3 e^{-c_4 p}}, \quad k_t = c_t k_n, \quad k_r = c_r k_n\]

长期稳定更新的攒劲资源： >>>点此立即查看<<<

其中，\(c_t\) 和 \(c_r\) 是根据几何-材料常数预先确定的已知量，真正需要拟合的是 \(c_1\) 到 \(c_4\) 这四个核心参数。难点在于，理论频率 \(f\) 并非这些 \(c_i\) 的直接显函数——它通过求解特征方程 \(\det(\mathbf{A}(\omega; c_i)) = 0\) 隐式给出。这意味着梯度类方法（如 curve_fit 或 minimize）容易陷入局部极小值或遭遇雅可比矩阵奇异导致计算失败。针对此类问题，推荐采用无梯度、鲁棒性强的全局优化算法——differential_evolution。

以下是一套经过实际验证、可直接用于同类问题的完整优化流程。

步骤 1：重构目标函数——避免符号计算效率瓶颈

原代码中每次调用 CBmodel_T 都会触发 SymPy 符号微分和行列式展开，单次计算尚可接受，但在优化循环中效率极低。正确的做法是：预先推导并固化数值化的特征行列式表达式，确保其支持向量化运算。这样每次 det 求值均转为纯数值计算，速度显著提升：

import numpy as np
from scipy.optimize import differential_evolution
from scipy.optimize import root_scalar

ct = ((np.pi/2) * (1 - 0.3)) / (2 - 0.3)
cr = (1/4) * ((0.01848/2)**2 + 0.005**2) / (0.0067**2)

def det_numeric(beta, kt, kr):
    # 注意：此处为示意结构，实际应使用 sp.lambdify 从符号推导结果自动生成
    term1 = -(2 * beta**6 * np.cos(beta) * np.cosh(beta) - 2 * beta**6) * kt * kr
    term2 = (2 * beta**7 * np.sin(beta) * np.cosh(beta) - 2 * beta**7 * np.cos(beta) * np.sinh(beta)) * kt
    term3 = (2 * beta**9 * np.sin(beta) * np.cosh(beta) + 2 * beta**9 * np.cos(beta) * np.sinh(beta)) * kr
    term4 = 2 * beta**10 * (np.cos(beta) * np.cosh(beta) - 1)
    term5 = 2 * beta**8 * p_val * np.sin(beta) * np.sinh(beta)  # p_val 需外部传入
    return term1 + term2 + term3 + term4 + term5

def first_frequency_numeric(p, c1, c2, c3, c4):
    kn = (c1 * np.tanh(c2 * p)) / (1 + c3 * np.exp(-c4 * p))
    kt = ct * kn
    kr = cr * kn
    def char_eq(beta):
        return det_numeric(beta, kt, kr)
    try:
        res = root_scalar(lambda b: char_eq(b), bracket=[0.1, 10], method='brentq')
        beta1 = res.root
        w1 = beta1 ** 2
        scaling = np.sqrt((E * I) / (rho * A * L**4)) / (2 * np.pi)
        return w1 * scaling
    except Exception:
        return np.inf

pdata = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
fdata = np.array([31786.90305, 33344.70228, 34007.47212, 34388.51794, 
                  34640.69891, 34822.00471, 34959.70061, 35068.44216, 
                  35156.86945, 35230.43423])

步骤 2：定义优化目标函数——最小化残差平方和

目标函数用于计算理论一阶频率与实验值之间的误差，具体采用最小化均方误差（MSE）策略。考虑到计算过程中可能出现发散或负频率等无效解，需为这些情况设置足够大的惩罚值，引导优化器绕过无效区域：

def objective(c):
    c1, c2, c3, c4 = c
    f_pred = np.array([first_frequency_numeric(p, c1, c2, c3, c4) for p in pdata])
    if np.any(np.isnan(f_pred)) or np.any(f_pred <= 0):
        return 1e12
    return np.mean((f_pred - fdata) ** 2)

bounds = [
    (1e2, 1e5),   # c1: 等效法向刚度系数（N/m）
    (1e-1, 1e2),  # c2: 预紧力敏感度
    (1e-2, 1e2),  # c3: 指数修正幅值
    (1e-2, 1e2)   # c4: 指数衰减率
]

步骤 3：执行全局优化并验证结果

result = differential_evolution(
    objective,
    bounds,
    strategy='best1bin',
    maxiter=1500,
    popsize=20,
    tol=1e-4,
    seed=42,
    disp=True,
    workers=-1
)

if result.success:
    print(" 优化成功！最优参数：")
    print(f"  c1 = {result.x[0]:.3e}, c2 = {result.x[1]:.3e}")
    print(f"  c3 = {result.x[2]:.3e}, c4 = {result.x[3]:.3e}")
    print(f"  最小MSE = {result.fun:.6f}")
    f_fit = np.array([first_frequency_numeric(p, *result.x) for p in pdata])
    import matplotlib.pyplot as plt
    plt.figure(figsize=(8, 5))
    plt.scatter(pdata, fdata, label='实验数据', color='red', s=50)
    plt.plot(pdata, f_fit, 'b--', label='拟合曲线', linewidth=2)
    plt.xlabel('预紧力 $p$ (N)')
    plt.ylabel('一阶固有频率 $f_1$ (Hz)')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.title('C参数优化后的频率响应拟合')
    plt.show()
else:
    print(" 优化未收敛，请检查边界或目标函数鲁棒性。")

关键注意事项

行列式表达式必须零误差：det_numeric 中的公式必须与原始 SymPy 推导结果完全一致。推荐使用 sp.lambdify 从符号表达式自动生成，而非手动编写，以避免引入错误。
初始搜索范围需合理设定：边界过宽会显著增加计算时间，过窄则可能错过全局最优。建议先运行若干轮粗粒度网格搜索，逐步收窄边界范围。
妥善处理病态情况：root_scalar 可能无法找到根或报错。目标函数必须在此类情况下返回足够大的惩罚值（如 1e12），而非直接抛出异常，以免中断优化过程。
性能优化选项：若调用次数极高，可考虑使用 Numba JIT 预编译 first_frequency_numeric，或缓存常见 \(p\) 值下的 \(\beta\) 解。
不确定性评估不可缺失：优化完成后，建议在最优解附近使用 scipy.optimize.least_squares 进行局部拟合，获取参数的协方差矩阵和置信区间，以实现完整的参数标定。