自动化2011级系统辨识大作业-王万秋

系统辨识大作业班级：自动化1101班姓名：王万秋学号：11052204第一题模仿index3，搭建如下的单输入－单输出系统的差分方程)()2()1()2()1()(2121kVkubkubkzakzakz)2()1()()(321kvckvckvckV取真值11.4a、20.8a、0.11b、20.5b、10.6c、2.12c和3.03c，输入信号采用4阶M序列，幅值为1。当)(kv的均值为0，方差分别为0.1和0.5的高斯噪声时，分别用一般最小二乘法、递推最小二乘法和增广递推最小二乘法估计参数。并通过对三种方法的辨识结果的分析和比较，说明上述三种参数辨识方法的优缺点。（15分）利用simulink搭建的模型框图如下：一般最小二乘法程序：u=UY(1:450,1)';%输入矩阵z=UY(1:450,2)';%输出矩阵H=zeros(400,4);fori=1:400H(i,1)=-z(i+1);H(i,2)=-z(i);H(i,3)=u(i+1);H(i,4)=u(i);endtheta=inv(H'*H)*H'*(z(3:402))'递推最小二乘法程序代码：u=UY(1:450,1)';%输入矩阵z=UY(1:450,2)';%输出矩阵P=100*eye(4);%估计方差Pstore=zeros(4,401);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];Theta=zeros(4,401);%参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];K=[10;10;10;10];fori=3:402h=[-z(i-1);-z(i-2);u(i-1);u(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];endi=1:401;theta=Theta(:,length(Theta(1,:)))subplot(2,1,1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:))title('递推最小二乘的估计值过度情况')legend('a1','a2','b1','b2')subplot(2,1,2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:))title('递推最小二乘估计方差过度情况')legend('a1','a2','b1','b2')增广递推最小二乘法程序代码：u=UY(1:450,1)';%输入矩阵z=UY(1:450,2)';%输出矩阵v=UY(1:450,3)';%噪声矩阵P=100*eye(6);%估计方差Pstore=zeros(6,300);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)];Theta=zeros(6,300);%参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3;3;3];K=[10;10;10;10;10;10];fori=3:300h=[-z(i-1);-z(i-2);u(i-1);u(i-2);v(i-1);v(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(6)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5),P(6,6)];endi=1:300;theta=Theta(:,length(Theta(1,:))-10)subplot(2,1,1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:),i,Theta(6,:))title('增广递推最小二乘估计值过渡情况')legend('a1','a2','b1','b2')subplot(2,1,2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5,:),i,Pstore(6,:))title('增广递推最小二乘估计方差过度情况')legend('a1','a2','b1','b2')辨识结果：噪声辨识方法一般最小二乘法递推最小二乘法增广递推最小二乘法噪声方差为0.1theta=1.11410.56360.98130.2415theta=1.11450.56400.98130.2420theta=1.40210.52540.98530.2438噪声方差为0.5theta=0.88100.41270.95870.0558theta=0.88120.41280.95870.0561theta=1.41820.56240.97110.2647噪声方差为0.1噪声方差为0.5由图表可得：一般最小二乘法和递推最小二乘法的辨识结果很接近，而增广递推最小二乘法的辨识结果和其他两种方法差距较大，尤其a1和b2。此外，噪声越大，对前两种方法辨识效果影响较大，对于增广递推最小二乘法，抗噪声能力较强。三种方法的优缺点：一般最小二乘法优点：白噪声可得无偏渐进无偏估计；算法简单、可靠；计算量最小；一次即可完成算法，适合于离线辨识。缺点：随着数据的增长，该方法将出现“数据饱和”现象。递推最小二乘法优点：可以减少数据存储量，避免矩阵求逆，兼骚运算量。缺点：为避免出现数据饱和现象，必须限制过去数据的影响。增广递推最小二乘法优点：比递推算法有着更高的准确度。增广最小二乘的递推算法对应的噪声模型为滑动平均噪声，扩充了参数向量和数据向量的维数，把噪声模型的辨识同时考虑进去。以上两种最小二乘法只能获得过程模型的参数估计，而增广最小二乘法同时又能获得噪声模型的参数估计。当数据长度较大时，辨识精度低于极大似然法。第二题模仿实验二，搭建对象，由相关分析法，获得脉冲响应序列ˆ()gk，由ˆ()gk，参照讲义，获得系统的脉冲传递函数()Gz和传递函数()Gs;应用最小二乘辨识,获得脉冲响应序列ˆ()gk;应用相关最小二乘法，拟合对象的差分方程模型；加阶跃扰动，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱)辨识二阶差分方程模型的参数，比较两种方法的辨识效果差异；采用不少于两种定阶方法，确定模型的阶次。（40分）搭建模型如下：1.相关分析法获得脉冲响应序列程序代码：Np=2^4-1;deta=2;a=5;U=zeros(Np,Np);fori=1:1:Npforj=1:1:NpU(i,j)=UY(Np+1-i+j,1);endendY=zeros(Np,1);fori=1:1:NpY(i,1)=UY(Np+i,2);endR=ones(Np,Np);fori=1:1:NpR(i,i)=2;endG=(1/(a*a*(Np+1)*deta))*R*U*Yt=0:2:28plot(t',G,'-*')获得的脉冲响应序列如下：G=0.01270.21050.30580.26900.18880.13300.07960.06030.05950.04190.02360.02280.01460.03370.0069响应序列的散点图：2.求解G(z)和G(s)求解G(z)程序代码：clcn=3;Ts=2;Np=15;u=UY(31:200,1);y=UY(31:200,2);f=zeros(120,6);f(:,1)=-1*y(31:150);f(:,2)=-1*y(30:149);f(:,3)=-1*y(29:148);f(:,4)=-1*u(31:150);f(:,5)=-1*u(30:149);f(:,6)=-1*u(29:148);yy=y(32:151);q=inv(f'*f)*f'*yy;a=q(1:3);b=q(4:6);a=a',b=b'G=tf(b,[1,a],2)则系统的离散传递函数：Transferfunction:-0.3946z^2-0.3493z-0.07706--------------------------------------z^3-0.7322z^2+0.02371z+0.01708Samplingtime:2因此，将离散传函转化为连续传函：Gs=d2c(G,'zoh')则有：Transferfunction:0.1374s^3-0.4647s^2-0.3989s-1.571----------------------------------------------s^4+3.063s^3+5.763s^2+3.892s+0.5904我觉得该连续传函不合理（阶次有点高）。3.应用最小二乘法获得脉冲响应序列Yy=UY(31:54,2)/2;fori=24:-1:1Uu(i,:)=UY(i+30:-1:16+i,1)';endGG=inv(Uu'*Uu)*Uu'*Yy;plot(0:2:29,GG,'b')holdonstem(0:2:29,GG,'b','filled')title('最小二乘法获得脉冲序列');获得的脉冲响应序列为：GG=0.05880.26170.37840.32430.24030.17860.14390.10480.08340.08110.07530.08530.07990.06830.0587响应序列的散点图：4.相关最小二乘法拟合对象的差分方程u=UY(1:500,1)';%输入矩阵z=UY(1:500,2)';%输出矩阵H=zeros(400,4);fori=1:400H(i,1)=-z(i+1);H(i,2)=-z(i);H(i,3)=u(i+1);H(i,4)=u(i);endTheta=inv(H'*H)*H'*(z(3:402))'计算结果：Theta=-0.96310.21450.38950.25925.加扰动后，用最小二乘法和带遗忘因子的最小二乘法辨识差分方程（并比较两种方法）一般最小二乘法辨识结果：theta=-1.21840.24410.44710.1591带遗忘因子的最小二乘法程序代码：na=2;nb=2;d=0;L=400;u=0.95;uk=UY(1:400,1)';yk=UY(1:400,2)';thetae_1=zeros(na+nb,1);P=10^6*eye(na+nb);fork=3:Lphi=[-yk((k-1):-1:(k-2))uk((k-1):-1:(k-2))];K=P*phi'/(u+phi*P*phi');thetae(:,k)=thetae_1+K*(yk(k)-phi*thetae_1);P=(eye(na+nb)-K*phi)*P/u;thetae_1=thetae(:,k);fori=d+nb:-1:2uk(i)=uk(i-1);endfori=na:-1:2yk(i)=yk(i-1);endendplot([1:L],thetae);theta=thetae(:,400)legend('a1','a2','b1','b2');title('带遗忘因子的最小二乘法便是参数')辨识结果：theta=-1.09070.10040.45980.2370比较两种辨识方法结果可得：a1和b1比较一致，而a2和b2差别较大，带遗忘因子辨识过程波动较大。6.两种定阶方法确定模型的阶次第一种方法——AIC定阶方法：程序代码：L=200;U=UY(:,1);Y=UY(:,2);forn=1:1:5N=100-n;yd(1:N,1)=Y(n+1:n+N);fori=1:Nforl