系统辨识课程报告.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+K

11、2*(y(i+2)-phi2,*theta);求解结果为：1231.48511.06701.14300.78200.46250.53710.49210.37600.44070.19560.19380.16745）用夏氏偏差修正法、夏氏改良法及其递推算法估计出夏氏偏差修正法:定义M阵：M=I-(),定义D阵：D=n(),=11f-(),Jn=M定义f:/=-D-,11,(,),ij+p,nj=D,n-(,)-,DIOTMy则有theta估计为：=(),1y(1)lf算法设计：定义M阵：T=InY(Phi*phi)*phi;M=eye(N)-phi*T;循环至收敛：while1e(3:N+2)=Y

12、phi*XSCtheta;omu=e(2：end-l),e(l:end-2);D=omu,*M*omu;f=Domu,*M*Y;theta=T*omu*f;ifabs(theta-theta)ZX12341.4881.38341.46280.79730.71460.77490.48220.38840.47110.20220.28730.225211夏氏改良法:将f用最小二乘重新估计：f=()-le再计算theta：=(),y(),11/算法设计：(修正上文中f的计算方法)e(3:N+2)=Y-phi*XSItheta;omu=-e(2:end-l),-e(l:end-2);f=(omu,*o

13、mu)omu,*e(3:N+2);XSItheta=LStheta-T*omu*f;结果：1231.48861.38341.46280.79730.71460.77490.48220.38840.47110.20220.28730.2252II递推夏氏法：定义beta阵：则beta的最小二乘估计为：(r)1y递推算法为：B、I=B+I(y+I-W1)P+I=PN-r、+1叭1P、N+=PWN+1(1+%+PI)1其中：W工=-y(%+N)-(N+1)M(7/+N+1)w(JV+1)一/(+N)e(+N+1-1)加1=+NT)算法设计：定义beta阵与P阵:p=(1E6)2*eye(6);bet

14、a=LStheta;O;0;迭代计算P,r,beta与thetafori-1:length(u)-2e=zeros(2+1,1);fai=-y(2:i+l),-y(l:i),u(2:i+l),u(l:i);e(3:i+2)=y(3:i+2)-fai*LStheta;phi=-y(i+l),-y(i),u(i+l),u(i),-e(i+l),-e(i)l,;r=p*phi(l+phi,*p*phi);p=p-r*phi,*p;beta=beta+r*(y(i+2)-phi,*beta);RXStheta=beta(1:4);end结果为：1231.48791.59021.86930.79700.

15、79520.96220.48220.39940.49450.20210.34180.33226)用增广矩阵法估计。；设theta为：0=aan60bnClePhi阵为：/T=一y(N1)一y(N)(+、)Z(.V)(?/+N-1)(N)则有：。i=+k、i(y*i*vziV)PH=P-K、，IKIP、K、产Px11(1+1PVi)1算法设计:fori=1:Nphi=-y(i+l),-y(i),u(i+l),u(i),-e(i+l),-e(i)l,;e(i+2)=y(i+2)-phi,*AMtheta;K=p*phi(l+phi,*p*phi);p=p-K*phi,*p;AMtheta=AMth

16、etaK*(y(i2)-phi,*AMtheta);end（P阵与前文定义相同）解算结果为：（仅取前四行为题目要求的theta）12311.48361.11741.122120.78640.51140.489730.48290.38420.431640.19910.18610.122655.8007e-04-0.0202-0.00996-0.0074-0.0285-0.01237）分析噪声式幻特性;对于Uyl数据：所有辨识结果如图:U2HlAMtheta|1.4836;0.7864;0.4829;0.199148王GYLStheta1.50350801604841;0.20473233IVth

17、eta1.5353;0.8332;0.4873;0.215632jLStheta1.4855;0.7869;0.4837;0.198232Hntheta1.9666J.2124;0.5702;0.283932王rankAIC115double120王rankerr115double120田rankF114double112王rankwhite115double120SRGYLStheta1.4851;0.7820;0.4921;0.1956322RlVtheta1.4861;0.7882;0.4843;0.1984320RLStheta1.4847;0.7865;0.4842;0.197832

