4. 迟滞数学模型（三）：Bouc-Wen模型及其参数辨识

好，咱们继续聊迟滞模型。前面讲了Preisach和Prandtl-Ishlinskii，这两个模型各有千秋，但都有一个共同点——它们属于“无记忆”的迟滞模型。什么意思呢？就是当前输出只取决于当前输入和历史极值，跟变化过程没关系。

但实际压电陶瓷不是这样的。你想想看，输入信号变化快慢不同，迟滞环的形状也会变。这就是所谓的“率相关”特性。那怎么办？这时候就该Bouc-Wen模型登场了。

4.1 Bouc-Wen模型长什么样？

Bouc-Wen模型最早是搞土木工程的人提出来的，用来描述建筑结构的滞回特性。后来被我们这帮做压电控制的人“借”过来用，发现效果还不错。

它的数学形式其实挺简洁的：

y(t) = k·u(t) + h(t)

其中h(t)满足：
dh/dt = A·du/dt - β·|du/dt|·|h|^(n-1)·h - γ·du/dt·|h|^n

这里y(t)是输出位移，u(t)是输入电压，k是线性增益系数，h(t)就是迟滞分量。A、β、γ、n是模型参数，决定了迟滞环的形状。

我个人习惯把Bouc-Wen看作一个“线性部分+非线性迟滞部分”的组合。线性部分好理解，就是压电陶瓷的压电效应；非线性部分嘛，就是那个让人头疼的迟滞。

4.2 参数到底有啥物理意义？

很多初学者问我：“这些参数看着好抽象，到底代表什么？”

我一般这么解释：

参数	物理意义	对迟滞环的影响
k	线性刚度系数	决定迟滞环的斜率
A	迟滞幅值系数	控制迟滞环的高度
β	迟滞形状参数	影响环的“胖瘦”
γ	迟滞平滑参数	决定环的尖锐程度
n	指数参数	控制从弹性到塑性过渡的平滑度

嗯，这里要注意：n一般取1或2。取1时模型相对简单，取2时能描述更复杂的迟滞行为。我在项目中通常先用n=1试试，效果不够好再调成2。

4.3 参数辨识怎么做？

参数辨识说白了就是“给模型找一套合适的参数，让它能拟合实验数据”。

我常用的方法分三步走：

1. 采集数据：给压电陶瓷施加一个正弦电压信号，记录输入电压和输出位移
2. 分离线性部分：用最小二乘法先估计k值
3. 优化非线性参数：用智能优化算法找A、β、γ、n的最优值

具体到第三步，我个人偏爱用粒子群优化（PSO）算法。为什么？因为它不需要求梯度，对初值不敏感，适合这种非线性参数辨识问题。

下面是我常用的MATLAB代码框架：

% Bouc-Wen模型参数辨识 - PSO算法
% 目标函数：最小化模型输出与实验数据的误差

function [params, error] = boucwen_identification(u_data, y_data, fs)
    % u_data: 输入电压数据
    % y_data: 输出位移数据  
    % fs: 采样频率
    
    % 第一步：估计线性增益k
    k_est = (y_data' * u_data) / (u_data' * u_data);
    
    % 第二步：PSO优化剩余参数
    lb = [0.1, 0.01, 0.01, 1];  % 参数下界 [A, beta, gamma, n]
    ub = [10, 10, 10, 2];       % 参数上界
    
    % PSO参数设置
    options = optimoptions('particleswarm', ...
        'SwarmSize', 50, ...
        'MaxIterations', 200, ...
        'Display', 'iter');
    
    % 定义目标函数
    fun = @(params) boucwen_error(params, u_data, y_data, k_est, fs);
    
    % 运行PSO
    [opt_params, error] = particleswarm(fun, 4, lb, ub, options);
    
    % 组合最终参数
    params = [k_est, opt_params];
end

function err = boucwen_error(params, u, y, k, fs)
    A = params(1); beta = params(2);
    gamma = params(3); n = params(4);
    
    % 用欧拉法求解微分方程
    dt = 1/fs;
    h = zeros(size(u));
    
    for i = 2:length(u)
        du = u(i) - u(i-1);
        dh = A*du - beta*abs(du)*abs(h(i-1))^(n-1)*h(i-1) ...
             - gamma*du*abs(h(i-1))^n;
        h(i) = h(i-1) + dh;
    end
    
    y_model = k*u + h;
    err = sum((y - y_model).^2);
end

4.4 避坑指南：我踩过的那些坑

做Bouc-Wen参数辨识，有几个坑我替你们踩过了：

• 初值敏感问题：我曾经用梯度下降法做优化，结果收敛到局部最优，迟滞环形状完全不对。后来改用PSO，问题就解决了。
• 过拟合陷阱：有一次我用高频噪声信号做辨识，模型参数拟合得特别好，但换到正弦信号就完全不能用。记住：训练数据要覆盖你实际工作频率范围。
• 数值稳定性：求解微分方程时，步长太大会导致数值发散。我建议采样频率至少是信号最高频率的100倍。

4.5 模型验证：怎么知道参数对不对？

参数辨识完了，怎么判断好坏？我一般看三个指标：

1. 均方根误差（RMSE）：小于位移幅值的2%就算合格
2. 迟滞环形状对比：肉眼观察模型输出和实验数据的迟滞环是否吻合
3. 泛化能力测试：用不同频率、不同幅值的信号验证模型

我个人最看重第三点。模型能拟合训练数据不算本事，能预测没见过的信号才算真功夫。

下面是我画的一个Bouc-Wen模型结构图，帮你理清思路：

从这张图你能看到，Bouc-Wen模型本质上是一个并联结构。输入信号同时走两条路：一条是线性通道，一条是带微分方程的非线性通道。两条路的输出加起来，就是最终的迟滞输出。

好了，Bouc-Wen模型就聊到这儿。这个模型虽然参数不多，但用好了能解决很多实际问题。下次你遇到率相关的迟滞问题，不妨试试它。

---

公众号:蓝海资料掘金营，微信deep3321