建模问题学习笔记——现代辨识法之最小二乘法

一、最小二乘法

1、推导

含噪声的被辨识系统模型为：
$$y(k)=a_1^{0}y(k-1)+a_2^{0}y(k-2)+...+a_n^{0}y(k-n)+b_1^{0}u(k-1)+b_2^{0}u(k-2)+...+b_n^{0}u(k-n)+v(k)$$
$$\begin{gathered} y(k)=\left[\begin{array}{cccccccc}y(k-1)& y(k-2)& ...& y(k-n)& u(k-1)& u(k-2)& ...& u(k-n)\end{array}\right]\left[\begin{array}{c}{a}_1^0\\ a_2^0\\ ...\\ a_n^0\\ b_1^0\\ b_2^0\\ ...\\ b_n^0\end{array}\right]+v(k) \\ y(k)={\psi }^T(k){\theta }^0+v(k) \end{gathered}$$
通过递推可以得到模型输出应当为：
$$\begin{gathered} \left[\begin{array}{c}y(n)\\ y(n+1)\\ ...\\ y(N)\end{array}\right]=\left[\begin{array}{cccccccc}y(n-1)& y(n-2)& ...& y(0)& u(n-1)& u(n-2)& ...& u(0)\\ y(n)& y(n-1)& ...& y(1)& u(n)& u(n-1)& ...& u(1)\\ ...& ...& ...& ...& ...& ...& ...& ...\\ y(N-1)& y(N-2)& ...& y(N-n)& u(N-1)& u(N-2)& ...& u(N-n)\end{array}\right]\left[\begin{array}{c}{a}_1^0\\ a_2^0\\ ...\\ a_n^0\\ b_1^0\\ b_2^0\\ ...\\ b_n^0\end{array}\right]+\left[\begin{array}{c}v(n)\\ v(n+1)\\ ...\\ v(N)\end{array}\right] \\ Y(N)=\Psi (N){\theta }^0+V(N) \end{gathered}$$

实际估计的参数中不仅包含噪声项，还应当包含估计误差：
$$Y(N)=\Psi (N)\theta +\epsilon (N,\theta )$$

为了使误差最小化，定义性能指标：
$$J(\theta )=(Y-\Psi \theta {)}^T(Y-\Psi \theta )$$

对性能指标求导得到：
$${\hat{\theta }}_{LS}=({\Psi }^T\Psi {)}^{-1}{\Psi }^{T}Y$$

其应用需要两个条件：

1. N>3n-1。n是系统的阶次，N是数据的长度。
2. n阶持续激励条件。输入必须采用随机序列或M序列。

2、举例代码验证

对于系统$y(k)=a_{1}y(k-1)+a_{2}y(k-2)+b_{1}u(k-1)+b_{2}u(k-2)+\zeta (k)$，其输入为幅值为1的伪随机二位式序列，噪声是一个方差可调的标准正态分布的随机序列，真实值$a_1=1.5,a_2=-0.7,b_1=1.0,b_2=0.5$。

clc;clear all
%% 生成M序列 这里n=2，我们取数据长度N=100
L=100;
y1=1;y2=1;y3=0;
for i=1:L
    x1=xor(y2,y3);      %异或
    x2=y1;
    x3=y2;
    y(i)=y3;
    if y(i)>0.5
        u(i)=1;
    else
        u(i)=-1;
    end
    y1=x1;y2=x2;y3=x3;
