最小二乘法参数估计
【2-1】设某物理量丫与XI 、X2、X3的关系如下：丫= 9 i X i + 9 2X 2+ 9 3X 3
由试验获得的数据如下表。 试用最小二 ?乘法确定模型参数9 1、 9 2 和9 3
X1: 0.62 0.4 0.42 0.82 0.66 0.72 0.38 0.52 0.45 0.69 0.55 0.36 X2: 12.0 14.2 14.6 12.1 10.8 8.20 13.0 10.5 8.80 17.0 14.2 12.8 X3: 5.20 6.10 0.32 8.30 5.10 7.90 4.20 8.00 3.90 5.50
3.80 6.20 Y: 51.6
49.9 48.5 50.6 49.7 48.8 42.6 45.9 37.8 64.8 53.4 45.3 解：MATLAB 程序为：
Clear all;
A= [0.6200 12.000 5.2000
0.4000 14.2000 6.1000
0.4200 14.6000 0.3200
0.8200 12.1000 8.3000
0.6600 10.8000 5.1000
0.7200 8.2000 7.9000
0.3800 13.0000 4.2000
0.5200 10.5000 8.0000
0.4500 8.8000 3.9000
0.6900 17.0000 5.5000
0.5500 14.2000 3.8000
0.3600 12.8000 6.2000
];
B=[51.6 49.9 48.5 50.6 49.7 48.8 42.6 45.9 37.8 64.8 53.4 45.3]'; C=i nv
(A'*A)*A'*B
=[0.62 12 5.2;0.4 14.2 6.1;0.42 14.6 0.32;0.82 12.1 8.3;
0.66 10.8 5.1;0.72 8.2 7.9;0.38 13 4.2;0.52 10.5 8;
0.45 8.8 3.9;0.69 17 5.5;0.55 14.2 3.8;0.36 12.8 6.2] 公式中的A 是
①N, B 是YN 运行M 文件可得结果：
在matlab 中的运行结果：
C=
29.5903 2.4466 0.4597
【2-3】 考虑如下模型
1 2 z 0.5z
1 2U (t) 1 1.3z 0.3z y(k),分别采用批处理最小二乘法、具有遗忘因子的最小二乘法(入 =0.95)和递
推最小二乘法估计模型参数(限定数据长度
N 为某一数值，如N=150或其它数 值），并将结果加以比较 解：
1、 批处理最小二乘法
M 文件如下：
y(t)
w(t) 其中w(t)为零均值、方差为 1的白噪声。根据模型生成的输入 /输出数据u(k)和
a=[1 -1.3 0.3]';b=[1 0.5]';d=1;% 对象参数
na=length(a)-1;nb=length(b)-1;%na,nb为 A,B 阶次
N=200;%观测数据组
uk=zeros(d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0
yk=zeros( na,1);%输出初值
x1=1;x2=1;x3=1;x4=0;%移位寄存器初值
S=1;%方波初值
W=randn(N,1);%产生均值为零方差为1的白噪声序列w(k) theta=[a(2:
na+1);b];%theta 为列向量，[-1.3 0.3 1 0.5]'，为对象参数真值for k=1:N phi(k,:)=[-yk;uk(d:d+nb)]';%phi 为行向量，组成phi 矩阵
y(k)=phi(k,:)*theta+W(k);% 采集输出数据
IM=xor(S,x4);%进行异或运算，产生逆M序列
if IM==0
u(k)=-1;
else
u(k)=1;
end
S=not(S);%产生方波
M=xor(x3,x4);% 进行异或运算，产生M 序列x4=x3;x3=x2;x2=x1;x1=M;% 寄存器移位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);
end
thetae=inv(phi'*phi)*phi'*y'% 计算参数估计值thetae
J=y*y'-2*(phi'*y')'*thetae+thetae'*phi'*phi*thetae% 求最小估计误差平方和
clear all;close all;
a=[1 -1.3 0.3]';b=[1 0.5]';d=1;% 对象参数
na=length(a)-1;nb=length(b)-1;%na,nb为A,B 阶次
N=200;%观测数据组
uk=zeros(d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0
yk=zeros( na,1);%输出初值
x1=1;x2=1;x3=1;x4=0;% 移位寄存器初值
S=1;%方波初值
W=randn(N, 1 ) ; %产生均值为零方差为1 的白噪声序列w(k) theta=[a(2:
na+1);b];%theta 为列向量，[-1.3 0.3 1 0.5]'，为对象参数真值for k=1:N phi(k,:)=[-yk;uk(d:d+nb)]';%phi 为行向量，组成phi 矩阵
y(k)=phi(k,:)*theta+W(k);% 采集输出数据
IM=xor(S,x4);%进行异或运算，产生逆M序列if IM==0
u(k)=-1;
else u(k)=1;
end
S=not(S);%产生方波
M=xor(x3,x4);% 进行异或运算，产生M 序列x4=x3;x3=x2;x2=x1;x1=M;% 寄存器移位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);
end thetae=inv(phi'*phi)*phi'*y'% 计算参数估计值thetae J=y*y'-
2*(phi'*y')'*thetae+thetae'*phi'*phi*thetae% 求最小估计误差平方和运行结果如下：thetae =
-1.3022
0.2982
0.9709
0.5091
J =
186.0817 与
真值比较
参数a1 a2 b0 b1
真值-1.3 0.3
1 0.5
估计-1.3022 0.2982 0.9709 0.5091
2、递推最小二乘法
clear all;close all;
a=[1 -1.3 0.3]';b=[1 0.5]';d=1;% 对象参数
na=length(a)-1;nb=length(b)-1;%na,nb 为A,B 阶次
N=200; %观测数据组
uk=zeros(d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0 yk=zeros (n a,1);% 输出初值
u=ra ndn( N,1);%输入采用白噪声序列
w=randn(N,1);%产生均值为零方差为1的白噪声序列w(k)
theta=[a(2:na+1);b];%theta为列向量，[-1.3 0.3 1 0.5]',为对象参数真值thetae_1=zeros (na+n b+1,1);%thetae的chuzhi
P=10A6*eye( na+nb+1);%P 初值
for k=1:N
phi=[-yk;uk(d:d+nb)];%phi 为行向量，组成phi矩阵
y(k)=phi'*theta+w(k);% 采集输出数据
K=P*phi/(1+phi'*P*phi); thetae(:,k)=thetae_1+K*(y(k)-
phi'*thetae_1);
P=(eye(na+nb+1)-K*phi')*P; 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);
end
%画参数随采样时刻的估计结果
plot([1:N],thetae);
xlabel('k');ylabel(' 参数估计a,b');
legend('a1','a2','b0','b1');
axis([0 N -2 2]);
运行结果：
图1 递推最小二乘法参数估计
图中可以看出，辨识过程很不稳定，参数波动较大，因为在试验中将输入和噪声幅值一样大, 噪声干扰相对输入较大，所以应该调整输入和噪声的幅值
分别另：
u=25*ra ndn (N,1),50*ra ndn( N,1),100*ra ndn (N,1)所得结果如下图
1) u=25*randn(N,1)
20 dO GO80100 150 1A0 160 1S0 200
k
2)u=50*ra ndn
(N,1)
3)u=100*randn (N,1)
综上可以看出只有当输入幅值较大时系统辨识的效果会更好，参数波动也不是很大。
3、遗忘因子最小二乘法
clear all;close all;
a=[1 -1.3 0.3]：b=[1 0.5]';d=1;% 对象参数
na=length(a)-1;nb=length(b)-1;%na,nb 为A,B 阶次
N=200; %观测数据组
uk=zeros(d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0
yk=zeros (n a,1);% 输出初值
u=1*randn(N,1);% 输入采用白噪声序列
xi=randn(N,1);% 产生均值为零方差为1的白噪声序列w(k)
thetae_1=zeros (na+n b+1,1);%thetae的chuzhi
P=10A6*eye( na+nb+1);%P 初值
lambda=0.95%遗忘因子
for k=1:N theta(:,k)=[a(2:na+1);b];% 对象参数真值phi=[-yk;uk(d:d+nb)];%phi 为行向量，组成phi矩阵y(k)=phi'*theta(:,k)+xi(k);% 采集输出数据
K=P*phi/(lambda+phi'*P*phi); thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);
P=(eye(na+nb+1)-K*phi')*P/lambda;
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);
end
%画参数随采样时刻的估计结果
subplot(1,2,1);
plot([1:N],thetae(1:na,:));hold on;plot([1:N],theta(1:na,:),'k:');
xlabel('k');ylabel(' 参数估计a');
legend('a1','a2');
axis([0 N -2 2]);
subplot(1,2,2)
plot([1:N],thetae(na+1:na+nb+1,:));hold on;plot([1:N],theta(na+1:na+nb+1,:),'k:'); xlabel('k');ylabel(' 参数估计b');
legend('b0','b1');
axis([0 N -0.5 2]);
运行结果：