18、SRXStheta1.4879;0.7970;0.4822;0.202132田SyS18战5OO1double4000田Uyl5002double8000Suy25002double8000三uy35002double8000田XSCtheta1.4886;0.7973;0.4822;0.202232由XSItheta1.4886;0.7973;0.4822;0.202232Sy5OO1double4000不难发现，最小二乘辨识结果LStheta与其余辨识方法所得结果差距较小（除步骤1所得的解析法求方程组再取均值的ntheta结果）。而最小二乘的辨识方法当且仅当噪声是白噪声时，其估计结果才是无

19、偏的。所以，当噪声是有色噪声时，广义最小二乘法、夏氏法、辅助变量法、增广矩阵法应当与最小二乘法相差一个量。故可以认为uyl样本数据中的噪声是白噪声。对于uy2与uy3数据:名称11字节9/Btheta|1.1174;0.5114;0.3842;0.186148丑IVtheta1.3965;0.7072;0.3802;0.297032丑LStheta1.1113;0.4963;0.3791;0.1879322theta1.6130;0.2895;-0.2678;0.694732卫rankAIC115double120卫rankerr1x15double120卫rankF114double112力

20、rankwhite115double120BRGYLStheta1.0670;0.4625;0.3760;0.193832丑RIVtheta1.1371;0.5283;0.3778;0.201932二RLStheta1.1123;0.4963;0.3797;0.1878325RXStheta1.5902;0.7952;0.3994;0.341832丑SyS18Bu5OO1double4000丑UyI5002double8000Suy25002double8000任uy35002double8000BXSCtheta1.3834;0.7146;0.3884;0.2873325XSltheta1.

21、3834;0.7146;0.3884;0.287332By5OO1double4000AMtheta1.1221;0.4897;0.4316;0.1226;-0.48由IVtheta1.6172;0.8687;0.4325;0.3138325LStheta1.1158;0.4801;0.4254;0.124532Btheta-3.3522;-2.9492;-2.7898;0.395732JrankAIC115double120Jrankerr115double120士rankF114double112Jrankwhite115double120BRGYLStheta1.1430;0.5371;

22、0.4407;0.167432由RIVtheta1.1672;0.5238;0.4264;0.144932丑RLStheta1.1168;0.4815;0.4261;0.124832RXStheta1.8693;0.9622;0.4945;0.332232sys18u500x1double4000uy15002double8000uy25002double8000uy35002double8000由XSCtheta1.4628;0.7749;0.4711;0.225232由XSItheta1.4628;0.7749;0.4711;0.225232田y5001double4000可以发现，最小二

