最小二乘法参数估计

1、 【2-1】设某物理量Y与X1、X2、X3的关系如下：Y=1X1+2X2+3X3由试验获得的数据如下表。试用最小二乘法确定模型参数1、2和3X1:0.62 0.4 0.42 0.82 0.66 0.72 0.38 0.52 0.45 0.69 0.55 0.36X2:12.0 14.2 14.6 12.1 10.8 8.20 13.0 10.5 8.80 17.0 14.2 12.8X3:5.20 6.10 0.32 8.30 5.10 7.90 4.20 8.00 3.90 5.50 3.80 6.20Y:51.6 49.9 48.5 50.6 49.7 48.8 42.6 45.9 37.

2、8 64.8 53.4 45.3解：MATLAB程序为：Clear all;A= 0.6200 12.000 5.20000.400014.2000 6.10000.4200 14.6000 0.32000.8200 12.10008.30000.6600 10.8000 5.10000.7200 8.20007.90000.3800 13.00004.20000.5200 10.5000 8.00000.45008.80003.90000.690017.00005.50000.550014.20003.80000.360012.80006.2000;B=51.6 49.9 48.5 50.6

3、 49.7 48.8 42.6 45.9 37.8 64.8 53.4 45.3'C=inv(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.59032.44660.4597 【2-3】考虑如下模型其

4、中w(t)为零均值、方差为1的白噪声。根据模型生成的输入/输出数据u(k)和y(k),分别采用批处理最小二乘法、具有遗忘因子的最小二乘法（=0.95）和递推最小二乘法估计模型参数（限定数据长度N为某一数值，如N=150或其它数值），并将结果加以比较。解：1、 批处理最小二乘法M文件如下：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列的零列向量，给输入赋初值0yk=z

5、eros(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

6、);%产生方波 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);endthetae=inv(phi'*phi)*phi'*y'%计算参数估计值thetaeJ=y*y'-2*(phi'*y')'*thetae+thetae'*phi'*phi*thetae%求最小估计误差平方和cl

7、ear 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列的零列向量，给输入赋初值0yk=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

8、 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 y

9、k(1)=y(k);endthetae=inv(phi'*phi)*phi'*y'%计算参数估计值thetaeJ=y*y'-2*(phi'*y')'*thetae+thetae'*phi'*phi*thetae%求最小估计误差平方和运行结果如下：thetae = -1.3022 0.2982 0.9709 0.5091J = 186.0817与真值比较参数a1a2b0b1真值-1.30.310.5估计-1.30220.29820.97090.50912、 递推最小二乘法clear all;close all;a=1 -1.

10、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列的零列向量，给输入赋初值0yk=zeros(na,1);%输出初值u=randn(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+nb+1,1);%thetae的chuzhiP=

11、106*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(

12、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*randn(N,1),50*randn(N,1),100*

13、randn(N,1)所得结果如下图1）u=25*randn(N,1)2) u=50*randn(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列的零列向量，给输入赋初值0yk=zeros(na,1);%输出初值u=1*randn

14、(N,1);%输入采用白噪声序列xi=randn(N,1);%产生均值为零方差为1的白噪声序列w(k)thetae_1=zeros(na+nb+1,1);%thetae的chuzhiP=106*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+

15、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);