end 
a1=1.5;a2=-0.7;b1=1.0;b2=0.5;
Y=[];
Fai=[];
for k=3:L
    v(k)=random('Normal',0,0.1); %正态分布的噪声
    y(k)=a1*y(k-1)+a2*y(k-2)+b1*u(k-1)+b2*u(k-2)+v(k);
    fai=[y(k-1),y(k-2),u(k-1),u(k-2)]';
    Fai=[Fai;fai'];
    Y=[Y;y(k)];
end
theta=inv(Fai'*Fai)*Fai'*Y;

用最小二乘法系统辨识出的结果为$a_1=1.4977,a_2=-0.6986,b_1=0.9937,b_2=0.5020$。

二、加权最小二乘法

1、推导

性能指标增加加权系数后，定义为：
$$J(\theta )=(Y-\Psi \theta {)}^{T}W(Y-\Psi \theta )$$

通常W取为对角矩阵，常用的指数加权函数是$w(k)=a{r}^{N-k}$。

对性能指标求导得到：
$${\hat{\theta }}_{WLS}=({\Psi }^{T}W\Psi {)}^{-1}{\Psi }^{T}WY$$

其应用同样需要：

1. N>3n-1。n是系统的阶次，N是数据的长度。
2. n阶持续激励条件。输入必须采用随机序列或M序列。

2、举例代码验证

clc;clear all
%% 生成M序列 这里n=2，我们取数据长度N=100
L=100;
y1=1;y2=1;y3=0;
for i=1:L
    x1=xor(y2,y3);      %异或
    x2=y1;
    x3=y2;
    y(i)=y3;
    if y(i)>0.5
        u(i)=1;
    else
        u(i)=-1;
    end
    y1=x1;y2=x2;y3=x3;
end 
a1=1.5;a2=-0.7;b1=1.0;b2=0.5;
Y=[];
Fai=[];
a=1.1;r=0.9;ww=[];
for k=3:L
    v(k)=random('Normal',0,0.1); %正态分布的噪声
    y(k)=a1*y(k-1)+a2*y(k-2)+b1*u(k-1)+b2*u(k-2)+v(k);
    fai=[y(k-1),y(k-2),u(k-1),u(k-2)]';
    Fai=[Fai;fai'];
    Y=[Y;y(k)];
    w=a*r^(L-k);
    ww=[ww w];
end
W=diag(ww);
theta=inv(Fai'*W*Fai)*Fai'*W*Y;

用加权最小二乘法系统辨识出的结果为$a_1=1.5045,a_2=-0.6969,b_1=1.0097,b_2=0.4989$。

三、递推最小二乘法

1、推导

当增加一对新的输入输出时，
$$\begin{gathered} P(N+1)=P(N)-\frac{P(N)\psi (N+1){\psi }^T(N+1)P(N)}{w^{-1}(N+1)+{\psi }^T(N+1)P(N)\psi (N+1)} \\ L(N+1)=P(N+1)\psi (N+1)w(N+1) \\ \Delta (N+1)=y(N+1)-{\psi }^T(N+1){\hat{\theta }}_{WLS}(N) \\ {\hat{\theta }}_{WLS}(N+1)={\hat{\theta }}_{WLS}(N)+L(N+1)\Delta (N+1) \end{gathered}$$

一般令${\hat{\theta }}_{WLS}(0)=0$，$P(0)=(1000,2000,...,5000)u_{0}I$，$u_0$为输入信号的峰值，辨识参数有几个P矩阵就有几维，注意小写希腊字母为向量，大写希腊字母为矩阵。

2、举例代码验证

clc;clear all
%% 生成M序列 这里n=2，我们取数据长度N=100
L=100;
y1=1;y2=1;y3=0;
for i=1:L
    x1=xor(y2,y3);      %异或
    x2=y1;
    x3=y2;
    y(i)=y3;
    if y(i)>0.5
        u(i)=1;
    else
        u(i)=-1;
    end
    y1=x1;y2=x2;y3=x3;
end 
a1=1.5;a2=-0.7;b1=1.0;b2=0.5;
Y=[];
Fai=[];
a=1.1;r=0.99;
P=1000*eye(4);
theta=[1;1;1;1];
thetaa=zeros(L,4);
thetaa(1,1:4)=theta';thetaa(2,1:4)=theta';
for k=3:L
    v(k)=random('Normal',0,0.1); %正态分布的噪声
    y(k)=a1*y(k-1)+a2*y(k-2)+b1*u(k-1)+b2*u(k-2)+v(k);
    fai=[y(k-1),y(k-2),u(k-1),u(k-2)]';
    
    Fai=[Fai;fai'];
    Y=[Y;y(k)];
    w=a*r^(L-k);
    
    P=P-P*fai*fai'*P/(w^(-1)+fai'*P*fai);
    LL=P*fai*w;
    zeta=y(k)-fai'*theta;
    theta=theta+LL*zeta;
    
    thetaa(k,:)=theta';
end
figure(1)
plot(1:L,thetaa(:,1),'r');
hold on
plot(1:L,thetaa(:,2),'b');
plot(1:L,thetaa(:,3),'g');
plot(1:L,thetaa(:,4),'k');
legend('a1','a2','b1','b2');

用递推加权最小二乘法系统辨识出的结果为$a_1=1.5009,a_2=-0.7072,b_1=0.9997,b_2=0.4915$。

本文是作者在日常学习生活中所作，难免有遗漏或错误，遇到问题的读者请点击给我写信向我的邮箱反馈，不胜感激。