最优估计作业5.docx
《最优估计作业5.docx》由会员分享，可在线阅读，更多相关《最优估计作业5.docx（12页珍藏版）》请在冰豆网上搜索。
最优估计作业5
姓名朱凯歌学号**********
Homework1
（1）批处理：
仿真如图所示：
（2）递推算法：
取
的初值为
，
，得估计
仿真如图所示：
分析：
由仿真结果可得，观测数据长度增加，递推估计参数精度越高。
程序代码：
%Themovementparameterofaaircraftevolveswithtimeaccording
%toy（t）=a*t^2+b*t+c.Sampleineach1speriod,andthusobtain20
%datapair（1,y1）,...,（20,y20）.ObtaintheLSestimateofa,b
%andc.
%Dataofy:
2.98284.5256.11557.2348.43299.125910.180010.8600
%10.930011.141010.609010.480010.38309.58088.46117.46786.3942
%4.15923.00290.5503
%Requirement
%
（1）preparematlabprogramsforbatchandrecursivemethods
%
（2）submitthetechnicalreportincludingtheestimatesofa,bandc;
%plotthet-vs-estimatecurveofa,bandc,analyzetherelationship
%betweendatalengthandestimateaccuracy;
%plotmeasurementdataandactualparameters
clearall
closeall
%实验数据
Y=[2.98284.5256.11557.2348.43299.125910.180010.860010.9300...
11.141010.609010.480010.38309.58088.46117.46786.39424.1592...
3.00290.5503]';
H=[];
fori=1:
20
H=[H;i^2i1];
end
%批处理算法
E_abc_bat=inv（H'*H）*H'*Y;
Y_E_bat=H*E_abc_bat;
figure
（1）;
plot（1:
20,Y,'-g*',1:
20,Y_E_bat,'-bo'）;
legend（'量测数据曲线','估计曲线'）;
xlabel（'t'）
ylabel（'y'）
title（'最小二乘批处理算法'）
%递推算法
X0=zeros（3,1）;
P=1000000*diag（[111]）;
X_re=zeros（3,21）;
X_re（:
1）=X0;
fori=1:
20
P=P-P*H（i,:
）'*...
inv（H（i,:
）*P*H（i,:
）'+1）*H（i,:
）*P;
K=P*H（i,:
）';
X_re（:
i+1）=X_re（:
i）+K*（Y（i）-H（i,:
）*X_re（:
i））;
end
E_abc_re=X_re（:
21）;
Y_E_re=H*E_abc_re;
figure
（2）;
plot（1:
20,Y,'-g*',1:
20,Y_E_re,'-bo'）;
legend（'量测数据曲线','估计曲线'）;
xlabel（'t'）
ylabel（'y'）
title（'最小二乘递推算法'）
figure（3）;
plot（1:
20,X_re（1,2:
21）,'-g*',1:
20,X_re（2,2:
21）,'-bo',...
1:
20,X_re（3,2:
21）,'-rs'）;
legend（'a','b','c'）;
xlabel（'t'）
ylabel（'估计值'）
title（'数据长度与递推估计参数变化关系'）
gridon
Homework2
（1）30次独立量测
Z1=
Columns1through30
1.29954.24903.50423.13651.99583.0407
0.18964.4116-1.05332.90123.54836.8827
0.1460-1.3455-1.77763.1718-3.89301.1855
2.65012.51103.19032.43983.03683.0646
1.00532.72344.12135.54831.42691.3387
-1.39590.1632-2.87020.37640.71071.4107
3.28503.08092.18431.62052.12772.5052
-1.09664.8442-2.97516.97563.96611.4475
-0.45190.2434-0.4241-0.1912-0.87611.6778
2.77312.03443.24463.42544.06933.6951
2.85841.9805-2.0903-3.10110.7746-4.9242
0.50490.40840.0905-1.23770.03681.0777
2.51632.80710.99630.78024.01422.8124
0.84921.62673.01610.5196-4.72122.7393
0.5969-0.04091.67550.8411-1.0760-2.8060
（2）采用递推最小二乘估计，选取初值
，
，仿真如图：
分析：
随着时间的推移，X的估计趋于真实值且趋于稳定。
（3）100次仿真得到X的每一个分量元素估计样本方差曲线
分析：
X的每一个分量元素估计的抽样误差方差随着量测数据增多逐渐减少并趋于稳定值。
（4）这里，
被认为是理论方差，仿真如图所示：
分析：
X的每一个分量元素的理论方差随着量测数据增多逐渐减少并趋于稳定值。
程序代码：
%ConsidertheactualbutunknownvectorX=[320]',themeasurement
%matrixHisa3*3unitematrix,measurementnoiseiszero-mean
%Gaussianwithcovariancediag{1,9,2}.
%
（1）bringout30independentmeasurements
%
（2）plottheLSestimateofeachelementinXvsobservationtimes;
%Alsoplotmeasurements
%（3）plotsamplederrorvarianceofeachelementinXbasedon100simulations
%（4）plotthecorrespondingtheoreticalvariance
clearall
closeall
X=[320]';
COV=diag（[1,9,2]）;
H=eye（3）;
%30次独立观测样本
Z1=zeros（3,30）;
fori=1:
30
Z1（:
i）=H*X+sqrt（COV）*randn（3,1）;
end
%递推最小二乘估计
X0=zeros（3,1）;
P=1000000*diag（[111]）;
X_EST=zeros（3,31）;
X_EST（:
1）=X0;
fori=1:
30
P=P-P*H'*inv（H*P*H'+eye（3））*H*P;
K=P*H';
X_EST（:
i+1）=X_EST（:
i）+K*（Z1（:
i）-H*X_EST（:
i））;
end
X_EST=X_EST（:
2:
31）;
figure
（1）
plot（1:
30,X_EST（1,:
）,'*g-',1:
30,Z1（1,:
）,'-bo'）;
legend（'参数1估计曲线','量测数据曲线'）;
xlabel（'k'）
ylabel（'X1'）
title（'X1的递推最小二乘估计'）
gridon
figure
（2）
plot（1:
30,X_EST（2,:
）,'*g-',1:
30,Z1（2,:
）,'-bo'）;
legend（'参数2估计曲线','量测数据曲线'）;
xlabel（'k'）
ylabel（'X2'）
title（'X2的递推最小二乘估计'）
gridon
figure（3）
plot（1:
30,X_EST（3,:
）,'*g-',1:
30,Z1（3,:
）,'-bo'）;
legend（'参数3估计曲线','量测数据曲线'）;
xlabel（'k'）
ylabel（'X3'）
title（'X3的递推最小二乘估计'）
gridon
%100次仿真
forN=1:
100
Z=zeros（3,30）;
fori=1:
30
Z（:
i）=H*X+sqrt（COV）*randn（3,1）;
end
X0=zeros（3,1）;
P=1000000*diag（[111]）;
X_EST=zeros（3,31）;
X_EST（:
1）=X0;
fori=1:
30
P=P-P*H'*inv（H*P*H'+eye（3））*H*P;
K=P*H';
X_EST（:
i+1）=X_EST（:
i）+K*（Z（:
i）-H*X_EST（:
i））;
P_EST（:
:
i）=P;
end
X_EST=X_EST（:
2:
31）;
X_EST_N（:
:
N）=X_EST;
end
%计算样本方差
fori=1:
3
forj=1:
30
v=var（X_EST_N（i,j,:
））;
V（i,j）=v;
end
end
figure（4）
subplot221;
plot（1:
30,V（1,:
）,'-r*'）;
legend（'参数1的样本方差曲线'）;
xlabel（'k'）;
ylabel（'方差'）;
title（'100次仿真下参数1的样本方差'）;
subplot222;
plot（1:
30,V（2,:
）,'-go'）;
legend（'参数2的样本方差曲线'）;
xlabel（'k'）;
ylabel（'方差'）;
title（'100次仿真下参数2的样本方差'）;
subplot223;
plot（1:
30,V（3,:
）,'-bs'）;
legend（'参数3的样本方差曲线'）;
xlabel（'k'）;
ylabel（'方差'）;
title（'100次仿真下参数3的样本方差'）;
%参数的理论方差
V_T（1,:
）=P_EST（1,1,:
）;
V_T（2,:
）=P_EST（2,2,:
）;
V_T（3,:
）=P_EST（3,3,:
）;
figure（5）
subplot221;
plot（1:
30,V_T（1,:
）,'-r*'）;
legend（'参数1的理论方差曲线'）;
xlabel（'k'）;
ylabel（'方差'）;
title（'参数1的理论方差'）;
subplot222;
plot（1:
30,V_T（2,:
）,'-go'）;
legend（'参数2的理论方差曲线'）;
xlabel（'k'）;
ylabel（'方差'）;
title（'参数2的理论方差'）;
subplot223;
plot（1:
30,V_T（1,:
）,'-bs'）;
legend（'参数3的理论方差曲线'）;
xlabel（'k'）;
ylabel（'方差'）;
title（'参数3的理论方差'）;
升级会员 