23、乘与其他辨识方法的辨识结果存在较明显差距，可认为噪声为有色噪声。二.用极大似然法估计0。考虑上文分析，噪声很可能不是白噪声序列，故考虑带估参数为：(ImbQbnCCnJr则y的预测值为：Iy（八）=一一ally(k-)+5。（八）+bnu(k-)+ce(k-1)+c(A-TI)预测误差为：e(k)=y(k)-y(k)同时计算：J=e2(,Afi=da)计算hassian矩阵：现二号一迦In誓一(jTc,当口啥=7)一,叫户dbi筒idOitA2-e(ki)_VcgJl二阶偏导近似为：亚_e(k)e(k)a，幺28-，坨对theta进行估值：迭代至theta收敛。算法设计：j=u;ddJ=O;s

24、ig2=0;fori=3:length(u)sig2=sig2+e(i)2N;J=J+0.5*e(i)-2;epal(i)=y(i-l)-theta(5)*epal(i-l)-theta(5)*epaepa2(i)=y(i-2)-theta(5)*epa2(i-l)-theta(5)*epaepbl (i) = -u(i-l) epb2(i) = -u(i-2) epcl (i) = -e(i-l) epc2(i) = -e(i-2)-theta(5)*epbl(iel)-theta(5)*ep- theta(5)*epb2(i-l)-theta(5)*ep- theta(5)*epcl(i-

25、l)-theta(5)*ep- theta(5)*epc2(iel)-theta(5)*epde=epal(i),epa2(i),epbl(i),epb2(i),epcl(i),epc2(idj=dj+e(i)*de;ddj=ddj+de*de,;endtheta=theta-ddjdj;计算结果:123，1.46611.72231.7111;0.78850.87440.87900.47950.35140.34730.19700.44070.42780.24540.95590.97500.3514-0.00550.1152三.以上题的结果为例，进行：1 .分析比较各种方法估计的精度与计算量：a

26、N元方程解析法：精度极差，几乎没有参考性，没有综合总体数据计算，仅涉及4阶矩阵计算，计算量很小。b.最小二乘法及其递归算法：在噪声为白噪声时其精度较高，而当噪声为有色噪声时，最小二乘法会出现明显偏差。算法较为简单，当数据量大时涉及较大矩阵求逆运算，运算量增大，其递归算法可以避免大型矩阵求逆，减小计算量，运算结果与最小二乘几乎一致。c.辅助变量法及其递归算法：引入辅助变量修正有色噪声影响，但当噪声为白噪声时反倒导致了精度下降。同时其迭代特性使得运算量相对最小二乘法较大。d.广义最小二乘法及其递归算法：精度较高，对有色噪声与白噪声均能较好估计,*在算法迭代运行时，出现预测结果在多个值中反复跳变的

27、情况，使得结果不能收敛。同时，其迭代过程反复使用最小二乘法估计theta与f值，计算量很大。e.夏氏法及其递归算法：改善广义最小二乘的大计算量，去掉了f参数反复数据滤波，精度近似于广义最小二乘法。计算量较大f.增广矩阵法及其递归算法：精度较高，同时估计噪声参数，运算量一般。g.极大似然法：精度较高，计算量很大。2 .分析噪声方差的影响；噪声方差将很大程度影响不同辨识方法的辨识结果，当噪声的方差增大时，辨识结果与实际值的误差也会增大，甚至辨识结果将不具备可信度。而当噪声方差为O时，辨识结果即为真值。3 .比较白噪声和有色噪声对辨识的影响。白噪声易于辨识，此时最小二乘法具有无偏性与一致性，而当噪声

28、为有色噪声时，由于无法确定噪声的相关参数，对模型的辨识也逐渐困难，各类辨识方法也均采用对噪声估计、迭代的方法消除有色影响，辨识困难，结果误差大。四.系统模型阶次的辨识：1 .按残差方差定阶：取最小二乘估计：=()-lV计算残差e(k)=za(zl)y(k)b(zl)u(k)+NJ”=X。2(卜)A4I则，当n从1增大时，Jn随n的增大而减小，若no为正确的阶，则Jn图像将在no处为最后一次陡峭的下降，随后Jn图像将缓慢变化。算法设计：N=length(u)-i;Y(1:N)=y(i+l:N+i);forj=1:Nform=l:iphis(j,m)=-y(i+j-m);phXs(j,im)=u(

29、i+j-m);endendtheta=(phis,*phis)phis,*Y;err=Y-phis*theta;2.(i)=err,*err;结果作图为：其均在第二点完成一次陡峭的下降，故可认为模型阶次为2次。（由于噪声的影响，也可认为uy3的阶次为3）2.Akaike信息准则准则定义为：AIC=-21nL+2p计算噪声方差：InL三-yln2-yln?一步2()林=52()(k)=y(k)+5aiy(k-i)一b,u(k-i)-lcj(k-/)JI/0*IInL=-yIne+c选取不同阶数na,nb,nc,计算AIC：AIC=NlnJ+2(wu+nf)算法设计：theta=(phi,*phi

30、)phi,*y(n+l):500);e(n+l):500)=y(n+l):500)-phi*theta;j(n)=e*e,;AlC(n)=N*log(j(n)/N)+2*(na+nb);结果：3.按残差白色定阶:计算残差自相关函数:R(i)=Se(k)e(k+i)R()7)r(o)=ySj阳A)r()=算法设计：双重循环迭代，计算不同n不同i(算法中对应t)时的r值。r=zeros(10,1);r=err2(i)/N;fort=1:10fork=i+1:N+i-t-10r(t)=r(t)+err(k)*err(k+t);r=r/N;R(:,i)=abs(rr);结果；对于UyI数据：可以看到n

