简介：HammersteinToolbox是一套面向MATLAB/Simulink环境的非线性系统建模工具箱，主要服务于需要对Hammerstein模型进行辨识、参数优化与仿真的工程师和研究人员。该资源共9个文件，包含8个m源码文件和1个txt许可协议文件，压缩包仅11KB，其中m文件覆盖模型识别、非线性卷积、逆滤波计算等核心函数及完整示例程序，结构轻量且便于二次开发。基于Simulink基础，工具箱提供从输入输出数据中自动估计静态非线性部分与线性动态部分的算法，支持最小二乘、梯度下降等优化策略，帮助用户在电力系统控制、声学信号处理等场景下快速搭建并验证非线性模型。借助演示程序和可视化工具，读者可以连贯掌握模型识别、参数整定、残差与阶跃响应分析的完整流程，从而提升系统建模精度与仿真效率。该资源已有562人学习下载，适合正在研究非线性系统辨识或Simulink建模的MATLAB开发者。 这篇博客我会直接从一个实际的建模问题切入，全程使用中文，按照从业者的口吻来展开，不带任何平台痕迹。
1. Hammerstein模型为什么值得单独做个工具箱
前阵子接了一个pH中和过程的建模需求，一开始我直接拿MATLAB里的System Identification Toolbox去跑，结果线性模型怎么调都差一口气。后来才意识到，问题不是出在算法上，而是线性框架本身就表达不了这类系统的核心特征。这类系统最典型的描述方式，就是Hammerstein模型——一个静态非线性环节串联一个线性动态环节。而网上现成的MATLAB实现要么太零散，要么只覆盖了最简单的辨识场景，所以我花了些时间把一套可复用的Hammerstein工具箱整理了出来，这篇内容就把它的模型原理、模块设计、辨识流程和踩坑细节一次讲清楚。
Hammerstein模型的结构用一句话概括：先经过一个无记忆的静态非线性函数，再进入一个线性动态系统。数学上可以表达为：
v(k) = f( u(k) ) y(k) = G(q) v(k) + e(k)其中u(k)是输入，v(k)是中间不可测变量，G(q)是离散传递函数，e(k)是噪声。因为v(k)在物理上通常没有对应的传感器，所以不能直接测、只能靠算法去估计，这正是该模型辨识困难的根本原因。
我可以用一个生活中的类比帮助理解：快餐厅的出餐流程是典型的Hammerstein结构。顾客点单（输入u）后，厨房先做一道“翻译”——把订单内容映射为具体的食材和份量（静态非线性f），然后进入后厨生产线（线性动态G），最后出餐（输出y）。中间的食材分配环节（v）你不会直接写在订单上，但它真实存在，而且决定了出餐口看到的实际品质。
1.1 哪些工程场景能对上号
Hammerstein模型的工业场景覆盖比很多人想象中要广。pH中和过程是最经典的代表，酸碱滴定曲线本身就是强非线性，加上混合搅拌罐的动态特性，正好是静态非线性加线性动态的结构。蒸馏塔的塔板温度控制、伺服电机带摩擦负载的运动控制、生物发酵过程的底物消耗和菌体生长，也都具备类似特征。
实际建模时，很多人会在两个极端之间纠结：用线性模型，结构简单但误差大；用神经网络黑箱，拟合能力强但缺乏可解释性，而且对样本量和训练过程都很敏感。Hammerstein模型处在两者之间——用事先确定的函数结构去描述非线性，用传递函数去描述动态特性，参数数量少、物理意义明确、泛化能力也有保障。这正是它值得被做成可复用工具箱的原因。
1.2 与Wiener模型的区别
在块结构非线性模型中，Hammerstein模型还有一个经常被混淆的“兄弟”——Wiener模型。两者的区别在于非线性环节的位置：Hammerstein是“非线性在前、线性在后”，Wiener是“线性在前、非线性在后”。这个顺序不仅仅影响数学表达，也直接决定了辨识算法的选择。下表是两者的核心对比：
| 对比维度 | Hammerstein模型 | Wiener模型 |
|---|---|---|
| 结构顺序 | 静态非线性 f → 线性动态 G | 线性动态 G → 静态非线性 f |
| 典型场景 | 执行器非线性（阀门、功放） | 传感器非线性（测量变送器） |
| 中间变量含义 | 非线性环节的输出 | 线性环节的输出 |
| 辨识难度 | 相对成熟，过参数化+LS | 需要迭代或特殊输出误差法 |
| 工具箱覆盖 | 本工具箱主攻 | 可作为后续扩展 |
2. 工具箱的模块划分与调用关系
这个工具箱在设计上遵循一个原则：让使用者用最少的代码完成“数据准备→模型辨识→结果验证”的闭环。我按功能将其划分为三个模块层次，分别对应模型对象层、数据层和算法层。
2.1 模型对象层：把非线性环节和线性环节绑在一起
模型对象是整个工具箱的数据核心。非线性环节和线性环节原本是两块独立的数据结构，如果不加封装，辨识过程中很容易出现参数错配。我设计了一个结构体HammersteinModel，把两者绑定在一个对象里：
model = struct(); model.nl_type = 'poly'; % 非线性类型：'poly','spline','deadzone' model.nl_order = 3; % 多项式阶数 model.nl_coef = []; % 非线性环节系数，辨识后填充 model.G_num = []; % 线性环节分子系数 model.G_den = []; % 线性环节分母系数 model.Ts = 0.1; % 采样周期 model.fit_time = datetime('now'); % 辨识时间戳把模型配置和数据分离，在调参时非常方便。你可以在不改变数据的情况下反复修改模型结构参数，对比不同配置的辨识效果，而不需要复制数据、移动文件等。
2.2 数据层：激励信号生成才是辨识成败的一半
很多人做系统辨识时只关心辨识算法，忽略了激励信号设计。实际上，输入信号的品质直接决定了参数的可辨识性，在Hammerstein模型中尤其明显。
工具箱提供了三种默认激励信号：PRBS（伪随机二进制序列）、随机幅值序列、多正弦叠加。我第一次使用工具箱做仿真验证时，就用PRBS信号，结果非线性环节始终辨识不准，后来换成随机幅值序列才解决问题。原因在于PRBS只有两个幅值电平，最多激励非线性函数上的两个点，中间区域的特性完全没被激活。
% 生成激励信号 u1 = htSignalGen('prbs', 1000, 'levels', [-1, 1]); % 两电平PRBS u2 = htSignalGen('randamp', 1000, 'range', [-1, 1], 'levels', 9); % 9幅值随机序列2.3 辨识算法层：三种估计算法各管什么场景
辨识算法层是整个工具箱中代码量最大的部分，内置了三种估计算法：
过参数化最小二乘：实现简单、计算速度快，适用于噪声较小、模型阶次已知的场景，原理上先将非线性环节用基函数展开，再和线性环节参数合并成回归问题。
交替最小二乘迭代法：固定非线性参数求线性参数，再固定线性参数求非线性参数，往复迭代。适用于中等噪声、结构相对简单的场景，但收敛性依赖初值。
基于优化工具箱的非线性最小二乘：通过lsqnonlin同时优化全部参数，适用于小规模问题、精度要求高的场景，但初始化敏感、计算量最大。
实际使用中的选择建议很简单：先快速用第一种方法获得一个初值，再用第二种方法精修，如果精度还不够，才动用第三种方法。这三种算法的组合可以覆盖绝大多数Hammerstein辨识需求。
3. 跑通一个完整的辨识流程
接下来的内容是一个可以直接在MATLAB中复现的完整案例。我构造一个已知的Hammerstein系统，用工具箱辨识它，随机验证工具箱的辨识效果。这个流程也是普通用户最常用的路径。
3.1 第一步：生成仿真数据
构造一个标准的Hammerstein系统：非线性环节为f(u)=tanh(2u)，线性环节为离散传递函数G(z)=0.2z/(z^2-1.5z+0.7)，采样周期0.1秒。在此基础上叠加高斯白噪声，信噪比设为30dB。
% 步骤1：构建真实系统对象 trueModel = htCreateModel('poly', 9, [0, 2, 8/3, 0, 0, 0, 0, 0, 0], [0.2, 0], [1, -1.5, 0.7], 0.1); % 步骤2：生成激励信号并仿真 u = htSignalGen('randamp', 2000, 'range', [-1.2, 1.2], 'levels', 11); [y, v] = htsimulate(trueModel, u); % 步骤3：叠加测量噪声 rng(2024); y_noisy = y + 0.02 * randn(size(y));需要特别强调的是，这里的非线性环节使用了阶数为9的多项式，系数中只有三项非零，其余项为0。这样设计是为了验证工具箱能否在允许冗余参数的情况下辨识出真实的稀疏结构——实际工程中我们并不知道非线性函数的真实形式，所以工具箱必须有足够的表达容量。
3.2 第二步：配置模型结构并执行辨识
辨识的关键是给算法一个合理的先验配置。我选择了和真实系统相同的线性阶次（分母3阶、分子2阶），非线性基函数选择多项式且最高阶数设为9。如果对模型结构一无所知，可以通过后文第四节的方法进行阶次选择。
% 配置待辨识模型结构 config = struct(); config.nl_type = 'poly'; config.nl_order = 9; config.G_num_order = 2; config.G_den_order = 3; config.Ts = 0.1; % 执行辨识 estModel = htIdentify(u, y_noisy, config, 'method', 'alternating');htIdentify是工具箱的主入口函数，内部会根据method参数分派到具体的算法实现。这里选择alternating（交替最小二乘迭代法），因为它对噪声有中等程度的鲁棒性，并且不要求噪声必须是白噪声。
3.3 第三步：验证结果并判定模型质量
辨识结束后，必须用独立于辨识数据的验证集来检验模型质量。直观的方法是直接对比真实系统输出和模型预测输出，但更定量化的指标是BFR（Best Fit Rate，最佳拟合率）：
% 用新激励信号生成验证数据 u_val = htSignalGen('prbs', 1000, 'levels', [-1, 1]); y_val_true = htsimulate(trueModel, u_val); y_val_model = htsimulate(estModel, u_val); % 计算BFR拟合指标 BFR = max(0, 1 - norm(y_val_true - y_val_model) / norm(y_val_true - mean(y_val_true))); fprintf('BFR = %.2f%%\n', BFR * 100);在一次典型运行中，这个流程得到的BFR约为85%到92%。非线性环节估计结果和真实函数tanh(2u)在u的取值范围内几乎重叠，只在两端边界处有微小偏差，这是多项式在边界处的固有问题，不影响整体模型质量。
4. 算法实现里的关键细节
工具箱中的各个算法模块，单独拿出来都不算特别高深，但把它们组合成一个稳定、易用的工具箱，需要处理很多教科书上没有的工程细节。
4.1 过参数化怎么避免数值病态
过参数化方法的基本思路，是将多项式基函数展开后的各分量，分别通过线性动态系统，然后再叠加。以多项式非线性为例，模型可以写成：
y(k) = a1·G(q)·u(k) + a2·G(q)·u²(k) + ... + an·G(q)·uⁿ(k) + e(k)每个基函数的输出被同一个线性子系统滤波，而各基函数之间的相关性会导致信息矩阵近似奇异。这正是数值病态的根源。我在工具箱中采用了两种手段来处理：
一是输入信号预处理。在构造基函数之前，先把输入u归一化到[-1,1]区间，避免高次项数值过大。二是在最小二乘求解时加入Tikhonov正则化项，惩罚过大的参数值。这两种手段组合使用后，即使用到9阶多项式也不会出现数值发散。
4.2 迭代法收敛控制与初值策略
交替最小二乘虽然实现简单，但对初值比较敏感，容易陷入局部极小值。我的策略是：任何一个辨识任务，都先用过参数化方法快速得到一个粗糙解，用它来初始化迭代法，而不是随机生成初值。这个策略大幅提升了收敛到全局最优的概率。
迭代过程的停止条件也很关键，单纯设定最大迭代次数会导致次优解，单纯设定参数变化阈值又可能在早期就过早停止。工具箱采用双条件组合：参数变化相对值小于1e-6且连续保持5次迭代，才判定收敛。这避免了偶然的震荡造成的误判。
4.3 结构参数怎么选：阶次与基函数
线性部分阶次可以通过经典的AIC/BIC准则在工具箱中自动搜索，但非线性部分的阶次选择更微妙——它直接对拟合精度产生影响，阶次太低欠拟合，阶次太高不仅带来过拟合，还会引发数值病态。
我的实操建议是：不要试图在一次辨识中确定所有参数。先用低阶非线性多项式（2到3阶）跑通流程，看残差的非线性特征；如果残差中仍有明显的非线性趋势，再逐步提高阶次。交叉验证是判断阶次是否合适的最终标准：训练集上拟合度持续上升、验证集上拟合度开始下降，这个拐点往往就是合适的阶次。
基函数的选择上，我的默认推荐是多项式。它实现简单、导数连续、易于解析处理，但在动态范围较大的场景中，样条基函数比多项式更稳定，只是需要额外设置节点位置。工具箱预留了spline选项，使用augknt函数生成均匀节点。
5. 踩坑记录与排查思路
工具箱在反复使用的过程中，我积累了一些典型的失败案例，记录在这里，避免你再走一遍弯路。
5.1 激励信号幅值太小，非线性环节辨识不出来
这个问题的典型现象是：辨识结果在训练集上和真实输出几乎重合，但把非线性环节单独画出来，和真实函数相差很大。第一次遇到时我很困惑，因为整体模型拟合度很高，为什么会这样？
检查发现，PRBS信号的两电平幅值分别为-1和1，非线性函数在[-1,1]范围内确实有输出，但对于强非线性函数，关键的特征区域往往集中在幅值更大的区间。两电平激励只能激活两个点，无法还原中间区域的曲线形状。解决办法是改用多幅值随机序列，确保输入信号的电平数覆盖非线性函数的特征区间。信号幅值范围需要覆盖实际工作区间，不要贸然超过设备允许范围。
5.2 多项式阶数一高就发散
多项式阶数高到一定程度后，辨识结果开始出现大幅震荡，甚至参数值数量级达到10的10次方以上。这是典型的高阶多项式数值病态问题。基函数u、u²、u³等在高阶时数值差异极大，导致信息矩阵条件数爆炸。
解决方案有两层。第一层是对输入做归一化处理，将u缩放到[-1,1]区间，这个操作能显著改善条件数。第二层是改用正交多项式基函数，工具箱中预留了legendre选项，使用勒让德多项式作为基函数，在线性最小二乘求解时可显著降低条件数。如果两种方法都用了仍然发散，那就是阶次过高，应该降低阶次。
5.3 模型拟合好但验证集一塌糊涂
这是过拟合的典型表现。训练集BFR高达95%，但换一批数据之后BFR直接掉到50%以下，说明模型记住的是训练数据中的噪声细节，而不是系统的真实动态。
排查步骤是：先检查线性部分阶次是否过高，再利用交叉验证绘制“阶次-拟合度”曲线，观察训练集和验证集的分叉点。使用工具函数htModelSelection自动完成这个过程更好——它会循环不同阶次组合，计算BFR值，并输出推荐配置。从经验来看，Hammerstein模型的实际应用中，线性部分阶次超过5阶、非线性部分超过7阶，都属于高风险配置，除非有充分的物理依据，否则不建议使用。
5.4 有色噪声下的有偏估计问题
最后再提一个辨识领域的经典问题：当噪声是有色噪声而非白噪声时，普通最小二乘估计是有偏的。很多使用者忽略了这个前提条件，导致参数估计总差一点点。
工具箱目前的默认实现假设噪声是白噪声。如果遇到有色噪声场景，建议在辨识前先对数据做预白化滤波，或者改用辅助变量法。这个扩展也在工具箱的后续计划中——处理Hammerstein模型的辅助变量辨识，是把工具箱推向工业级别应用的关键一步。
我在实际使用这个工具箱时最大的体会是：工具箱最大的价值不在于它的算法有多高明，而在于它把“中间变量不可测”这个Hammerstein辨识的核心难点封装干净了。使用者不需要对付数学推导，只需要明白“数据质量决定辨识上限，算法只是逼近这个上限”这一件事，就能把工具用好。后续我打算继续扩展多变量Hammerstein模型的支持，以及Hammerstein-Wiener串联结构，如果你在实际使用中也碰到了新的坑点，欢迎交流，我来更新踩坑记录。
本文还有配套的精品资源，点击获取
