【Hammerstein模型的级联】快速估计构成一连串哈默斯坦模型的结构元素研究（Matlab代码实现）

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

在许多振动应用中，所研究的系统略微非线性。Hammerstein模型的级联可以方便地描述这样的系统。Hammerstein提供了一种基于指数正弦扫描相位属性的简单方法。

构成一连串Hammerstein模型的结构元素可以在快速估计中起到关键的作用。Hammerstein模型由级联的非线性静态函数和线性动态函数组成。以下是研究Hammerstein模型结构元素的步骤：

1. 数据收集：首先，收集用于建立Hammerstein模型的数据。这些数据应包括系统的输入和输出信号，以便进行模型参数估计和验证。

2. 静态非线性函数选择：选择适当的静态非线性函数作为Hammerstein模型的非线性部分。常见的选择包括多项式函数、幂函数、指数函数、Sigmoid函数等。根据系统的特性和预期的非线性行为，选择最能表示系统的非线性特点的函数。

3. 参数估计：对选择的静态非线性函数进行参数估计。参数估计的方法可以根据函数的性质灵活选择，例如最小二乘法、最大似然估计法等。根据所选方法，使用数据集中的输入和输出信号优化非线性函数的参数。

4. 线性动态函数选择：选择适当的线性动态函数作为Hammerstein模型的动态部分。常见的选择包括传递函数、状态空间模型等。根据系统的动态特性，选择最适合描述系统响应的线性动态函数。

5. 参数估计：对选择的线性动态函数进行参数估计。使用数据集中的输入和输出信号，在模型的非线性部分和线性动态部分之间优化参数。

6. 模型验证：使用建立的Hammerstein模型对独立数据集进行验证。计算预测输出与真实输出之间的误差，评估模型的准确性和可靠性。如果有必要，可以对模型进行进一步调整和改进。

7. 性能分析：对Hammerstein模型的性能进行分析。例如，可以通过计算模型的拟合优度（如均方根误差）来评估模型的准确性。此外，还可以进行稳定性分析、系统辨识度评估等进一步分析。

需要注意的是，构建Hammerstein模型需要对非线性和线性组成部分的选择和参数估计进行适当的判断和调整。根据具体问题的复杂性和数据的可用性，可以采用各种方法和技术来加快估计和验证过程。

📚2 运行结果

 部分代码：

function hhat = Hammerstein_ID(input_sig,output,duration,f1,f2,fs,N,opt_meth,opt_filt)

%---------------------------------------------------------
%
% hhat = Hammerstein_ID(input_sig,output,f1,f2,fs,N,opt_meth,opt_filt)
%
% Estimates the Kernels "h" of the cascade of Hammerstein model of order N fed with
% the input signal "input" and where the corresponding output signal "output"
% has been recorded. "input" has to be an exponential sine sweep going from
% f1 to f2.
%
% Input parameters:
%   input_sig  : input exponential sine sweep
%   output : output of the system fed with the input signal
%   f1 : starting frequency of the sweep
%   f2 : end frequency of the sweep
%   fs : sampling frequency
%   N  : Order of the model to be identified
%   opt_meth : Method to use for the estimation (string expected)
%       - 'Reb': Method proposed by R閎illat et al. in [1]
%       - 'Nov': Method proposed by Novak et al in [2]
%   opt_filt : Specifies the method to use to compute the inverse filter 
%       (string expected). By default 'TFB_linear' is chosen.
%       - 'TFB_square': FTT based filter with a square window and
%         regularization (see [1])
%       - 'TFB_linear': FTT based filter with a square window with continuous
%         linear borders and regularization (see [1])
%       - 'TFB_gevrey': FTT based filter with a square window with infinitely
%         continuous gevrey borders and regularization (see [1])
%       - 'Nov' : Filter based on the analytical formulation using aymptotic 
%         signals (see [2]).
%
% Output:
%   h : 2D matrix containing the pseudo impulse responses (temporal domain) 
%   of the estimated kernels.

display('--> Hammerstein Identification in progress ...')

% Check arguments
if nargin<6
    display('    => ERROR : Incorrect number of arguments')
    return
elseif nargin<7
    display('    => No method option and filtering option specified. ')
    display('    => Method option = ''Reb'' chosen by default.')
    display('    => Filtering option = ''TFB_linear'' chosen by default.')
    opt_meth = 'Reb' ;
    opt_filt = 'TFB_linear';
elseif nargin<8
    display(['    => Method ' opt_meth ' chosen'])
    display('    => No filtering option specified. ')
    if strcmp(opt_meth,'Reb')
        opt_filt = 'TFB_linear';
        display('    => Filtering option = ''TFB_linear'' chosen by default.')
    elseif strcmp(opt_meth,'Nov')
        opt_filt = 'Nov';
        display('    => Filtering option = ''Nov'' chosen by default.')
    else 
        display('    => ERROR : Unknown method option')
        display('    => Select ''Reb'' or ''Nov''')
        return
    end
else
    
    if ( strcmp(opt_meth,'Nov') || strcmp(opt_meth,'Reb'))
        display(['    => Method ' opt_meth ' chosen'])
    else
        display('    => ERROR : Unknown method option')
        display('    => Select ''Reb'' or ''Nov''')
        return
    end
    
    if ( strcmp(opt_filt,'TFB_square') || strcmp(opt_filt,'TFB_linear') || strcmp(opt_filt,'TFB_gevrey') || strcmp(opt_filt,'Nov'))
        display(['    => Filtering ' opt_filt ' chosen'])
    else
        display('    => ERROR : Unknown filtering option')
        display('    => Select ''TFB_square'', ''TFB_linear'', ''TFB_gevrey'' or ''Nov''')
        return
    end
end

% Equivalent pulsations
w1 = f1/fs*2*pi;
w2 = f2/fs*2*pi;

% Convolution of the response with the inverse of the sweep
if strcmp(opt_meth,'Reb')
    inverse_input_sig = compute_inverse_filter(input_sig,f1,f2,fs,opt_filt) ;
    gToCut = convq(output,inverse_input_sig);
elseif strcmp(opt_meth,'Nov')
    % Nonlinear convolution in the spectral domain
    gToCut = nonlinear_convolution(output,duration,f1,f2,fs);        
    gToCut = [gToCut; gToCut];
end

% Computation of the delay of the pseudo RI
if strcmp(opt_meth,'Reb')
    T = length(input_sig); % Actual length of the sweep (in samples)
    deltaT = T*log(1:N)/log(w2/w1);
elseif strcmp(opt_meth,'Nov')
    T = length(output); % Actual length of the output (in samples)
    L = 1/f1*round( (duration*f1)/(log(f2/f1)) ); 
    deltaT = L*log(1:N)*fs;
end

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1] M. Rébillat, R. Hennequin, E. Corteel, B.F.G. Katz, "Identification of cascade of Hammerstein models for the description of non-linearities in vibrating devices", Journal of Sound and Vibration, Volume 330, Issue
5, Pages 1018-1038, February 2011.

[2] A. Novak, L. Simon, F. Kadlec, P. Lotton, "Nonlinear system identification using exponential swept-sine signal", IEEE Transactions on Instrumentation and Measurement, Volume 59, Issue 8, Pages 2220-2229, August 2010.

🌈4 Matlab代码实现