基于最小二乘法和最大似然估计法的系统参数辨识MATLAB仿真

目录

一、理论基础

1.1、最小二乘法

1.2、最大似然估计法

1.3、基于最小二乘法和最大似然估计法的系统参数辨识理论研究

二、MATLAB仿真程序

三、仿真结果

---

一、理论基础

        系统辨识源于工业工程控制，随着计算机技术的发展，出现了许多辨识软件可以用于辅助辨识理论的研究。它不仅在工业领域中有着广泛的应用，而且有经济、社会和环境也有重要的应用价值。1）用于控制系统的设计和分析2）用于在线控制3）用于预报预测4）用于监视系统参数并实现故障诊断。
       系统辨识是研究建立系统数学模型的一种理论和方法。所谓辨识就是从含有噪声的输入输出数据中提取被研究对象的数学模型。一般来说，这个模型只是对象的输入输出特性在某种准则意义下的一种近似，挖的程度取决于人们对系统先验知识的认识深化程度和对数据集合性质的了解，以及所选用的辨识方法是否合理。或者说，辨识技术帮助人们在表征被研究系统对象、现象或过程的复杂因果关系时，尽可能准确地确立它们之间的定量依存关系。基于辨识原理和方法的基础上，需要设计如何将辨识流程转换成计算机实现。这就需要仔细分析辨识流程的总体结构以及其中的各个步骤的先后顺序，互相衔接等。在采用离线或在线算法时，所对应的流程也是完全不同的。
       利用系统辨识方法建立系统数学模型，不仅是由于需要提高系统的性能或控制、增加对系统的认识或仿真，而且也反映出发展与需要。例如，社会与环境系统、经济与企业系统、生物与生态系统和其他高精工程系统的需要是现代系统辨识得到发展的动力。教学与优化理论的引入用计算机的调整发展，给予系统辨识理论与应用注入了新的活力。

       系统参数辨识是一种通过对系统输入和输出数据进行处理，估计系统模型参数的方法。最小二乘法和最大似然估计法是常用的系统参数辨识方法的基本原理如下：

1.1、最小二乘法

       最小二乘法是一种通过最小化预测输出与实际输出之间的误差平方和来估计系统模型参数的方法。假设系统模型的输出为y(t)，实际输出为d(t)，预测输出为y^(t)，则误差e(t)可以表示为：

e(t) = d(t) - y^(t)

我们的目标是找到一组模型参数，使得误差平方和J最小，即：

J = ∑[e(t)]^2 = ∑[d(t) - y^(t)]^2

       为了使J最小，我们可以采用梯度下降法或其他优化算法来求解模型参数。在最小二乘法中，我们通常假设系统模型是线性的，即：

y^(t) = ∑[phi(t,i) * θ(i)] + e(t)

其中，phi(t,i)是基函数，θ(i)是模型参数，e(t)是噪声。将上式代入误差平方和J中，可以得到：

J = ∑[d(t) - ∑[phi(t,i) * θ(i)]]^2

       为了最小化J，我们可以对其求导，并令导数为0，从而得到模型参数的估计值。具体地，我们可以采用批量最小二乘法或递推最小二乘法来估计模型参数。

1.2、最大似然估计法

       最大似然估计法是一种通过最大化似然函数来估计系统模型参数的方法。假设系统模型的输出为y(t)，实际输出为d(t)，则似然函数L可以表示为：

L = ∏[p(d(t)|y(t))]

我们的目标是找到一组模型参数，使得似然函数L最大，即：

argmax[L] = argmax[∏[p(d(t)|y(t))]]

       为了使L最大，我们可以采用梯度上升法或其他优化算法来求解模型参数。在最大似然估计法中，我们通常假设噪声是服从某种概率分布的，例如高斯分布。假设噪声e(t)服从均值为0、方差为σ^2的高斯分布，则似然函数L可以表示为：

L = ∏[(1/(σ√(2π))) * exp(-(e(t)^2)/(2σ^2))]

其中，e(t) = d(t) - y^(t)。将上式取对数，可以得到：

ln[L] = -N/2 * ln(2πσ^2) - ∑[(e(t)^2)/(2σ^2)]

         其中，N是数据长度。为了使ln[L]最大，我们可以对其求导，并令导数为0，从而得到模型参数的估计值。具体地，我们可以采用批量最大似然估计法或递推最大似然估计法来估计模型参数。需要注意的是，在最大似然估计法中，我们需要知道噪声的概率分布及其参数，这通常需要对系统进行一些先验知识或假设。

