系统辨识课程报告.docx

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

1、前肛7去上当系统辨识课程报告学 号:姓 名:西北工业大学研究生院一.设SISO系统差分方程为y (k) = - a1 y(k -1) - a2 y(k -2) + b1u(k -1) + b2u(k -2) + g(k)辨识参数向量为 =a1 a2 b1 b2,输入输出数据详见数据uyl. txtuy3. txt。久左)为噪声方差各异的白噪声或有色噪声。试求解：D用n元一次方程解析法，再求其平均值方法估计夕分析：方程解析法要求n元一次方程对应n个方程，由于噪声未知；辨识向 量夕=U&b1 b2V,故选择每四组输入输出数据为一组方程组求解，最后 求取均值；考虑到本文给出的系统为2阶系统，构建理论

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);0 for i = 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);thelja(:, i)

3、= AB;end-theta = mean(theta, 2);算法将依次选取a, b, c, d中连续四点数据构建系数矩阵A,同时选取对应四点 输出构建结果矩阵B,则对应每次的夕 求取结果为：thetaO(:,i) = AB;最终对theta按行求取均值得到theta估计值：-123411.96661.61303352221.21240.2895-2.949230.5702-0.2678-2.789840.28390.69470.3957按列依次为uyl, uy2, uy3数据。2)用最小二乘及递推最小二乘法估计6;最小二乘原理: 构建参数矩阵一() 一 y( + 1)y( n I N -

4、1)一3一八2)-y(N)(H + 1) (n 2)(1)(2)则有最小二乘估计为G = (1),ry算法实现： function theta = 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);结果:I LStheta_i二 4x3 double12341.48551.11131.11580.78690.49630.480,0.4837037910.42540.19820.1879工 1245递推最小二乘：构建PO与thetaO：P (

5、yo ) 1 0 VO = P、0YYO则每次数据更新后的递推算法为：07 * I = V K4 16、+ I - VzN+ I 、)Kv+ = Pn + i(1 + 档 7 PW、+ I)TPNrl = Pn -PWNr(I + W+PWNT)TW工 R算法实现：计算初始P, theta (选取前50点计算)phi = E-y(2:49),-y(l :48), u(2:49),u(l :48);P = inv (phi, * phi);theta =P* Phi*y(3:50);数据更新计算每次theta：Ifor i - 50 : length(u) - 2phis =-y(i+l),-y

6、(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结果：1231.48471.1160.78650.49630.48150.48420.37970.42610.19780.18780.12483)用辅助变量法及其递推算法估计e；辅助变量：构建辅助变量矩阵Z：一.v()-y(l) / (F/ 1)M (1)-y(+ 1) 一y(2) (十2)”(2) -.v( + N-

7、l) 一y(N)u( w N) U(N)-则有v =(Z) ,Zj由于Z无法先验得知，故需先通过最小二乘估计theta,计算得到Z再通过辅助变量估计theta,循环迭代至收敛算法设计：输入上步通过最小二乘所得theta,计算先验Z阵，迭代求解至theta收敛 function IVtheta = IV(u, y, theta)phi = -y(2:end-1), -y(1:end-2), u(2：end-l), u(1:end-2);IVtheta = theta;E3while 1Z = phi;yl = Z * IVthetaq Z (2: length (u) - 2,1)Z (3:le

8、ngth (u) - 2,2)IVtheta = (Z, * phi)=-yl(1:length(u)-3,1);=-yl (1: length(u) -4, 1); Z, * y (3:end);if abs(IVtheta theta) r( - 1) - w )-一?%(-r),A = t n + 1,，+ N用估计残差代替模型噪声影响，计算。IllU阵与f函数:-W (-1) -e (-加+1)+e(V) -C(-m+2). e( + N- 1)-6)( +N-2)-e(n + N - m)-求得f函数估值A) = (OT (O)计算u, y修正值、(A) = yU) + 和、(左 -

9、 1) + + 为)- m)()=(女)+ /)(4 - I) + 7%(6-m)重新使用最小二乘法估计theta = (,) 1y循环迭代至theta收敛算法设计：计算残差：e(l) = 0;e(2) = 0;e (3:long) = y (3:end)- phi*theta;定义OmU阵与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)

10、= u(3:end) + f(l) * u(2:end 1)+ f(2) * u(l:end -2); phi = -Y(2:endl), Y(1:end-2), U(2：end-l), U(l: end-2);theta = (phi, * phi) phi, *Y (3: end);循环迭代至收敛：if abs(theta - theta)0.OOOOl break;end结果：递推广义最小二乘法：1 = On +- 飞/ 份、)PMl = PS)-GLK.必)Ki = P(lrP)1p9)=(哂八，)7n+ I = /n + KGlI (e.7 + I - 3工 + |八)P%= Pg-

11、KHMhPM)KHl = P9*N+l(l + +iP*Ln.)TP = (C2n)T算法设计：求P, theta初值Y 二 y;U = u；phil = E-y(2:49),-y(l:48),u(2:49),u(l:48);P = inv(phi * phil);P2 二 P;theta =E* phi *y(3:50);循环计算新theta值for i = 50 : length (u) - 2phi = -y(2i+l), -y(l:i),u(2:i+l), u(l:i);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) - phi