系统辨识课程报告.docx

《系统辨识课程报告.docx》由会员分享，可在线阅读，更多相关《系统辨识课程报告.docx（24页珍藏版）》请在优知文库上搜索。

1、雷/7n/挈系统辨识课程报告学号:姓名:西北工业大学研究生院一.设SISO系统差分方程为y(k)=-aly(k-1)-a2y(k-2)+bxuk-1)+b2uk-2)+(k)辨识参数向量为=qa2ab2f输入输出数据详见数据uyl.txtuy3.txt。察Q为噪声方差各异的白噪声或有色噪声。试求解：D用n元一次方程解析法，再求其平均值方法估计。分析：方程解析法要求n元一次方程对应n个方程，由于噪声未知；辨识向量。=aa2Ub2,故选择每四组输入输出数据为一组方程组求解，最后求取均值；考虑到本文给出的系统为2阶系统，构建理论上输出y仅与前两时刻的y与U相关，在不考虑噪声情况下，设置n元一次方程的

2、系数矩阵与结果矩阵为：a=-y(2:499);b=-y(l:498);c=u(2:499);d=u(1:498);X=y(3:500);Hfori=1:495A=a(i),b(i),c(i),d(i).a(i+l),b(i+l),c(i+l),d(i+l).a(i+2),b(i+2),c(i+2),d(i+2);.a(i+3),b(i+3),c(i+3),d(i+3);B=x(i);x(i+l);x(i+2);x(i+3);thetaQ(:,i)=AB;end-theta=mean(theta,2);算法将依次选取a,b,c,Cl中连续四点数据构建系数矩阵A,同时选取对应四点输出构建结果矩阵B

3、,则对应每次的6求取结果为：thetaO(:,i)=AB;最终对theta按行求取均值得到theta估计值:123411.96661.6130-3.352221.21240.2895-2.949230.5702-0.2678-2.789840.28390.69470.3957按列依次为uyl,uy2,uy3数据。2）用最小二乘及递推最小二乘法估计最小二乘原理:构建参数矩阵y(n)7W(M+1)(1)-(w+1)-y(2)u(n+2)(2)y(IN1)-V(N)(力+N)u(N)-则有最小二乘估计为=(l)-,ij算法实现：functiontheta=LS(u,y)phi=-y(2:end-l)

4、,-y(1:end-2),u(2:end-l),u(1:end-2);|theta=(phi,*phi)phi,*y(3:end);结果：LStheta,4x3double12341.48551.11131.11580.78690.49630.48010.48370.37910.42540.19820.18790.1245递推最小二乘：构建PO与thetaO：P、0-(yf0o)，5、O=PNO0YM)则每次数据更新后的递推算法为：。74I=XJ+KN4(X、.I-V,vI)K十I=P+(l+V,*IP+1)1P.=P-PWy.(1+Wt+1PW、.|)”马，P算法实现：计算初始P,theta

5、(选取前50点计算)phi=-y(2:49),-y(1:48),u(2:49),u(l:48);P=inv(phi,*phi);theta=P*phi*y(3:50);数据更新计算每次theta：fori=50：length(u)-2phis=-y(i+l),-y(i),u(i+l),u(i)l,;K=P*phis/(1+phis,*P*phis);P=P-K*phis,*P;theta=theta+K*(y(i+2)-phis,*theta);theta=theta,;end结果：12341.48471.11231.11680.78650.49630.4815|0.48420.37970.4

6、2610.19780.18780.12483)用辅助变量法及其递推算法估计e；辅助变量：叫I构建辅助变量矩阵Z：.v()一、(1)+11(1)y(ni1)-y(2)(+2)(2) -wv(n+N1)一&(N)+N)“(N)则有OIV=(ZT0)7ZTy由于Z无法先验得知，故需先通过最小二乘估计theta,计算得到Z再通过辅助变量估计theta,循环迭代至收敛算法设计：输入上步通过最小二乘所得theta,计算先验Z阵，迭代求解至theta收敛functionIVtheta=IV(u,y,theta)phi=-y(2:end-l),-y(1:end-2),u(2:end-l),u(1:end-2)

7、;IVtheta=theta;while1Z=phi;yl=Z*IVtheta;|Z(2:length(u)-2,1)=-yl(1:length(u)-3,1);Z(3:length(u)-2,2)=-yl(1:length(u)-4,1);IVtheta=(Z,*phi)Z,*y(3:end);ifabs(IVtheta-theta)e计算U,y修正值.Q)=y(6)+7V)y(1)+2(-m)tttn()=w()+lu(-1)+fiu(km)重新使用最小二乘法估计theta=(l)lry循环迭代至theta收敛算法设计：计算残差：e(l)=0;e(2)=0;e(3:long)=y(3:en

8、d)-phi*theta;定义OnIU阵与f函数，更新u,y重新估计thetaoum=-e(2:end-l),-e(l:end-2),;f=(oum,*oum)oum,*e(3:end);Y(3:long)=y(3:end)+f(l)*y(2:end-1)+f(2)*y(l:end-2);U(3:long)=u(3:end)+f(l)*u(2:end-1)+f(2)*u(l:end-2);phi=-Y(2:end-1),-Y(1:end-2),U(2:end-1),U(1:end-2);theta=(phi,*phi)phi,*Y(3:end);循环迭代至收敛：ifabs(theta-thet

9、a)l.v)Pa=P华-K知成TP伊KaI=P*%N+l(1+t.Pg)N.I)1P0)=(CSn)T算法设计：求Ptheta初值Y=y；U=u;phil=-y(2:49),-y(l:48),u(2:49),u(l:48);P=inv(phi*phil);P2=P;theta=P*phi*y(3:50);循环计算新theta值fori=50：length(u)-2phi=-y(2:i+l),-y(l:i),u(2:i+l),u(l:i)J;phis=-y(i+l),-y(i),u(i+l),u(i)l,;K=P*phis/(1+phis,*P*phis);P=P-K*phis,*P;theta

10、=theta+K*(y(i+2)-phis,*theta);e(l)=0;e(2)=0;e(3:i+2)=y(3:i+2)-phi*theta;oum=-e(2:i+1),-e(l:i),;f=(oum,*oum)oum,*e(3:i+2),;Y(3:1+2)=y(3:i+2)+f(l)*y(2:i+l)+f(2)*y(l:i);U(3:i+2)=u(3:i+2)+f(l)*u(2:i+l)+f(2)*u(l:i);phi2=-Y(i+l),-Y(i),U(i+l),U(i);K2=P2*phi2/(1+phi2,*P*phi2);P2=P2-K2*phi2,*P2;theta=theta+K2*(y(i+2)-phi2,*theta);求解结果为：1231.48511.06701.14300.78200.46250.53710.49210.37600.44070.19560.19380.16745）用夏氏偏差修正法、夏氏改良法及其递推算法估计出夏氏偏差修正法:定义M阵：