31、2时即白噪声化明显，自相关性迅速下降，故n=2对于uy2数据：同样可以看到n=2时即白噪声化明显，自相关性迅速下降，故n=2,当n=3时自相关函数反而上升，噪声更倾向有色。04对于uy3数据：同样可以看到n=2时即白噪声化明显，自相关性迅速下降，故n=2o即全部为2阶系统。2 .分析噪声对定阶的影响；可以发现，噪声的存在将会影响不同方法对模型阶次的判断，在残差方差定阶与AIC定阶中,uy3数据似乎也可以看作3阶系统，说明了噪声的影响不容忽略。3 .比较所用三种方法的优劣及有效性；残差方差定阶：数据陡峭凭主观判断AlC定阶：取AlC最小处，无歧义；残差白色定阶：直观清晰，但也存在一定主观因素四

32、给出由正弦输入求取系统开环频率响应特性曲线的辨识方法。原理：G(j) =U(j)Y为输出的傅里叶变换，U为输入的傅里叶变换。流程：采取组合频率的正弦信号作为测试信号，测量系统稳态下的输出波形，利用傅叶变换对输出的组合波进行分解。取输出y的一个周期作傅里叶级数展开：y(t)=A0+Z(A“cosnt+Bnsinnt)n=i其中：(由于实验通得到采样信号通常为离散点，故选取离散点的傅里叶变换式)1 N2 W2n.、Va)CoSCz1)NyNBn=VX0sin(/)NTN计算各个谐波振幅与相位:=arctan(Bz,AJ则有：U(j)=1u(t)ejtdt=(r)coscotdt-w(r)sint

33、diY(j)=y(r)e7ft*rfr=y(t)costdt-y(t)sintdt其中计算COSWt积分时：UC计算SinWt积分时：八1.+5)石采用数值积分法得到四个积分值：(示例一个)w(r)sin(OtdtJ-J。M,办(i)Kf)+K)T=IK；?=sinco(i-1)m5-2sincoitus+sin(i+I)rnK：?=sin(n-I)MS-sincuntlsAt=LUSn令：1P二LT1y(i)K2)+y()KWWya)Kz+M)K7I=IZwa)Kr+()?=WUa)KV+()?r=l则可得实频特性与虚频特性分别为:Re()=pr+qsr252ps-qr-22r+5六.提出一

34、种自己创造的辨识新方法，并用所给数据进行辨识验证。那我就命个名叫“好编程”法了首先通过依次循环4个方程求解theta的估计值（4元一次方程解析法）共计得到495个theta值theta4495double考虑问题1时，对这个495个theta取了均值，但效果及其不理想；考虑对theta每一行进行从小到大排序：（第一行数据示例）可以看到，其真值应当处在中间位置，而非数学意义上的按行平均值，算数平均使得两端的异常值也被计算在内。考虑剔除theta排序后第一行的前30%与后30%对剩余求平均即得theta估计值：田4x3double123411.54141.39971.447420.91810.87

35、030.890930.50440.36770.528140.24490.33160.24995对比极大似然法：与最小二乘法旧6x3double123411.46611.72231.711120.788三0.871J0.879C30.4795三0.3514H03473匚40.197C)|0.44070.427850.245/L0.955S0.975C60.35INI-0.00550.115243double1231.48551.11131.11580.78690.49630.48010.48370.37910.42540.19820.18790.1245-PS虽然没什么道理，不过感觉效果还行？.

36、像是最小二乘和非白色噪声极大似然的结合？附录：所有程序、子函数文件:Rl1_uy1.txt31_uy2.txt01uy3.txt2021/1/59:082021/1/59:082021/1/59:0819KB19KB19KB隹|L系统辩识作业题.doc2021/1/59:08MicrosoftWord.44KBDoIlO92083帚书签PdfAM.mConfirmRankAIC.mConfirmRanke.measy.mGYLS.mIV.mInL.m1.S.mmain.mn.mRGYLS.mRIV.mRLS.mRXS.mSC.mXSI.m2021/1/59:552021/1/616:092021/1/617:352021/1/610:122021/1/621:062021/1/615:442021/1/611:512021/1/616:392021/1/611:142021/1/621:062021/1/611:102021/1/615:492021/1/611:592021/1/611:192021/1/517:102021/1/516:412021/1/517:00MicrosoftEdge.M文件M文件M文件M文件M文件M文件M文件M文件M文件M文件M文件M文件M文件M文件M文件M文件9,460KB1KB1KB1KB1KB1KB1 KB