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

一、最小二乘法

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 ) y(k)=a_1^0y(k-1)+a_2^0y(k-2)+...+a_n^0y(k-n)+b_1^0u(k-1)+b_2^0u(k-2)+...+b_n^0u(k-n)+v(k) y(k)=a10y(k1)+a20y(k2)+...+an0y(kn)+b10u(k1)+b20u(k2)+...+bn0u(kn)+v(k)
y ( k ) = [ y ( k − 1 ) y ( k − 2 ) . . . y ( k − n ) u ( k − 1 ) u ( k − 2 ) . . . u ( k − n ) ] [ a 1 0 a 2 0 . . . a n 0 b 1 0 b 2 0 . . . b n 0 ] + v ( k ) y ( k ) = ψ T ( k ) θ 0 + v ( k ) y(k)=\begin{bmatrix}y(k-1)&y(k-2)&...&y(k-n)&u(k-1)&u(k-2)&...&u(k-n)\end{bmatrix}\begin{bmatrix} a_1^0\\a_2^0\\...\\a_n^0\\b_1^0\\b_2^0\\...\\b_n^0 \end{bmatrix}+v(k)\\y(k)=\psi^T(k)\theta^0+v(k) y(k)=[y(k1)y(k2)...y(kn)u(k1)u(k2)...u(kn)]a10a20...an0b10b20...bn0+v(k)y(k)=ψT(k)θ0+v(k)
通过递推可以得到模型输出应当为:
[ y ( n ) y ( n + 1 ) . . . y ( N ) ] = [ 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 ) ] [ a 1 0 a 2 0 . . . a n 0 b 1 0 b 2 0 . . . b n 0 ] + [ v ( n ) v ( n + 1 ) . . . v ( N ) ] Y ( N ) = Ψ ( N ) θ 0 + V ( N ) \begin{bmatrix} y(n)\\y(n+1)\\...\\y(N) \end{bmatrix}=\begin{bmatrix} 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{bmatrix}\begin{bmatrix} a_1^0\\a_2^0\\...\\a_n^0\\b_1^0\\b_2^0\\...\\b_n^0 \end{bmatrix}+\begin{bmatrix} v(n)\\v(n+1)\\...\\v(N) \end{bmatrix}\\ Y(N)=\Psi(N)\theta^0+V(N) y(n)y(n+1)...y(N)=y(n1)y(n)...y(N1)y(n2)y(n1)...y(N2)............y(0)y(1)...y(Nn)u(n1)u(n)...u(N1)u(n2)u(n1)...u(N2)............u(0)u(1)...u(Nn)a10a20...an0b10b20...bn0+v(n)v(n+1)...v(N)Y(N)=Ψ(N)θ0+V(N)

实际估计的参数中不仅包含噪声项,还应当包含估计误差:
Y ( N ) = Ψ ( N ) θ + ε ( N , θ ) Y(N)=\Psi(N)\theta+\varepsilon(N,\theta) Y(N)=Ψ(N)θ+ε(N,θ)

为了使误差最小化,定义性能指标:
J ( θ ) = ( Y − Ψ θ ) T ( Y − Ψ θ ) J(\theta)=(Y-\Psi\theta)^T(Y-\Psi\theta) J(θ)=(YΨθ)T(YΨθ)

对性能指标求导得到:
θ ^ L S = ( Ψ T Ψ ) − 1 Ψ T Y \hat \theta_{LS}=(\Psi^T\Psi)^{-1}\Psi^TY θ^LS=(ΨTΨ)1ΨTY

其应用需要两个条件:

  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 ) + ζ ( k ) y(k)=a_1y(k-1)+a_2y(k-2)+b_1u(k-1)+b_2u(k-2)+\zeta(k) y(k)=a1y(k1)+a2y(k2)+b1u(k1)+b2u(k2)+ζ(k),其输入为幅值为1的伪随机二位式序列,噪声是一个方差可调的标准正态分布的随机序列,真实值 a 1 = 1.5 , a 2 = − 0.7 , b 1 = 1.0 , b 2 = 0.5 a_1=1.5,a_2=-0.7,b_1=1.0,b_2=0.5 a1=1.5,a2=0.7,b1=1.0,b2=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 a_1=1.4977,a_2=-0.6986,b_1=0.9937,b_2=0.5020 a1=1.4977,a2=0.6986,b1=0.9937,b2=0.5020

二、加权最小二乘法

1、推导

性能指标增加加权系数后,定义为:
J ( θ ) = ( Y − Ψ θ ) T W ( Y − Ψ θ ) J(\theta)=(Y-\Psi\theta)^TW(Y-\Psi\theta) J(θ)=(YΨθ)TW(YΨθ)

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

对性能指标求导得到:
θ ^ W L S = ( Ψ T W Ψ ) − 1 Ψ T W Y \hat \theta_{WLS}=(\Psi^TW\Psi)^{-1}\Psi^TWY θ^WLS=(ΨTWΨ)1ΨTWY

其应用同样需要:

  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 a_1=1.5045,a_2=-0.6969,b_1=1.0097,b_2=0.4989 a1=1.5045,a2=0.6969,b1=1.0097,b2=0.4989

三、递推最小二乘法

1、推导

当增加一对新的输入输出时,
P ( N + 1 ) = P ( N ) − P ( N ) ψ ( N + 1 ) ψ T ( N + 1 ) P ( N ) w − 1 ( N + 1 ) + ψ T ( N + 1 ) P ( N ) ψ ( N + 1 ) L ( N + 1 ) = P ( N + 1 ) ψ ( N + 1 ) w ( N + 1 ) Δ ( N + 1 ) = y ( N + 1 ) − ψ T ( N + 1 ) θ ^ W L S ( N ) θ ^ W L S ( N + 1 ) = θ ^ W L S ( N ) + L ( N + 1 ) Δ ( N + 1 ) 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) P(N+1)=P(N)w1(N+1)+ψT(N+1)P(N)ψ(N+1)P(N)ψ(N+1)ψT(N+1)P(N)L(N+1)=P(N+1)ψ(N+1)w(N+1)Δ(N+1)=y(N+1)ψT(N+1)θ^WLS(N)θ^WLS(N+1)=θ^WLS(N)+L(N+1)Δ(N+1)

一般令 θ ^ W L S ( 0 ) = 0 \hat \theta_{WLS}(0)=0 θ^WLS(0)=0 P ( 0 ) = ( 1000 , 2000 , . . . , 5000 ) u 0 I P(0)=(1000,2000,...,5000)u_0I P(0)=(1000,2000,...,5000)u0I u 0 u_0 u0为输入信号的峰值,辨识参数有几个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 a_1=1.5009,a_2=-0.7072,b_1=0.9997,b_2=0.4915 a1=1.5009,a2=0.7072,b1=0.9997,b2=0.4915

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

  • 0
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值