最小二乘法参数估计.docx
- 文档编号:15956808
- 上传时间:2023-07-09
- 格式:DOCX
- 页数:14
- 大小:125.39KB
最小二乘法参数估计.docx
《最小二乘法参数估计.docx》由会员分享,可在线阅读,更多相关《最小二乘法参数估计.docx(14页珍藏版)》请在冰点文库上搜索。
最小二乘法参数估计
【2-1】设某物理量Y与X1、X2、X3的关系如下:
Y=θ1X1+θ2X2+θ3X3
由试验获得的数据如下表。
试用最小二乘法确定模型参数θ1、θ2和θ3
X1:
0.620.40.420.820.660.720.380.520.450.690.550.36
X2:
12.014.214.612.110.88.2013.010.58.8017.014.212.8
X3:
5.206.100.328.305.107.904.208.003.905.503.806.20
Y:
51.649.948.550.649.748.842.645.937.864.853.445.3
解:
MATLAB程序为:
Clearall;
A=[0.620012.0005.2000
0.400014.20006.1000
0.420014.60000.3200
0.820012.10008.3000
0.660010.80005.1000
0.72008.20007.9000
0.380013.00004.2000
0.520010.50008.0000
0.45008.80003.9000
0.690017.00005.5000
0.550014.20003.8000
0.360012.80006.2000
];
B=[51.649.948.550.649.748.842.645.937.864.853.445.3]';
C=inv(A'*A)*A'*B
=[0.62125.2;0.414.26.1;0.4214.60.32;0.8212.18.3;
0.6610.85.1;0.728.27.9;0.38134.2;0.5210.58;
0.458.83.9;0.69175.5;0.5514.23.8;0.3612.86.2]
公式中的A是ΦN,B是YN,运行M文件可得结果:
在matlab中的运行结果:
C=
29.5903
2.4466
0.4597
【2-3】考虑如下模型
其中w(t)为零均值、方差为1的白噪声。
根据模型生成的输入/输出数据u(k)和y(k),分别采用批处理最小二乘法、具有遗忘因子的最小二乘法(λ=0.95)和递推最小二乘法估计模型参数(限定数据长度N为某一数值,如N=150或其它数值),并将结果加以比较。
解:
1、批处理最小二乘法
M文件如下:
clearall;closeall;
a=[1-1.30.3]';b=[10.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.30.310.5]',为对象参数真值
fork=1:
N
phi(k,:
)=[-yk;uk(d:
d+nb)]';%phi为行向量,组成phi矩阵
y(k)=phi(k,:
)*theta+W(k);%采集输出数据
IM=xor(S,x4);%进行异或运算,产生逆M序列
ifIM==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;%寄存器移位
fori=d+nb:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=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%求最小估计误差平方和clearall;closeall;
a=[1-1.30.3]';b=[10.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.30.310.5]',为对象参数真值
fork=1:
N
phi(k,:
)=[-yk;uk(d:
d+nb)]';%phi为行向量,组成phi矩阵
y(k)=phi(k,:
)*theta+W(k);%采集输出数据
IM=xor(S,x4);%进行异或运算,产生逆M序列
ifIM==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;%寄存器移位
fori=d+nb:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=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
与真值比较
参数a1a2b0b1
真值-1.30.310.5
估计-1.30220.29820.97090.5091
2、递推最小二乘法
clearall;closeall;
a=[1-1.30.3]';b=[10.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);%输出初值
u=randn(N,1);%输入采用白噪声序列
w=randn(N,1);%产生均值为零方差为1的白噪声序列w(k)
theta=[a(2:
na+1);b];%theta为列向量,[-1.30.310.5]',为对象参数真值
thetae_1=zeros(na+nb+1,1);%thetae的chuzhi
P=10^6*eye(na+nb+1);%P初值
fork=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);
fori=d+nb:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=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([0N-22]);
运行结果:
图1递推最小二乘法参数估计
图中可以看出,辨识过程很不稳定,参数波动较大,因为在试验中将输入和噪声幅值一样大,噪声干扰相对输入较大,所以应该调整输入和噪声的幅值
分别另:
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、遗忘因子最小二乘法
clearall;closeall;
a=[1-1.30.3]';b=[10.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);%输出初值
u=1*randn(N,1);%输入采用白噪声序列
xi=randn(N,1);%产生均值为零方差为1的白噪声序列w(k)
thetae_1=zeros(na+nb+1,1);%thetae的chuzhi
P=10^6*eye(na+nb+1);%P初值
lambda=0.95%遗忘因子
fork=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);
fori=d+nb:
-1:
2
uk(i)=uk(i-1);
end
uk
(1)=u(k);
fori=na:
-1:
2
yk(i)=yk(i-1);
end
yk
(1)=y(k);
end
%画参数随采样时刻的估计结果
subplot(1,2,1);
plot([1:
N],thetae(1:
na,:
));holdon;plot([1:
N],theta(1:
na,:
),'k:
');
xlabel('k');ylabel('参数估计a');
legend('a1','a2');
axis([0N-22]);
subplot(1,2,2)
plot([1:
N],thetae(na+1:
na+nb+1,:
));holdon;plot([1:
N],theta(na+1:
na+nb+1,:
),'k:
');
xlabel('k');ylabel('参数估计b');
legend('b0','b1');
axis([0N-0.52]);
运行结果:
改变输入信号的幅值,使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)。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 最小二乘法 参数估计