《最小二乘法参数估计.doc》由会员分享，可在线阅读，更多相关《最小二乘法参数估计.doc（12页珍藏版）》请在三一文库上搜索。
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.8
2、 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】考虑如下模型其中w(t)为零均值、方差为1的
4、白噪声。根据模型生成的输入/输出数据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=zeros(na,1);%输出初值x1=1;x2
5、=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
6、序列 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%求最小估计误差平方和clear all;close all;a=1 -1.3 0.3;b=1 0.5;d=1;%对象参数na=length(a)-1;nb=length(b)-1
7、;%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 k=1:N phi(k,:)=-yk;uk(d:d+nb);%phi为行向量，组成phi矩阵 y(k)=phi(k,:)*theta+W(k);%采集输出数据 IM=xor(S,
8、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);endthetae=inv(phi*phi)*phi*y%计算参数估计值thetaeJ=y*y-2*(phi*y)*thetae+thetae*phi*phi*thet
9、ae%求最小估计误差平方和运行结果如下：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.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=
10、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=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*t
11、hetae_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 递推最小二乘法参数估计图中可以看出，辨识过程很不稳定，参数波动较大，因为在试验中将输入和噪声幅
12、值一样大，噪声干扰相对输入较大，所以应该调整输入和噪声的幅值分别另：u=25*randn(N,1),50*randn(N,1),100*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(
13、d+nb,1);%创建d+nb行1列的零列向量，给输入赋初值0yk=zeros(na,1);%输出初值u=1*randn(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);%采集输出数据
14、 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:
15、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);运行结果：改变输入信号的幅值，使u=12.5*randn(N,1),25*randn(N,1),50*randn(N,1),100*randn(N,1)1) u=12.5*randn(N,1)2) u=25*randn(N,1)3) u=50*randn(N,1)4) 4)u=100*randn(N,1)比较可知：只有当输入幅度较大时系统辨识结果更好【3-2】设有被控过程：给定期望传递函数的分母多项式为，试按照极点配置方法设计控制系统，使期望输出无稳态误差，并写出控制表达式u(k)。