系统辨识课程报告.docx

雷/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元一次方程的系数矩阵与结果矩阵为：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,则对应每次的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),-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=PNO>0YM)则每次数据更新后的递推算法为：。74I=XJ+KN4(X、.I-V,v>I)K十I=P+(l+V,*IP+1)1P.=P-PWy.(1+Wt+1PW、.|)”马，P算法实现：计算初始P,theta(选取前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.42610.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);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)<0.00001break;endtheta=IVtheta;end结果：击4x3double12311.53531.39651.617220.83320.70720.868730.48730.38020.432540.21560.29700.3138递推辅助变量法:同RLS,取PN为Pv=(Z;v),VI-V÷K71(yV+1-QI)P+1=P.-K+1(1PKI=PZ、1I(1+WI.IPZJI)1算法设计:functionRIVtheta=RIV(u,y,RIVtheta)phi=-y(2:49),-y(l:48),u(2:49),u(l:48);Z=phi;yl=Z*RIVtheta;Z(2:48,1)=-yl(l:47,1);Z(3:48,2)=-yl(l:46,1);RIVtheta=(Z'*phi)Z'*y(3:50);P=inv(Z,*phi);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;RIVtheta=RIVtheta+K*(y(i+2)-phis,*RIVtheta);-end结果:1231.48611.13711.16720.78820.52830.52380.48430.37780.42640.19840.20190.14494）用广义最小二乘法及其递推算法估计。；广义最小二乘法：利用最小二乘估计计算残差：e（Q=G（之7）、（笈）-5（十九e"）（4）=y（八）+d"y（A-1）+-7?）-瓦1）“（）BTLa几），A=r,n+1,，"+N用估计残差代替模型噪声影响，计算。mu阵与f函数:-C(-1)-C(一加+1)一+e()一2(一加+2).-e(,)(+/V-1)-e(+N-2)-e(l)(n÷N-m)J求得f函数估值A)=(。g>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:end)-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-theta)<0.OOOOlbreak;end结果：递推广义最小二乘法：,v+I=ON+KeL(y、r1-WtPMl=PS)-KSI=p86n“(i+吗+PG)WN+)'p=(iv)-1fN+I/n+I(°n+1-3工>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=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阵：