一、最小二乘法
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(k−1)+a20y(k−2)+...+an0y(k−n)+b10u(k−1)+b20u(k−2)+...+bn0u(k−n)+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(k−1)y(k−2)...y(k−n)u(k−1)u(k−2)...u(k−n)]⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡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(n−1)y(n)...y(N−1)y(n−2)y(n−1)...y(N−2)............y(0)y(1)...y(N−n)u(n−1)u(n)...u(N−1)u(n−2)u(n−1)...u(N−2)............u(0)u(1)...u(N−n)⎦⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡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
其应用需要两个条件:
- N>3n-1。n是系统的阶次,N是数据的长度。
- 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(k−1)+a2y(k−2)+b1u(k−1)+b2u(k−2)+ζ(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)=arN−k。
对性能指标求导得到:
θ
^
W
L
S
=
(
Ψ
T
W
Ψ
)
−
1
Ψ
T
W
Y
\hat \theta_{WLS}=(\Psi^TW\Psi)^{-1}\Psi^TWY
θ^WLS=(ΨTWΨ)−1ΨTWY
其应用同样需要:
- N>3n-1。n是系统的阶次,N是数据的长度。
- 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)−w−1(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。
本文是作者在日常学习生活中所作,难免有遗漏或错误,遇到问题的读者请点击给我写信向我的邮箱反馈,不胜感激。