1.3、基于最小二乘法和最大似然估计法的系统参数辨识理论研究

        辨识就是在输入和输出数据的基础上，从一组给定的模型类中，确定一个与所测系统等价的模型。辨识的三大要素是输入和输出数据、模型类和等价准则。其中，数据是辨识的基础；准则是辨识的优化目标；模型类是寻找模型的范围。辨识有三个要素，数据、模型类和准则。辨识就是按照一个准则在一模型类中选择一个与数据拟合得最好的模型。辨识的实质就是从一模型类中选择一个模型，按照某种准则，使其能最好地拟合所关心的实际过程的动态特性。

        辨识的方法多种多样，在经典辨识方法中有频率响应法、脉冲相应法等，得到的模型是非参数模型，通常是系统在时域或频域的一个数值描述。另一类是最小二乘类辨识方法，在这种辨识方法中，首先给出模型的结构，在结构框架下确定系统模型的最优参数。这类具有格式规范的辨识方法可以演绎成递推形式，已经渐渐替代了经典的辨识算法，成为多数研究者所采用的方法。 下面给出一种最小二乘类辨识的模型： 系统的输入输出关系可写成如下的最小二乘格式： 

式中，adm和adg是校正因子。 最小二乘原理可适用于许多模型类，一般的模型格式可写成： 

拟使用MATLAB主要以最小二乘法为例进行系统辨识仿真。

二、MATLAB仿真程序

function [Cs,Er] = func_RML_system_identification(Dout,Din,L);

c1      = [0.0001,0.0001,0.0001,0.0001,0.0001,0.0001]';
p0      = eye(6,6);       
Doutf(1)= 0.1;
Doutf(2)= 0.1;
Dinf(2) = 0.1;
Dinf(1) = 0.1;
vf(2)   = 0.1;
vf(1)   = 0.1;
v       = zeros(1,L);
a1      = zeros(1,L);
a2      = zeros(1,L);
b1      = zeros(1,L);
b2      = zeros(1,L);
ea1     = zeros(1,L);
ea2     = zeros(1,L);
eb1     = zeros(1,L);
eb2     = zeros(1,L);
for k=3:L
    h   = [-Dout(k-1);-Dout(k-2);Din(k-1);Din(k-2);v(k-1);v(k-2)];
    hf  = h;
    K   = p0*hf*inv(hf'*p0*hf+1); 
    p   = [eye(6,6)-K*hf']*p0; 
    v(k)= Dout(k)-h'*c1;
    c   = c1+K*v(k) ;
    p0  = p;
    c1  = c;
    a1(k)=c(1);
    a2(k)=c(2);
    b1(k)=c(3);
    b2(k)=c(4);
    ea1(k)=abs(a1(k)+1.5);
    ea2(k)=abs(a2(k)-0.7);
    eb1(k)=abs(b1(k)-1.0);
    eb2(k)=abs(b2(k)-0.5);
    Doutf(k) = Dout(k)-c(5)*Doutf(k-1)-c(6)*Doutf(k-2);
    Dinf(k)  = Din(k)-c(5)*Dinf(k-1)-c(6)*Dinf(k-2);    
    vf(k)    = v(k)-c(5)*vf(k-1)-c(6)*vf(k-2);  
    hf=[-Doutf(k-1);-Doutf(k-2);Dinf(k-1);Dinf(k-2);vf(k-1);vf(k-2)];    
end 

Cs = [a1;a2;b1;b2];
Er = [ea1;ea2;eb1;eb2];

function [c,e] = func_Least_squares_system_identification(Dout,Din,L);

%用最小二乘算法辨识参数
%被辨识参数的初始值采用直接取方式
c0 = [0.00001 0.00001 0.00001 0.00001]'; 
p0 = 10^7*eye(4,4); 
%被辨识参数矩阵的初始值及大小
c  = [c0,zeros(4,L-1)];  
%相对误差的初始值及大小
e  = zeros(4,L); 
n  = 0;
for k=3:L; %开始递推运算
    k
    %求h(k)
    hk    = [-Dout(k-1),-Dout(k-2),Din(k-1),Din(k-2)]';  
    %求K(k)
    K     = p0*hk*inv(hk'*p0*hk+1); 
    %求θ(k)
    c1    = c0+K*[Dout(k)-hk'*c0]; 
    %求相对误差
    ee    =(K*(Dout(k)-hk'*c0)); 
    %新获得的参数作为下一次递推的旧参数
    e(:,k)= ee;  
    c0    = c1; 
    c(:,k)= c1;
    %求p(k)值
    pk    = p0-0.5*K*K'*[hk'*p0*hk+1]; 
    %把当前的p(k)值给下次用
    p0    = pk; 
    n     = n+1;
end

三、仿真结果

A27-05