极大似然参数估计法

ayuan0211

已于 2022-11-25 17:43:02 修改

阅读量2k

点赞数 3

分类专栏：基于MATLAB控制系统辨识系列文章标签：算法 matlab 学习

于 2022-11-25 17:42:50 首次发布

本文链接：https://blog.csdn.net/weixin_38342580/article/details/128007289

版权

基于MATLAB控制系统辨识系列专栏收录该内容

8 篇文章 6 订阅

订阅专栏

一、原理

极大似然参数估计法需要构造一个以观测数据和未知参数为自变量的似然函数，使这个函数达到极大参数值，就是模型的参数估计值。

通常噪声的概率密度函数作为似然函数，所以极大似然函数法需要已知噪声的分布。在最简单的情况下，可假定噪声具有正态分布。

优点：具有良好的渐进性质

缺点：计算量大

考虑控制系统模型简化为CARMA模型：

$A(z^{-1})y(k)=B(z^{-1})u(k)+C(z^{-1})\xi(k)$

则递推极大似然参数估计算法公式为：

$\left\{\begin{matrix} \hat{\theta}(k)= \hat{\theta}(k-1)+K(k) \hat{\xi}(k) \\K(k)=P(k-1)\phi_f(k)[1+\phi_f^T(k)P(k-1)\phi_f(k)]^{-1} \\P(k)=[I-K(k)\phi_f^T(k)]P(k-1) \\ \hat\xi(k)=y(k)-\phi^T(k)\hat{\xi}(k-1) \\ \phi(k)=[-y(k-1),\dots,-y(k-n_a),u(k-d),\dots,u(k-d-n_b),\hat\xi(k-1),\dots,\hat\xi(k-n_c)]^T \\ \phi_f(k)=[-y_f(k-1),\dots,-y_f(k-n_a),u_f(k-d),\dots,u_f(k-d-n_b),\hat\xi_f(k-1),\dots,\hat\xi_f(k-n_c)]^T \\ y_f(k)=y(k) -\hat c_1(k)y_f(k-1)-\dots-\hat c_{n_c}(k)y_f(k-n_c) \\ u_f(k)=u(k) -\hat c_1(k)u_f(k-1)-\dots-\hat c_{n_c}(k)u_f(k-n_c) \\ \hat\xi_f(k)=\hat\xi(k) -\hat c_1(k)\hat\xi_f(k-1)-\dots-\hat c_{n_c}(k)\hat\xi_f(k-n_c) \end{matrix}\right.$

此方法初值一般取 $\hat\theta(0)=\vec{0}$ 、 $P(0)=I$

算法步骤：

Setp1:设置初值 $\hat\theta(0)=\vec{0}$ 、 $P(0)=I$ ，输入初始数据；

Step2：采样当前输出y(k)和输入u(k);

Step3：计算 $\hat\xi(k)$ ,并构造 $\phi_f(k)$ .

Step4：利用上述公式，计算K(k）、 $\hat\xi(k)$ 和P(k);

Step5：计算 $y_f(k)$ 、 $u_f(k)$ 、 $\hat\xi_f(k)$ ;

Step6: $k\rightarrow k+1$ ，返回Step2，继续循环。

二、程序代码

系统为：

$y(k)-1.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+\xi(k)-0.5\xi(k-1)$

式中， $\xi(k)$ 是方差为1的白噪声。

取初值为 $\hat\theta(0)=\vec{0}$ 、 $P(0)=I$

选择方差为1的白噪声作为输入信号u(k),采用极大似然估计法进行参数估计。

%递推极大似然参数估计

a=[1 -1.5 0.7]';
b=[1 0.5]';
c=[1 -0.5]';
d=1;%对象参数

na=length(a)-1;nb=length(b)-1;nc=length(c)-1;%na、nb、nc为A、B阶次
nn=max(na,nc);%用于yf(k-i) uf(k-i)更新
L=1000;%仿真长度
uk=zeros(d+nb,1);%输入初值：uk(i)表示u(k-i)
yk=zeros(na,1);%输出初值
xik=zeros(nc,1);%白噪声初值 \xi
xiek=zeros(nc,1);%白噪声估计初值
yfk=zeros(nn,1);%yf(k-i)
ufk=zeros(nn,1);%uk(k-i)
xiefk=zeros(nc,1);%\xif(k-i)
u=sqrt(1.1)*randn(L,1);%输入采用白噪声序列
xi=sqrt(1)*randn(L,1);%白噪声序列 \xi

thetae_1=zeros(na+nb+1+nc,1);%参数估计初值
P=eye(na+nb+1+nc);

for k=1:L
    y(k)=-a(2:na+1)'*yk+b'*uk(d:d+nb)+c'*[xi(k);xik];
    
    %构造向量
    phi=[-yk;uk(d:d+nb);xiek];
    xie=y(k)-phi'*thetae_1;
    phif=[-yfk(1:na);ufk(d:d+nb);xiefk];
    
    %递推极大似然参数估计算法
    
    K=P*phif/(1+phif'*P*phif);
    thetae(:,k)=thetae_1+K*xie;
    P=(eye(na+nb+1+nc)-K*phif')*P;
    
    yf=y(k)-thetae(na+nb+2:na+nb+1+nc,k)'*yfk(1:nc);%yfk
    uf=u(k)-thetae(na+nb+2:na+nb+1+nc,k)'*ufk(1:nc);%ufk
    xief=xie-thetae(na+nb+2:na+nb+1+nc,k)'*xiefk(1:nc);%xiefk
    
    
    %更新数据
    thetae_1=thetae(:,k);
    
    for i=d+nb:-1:2
        uk(i)=uk(i-1);
    end
    uk(1)=u(k);
    
    for i=na:-1:2
        yk(i)=yk(i-1);
    end
    yk(1)=y(k);
    
    for i=nc:-1:2
        xik(i)=xik(i-1);
        xiek(i)=xiek(i-1);
        xiefk(i)=xiefk(i-1);
    end
    xie(1)=xi(k);
    xiek(1)=xie;
    xiefk(1)=xief;
    
    for i=nn:-1:2
        yfk(i)=yfk(i-1);
        ufk(i)=ufk(i-1);
    end
    yfk(1)=yf;
    ufk(1)=uf;     

end

figure(1)
plot([1:L],thetae(1:na,:),[1:L],thetae(na+nb+2:na+nb+1+nc,:));
xlabel('k');
ylabel('参数估计a、c');
legend('a_1','a_2','c_1');
axis([0 L -2 2]);
figure(2)
plot([1:L],thetae(na+1:na+nb+1,:));
xlabel('k');
ylabel('参数估计b');
legend('b_0','b_1');
axis([0 L 0 1.5]);

三、运行结果

a1=-1.49296309858392
a2=0.699476918865982
b0=1.02141798797169
b1=0.481829921238559
c1=0.00780319291205411

估计参数与理论值有一定差距，系统辨识认为误差在10%以内都是良好的模型复现。

附录：名词解释

1.似然函数：

假设一个袋子里有若干标号1、2的小球，从袋子中取出1的概率为 $\theta$ ，取出2的概率为 $1-\theta$

X	1	2
P	$\theta$	$1-\theta$

假设抽取5次，结果为 1 1 2 1 2

那么似然函数就是这5次结果对应概率的乘积即：

$\L(\theta)=\theta^3(1-\theta)^2$

最大似然估计则为函数L取最大时， $\theta$ 的值。

降低次数，取ln，得 $lnL(\theta)=3ln(\theta)+2ln(1-\theta)$

求关于 $\theta$ 的导数。

$\frac{Ln(\theta)}{d\theta}=\frac{3}{\theta}-\frac{2}{1-\theta}=0\rightarrow \hat{\theta}=\frac{3}{5}$

由于样本量少，所以 $\theta$ 与理论值有差别，这里只举例。

本文参考：1.《系统辨识与自适应控制MATLAB仿真》庞中华崔红

2. 十分钟搞定最大似然估计_哔哩哔哩_bilibili

ayuan0211

关注

3
点赞
踩
25

收藏

觉得还不错? 一键收藏
打赏
0
评论
极大似然参数估计法

极大似然参数估计法需要构造一个以观测数据和未知参数为自变量的似然函数，使这个函数达到极大参数值，就是模型的参数估计值。通常噪声的概率密度函数作为似然函数，所以极大似然函数法需要已知噪声的分布。在最简单的情况下，可假定噪声具有正态分布。优点：具有良好的渐进性质
复制链接

扫一扫