2022年2003版系统辨识最小二乘法大作业

1、西 北 工 业 大 学系统辩识 大作业题目：最小二乘法系统辨识可编辑资料 - - - 欢迎下载一， 问题重述：用递推最小二乘法，加权最小二乘法，遗忘因子法，增广最小二乘法，广义最小二乘法，帮忙变量法辨识如下模型的参数GYUs4s46.5s320s314s2120s211.5s3232s320离散化有z4 - 3.935 z3 + 5.806 z2 - 3.807 z + 0.9362=z4 - 3.808 z3 + 5.434 z2 - 3.445 z + 0.8187噪声的成形滤波器NEWs420s31120s2232s320离散化有4.004e-010 z3 + 4.232e-009 z2

2、 + 4.066e-009 z + 3.551e-010=z4 - 3.808 z3 + 5.434 z2 - 3.445 z + 0.8187采样时间 0.01s要求： 1. 用 Matlab 写出程序代码.2. 画出实际模型和辨识得到模型的误差曲线.3. 画出递推算法迭代时各辨识参数的变化曲线.最小二乘法 :在系统辨识领域中, 最小二乘法是一种得到广泛应用的估量方法, 可用于动态 , 静态 ,线性 , 非线性系统. 在使用最小二乘法进行参数估量时, 为了实现实时把握, 必需优化成参数递推算法 , 即最小二乘递推算法.这种辨识方法主要用于在线辨识.MATLAB是一套高性能数字运算和可视化软件

3、, 它集成概念设计, 算法开发 , 建仿照真 , 实时实现于一体, 构成了一个使用便利， 界面友好的用户环境, 其强大的扩展功能为各领域的应用供应了基础.对可编辑资料 - - - 欢迎下载于一个简洁的系统, 可以通过分析其过程的运动规律, 应用一些已知的定理和原理, 建立数学模型 , 即所谓的“白箱建模 ”.但对于比较复杂的生产过程, 该建模方法有很大的局限性.由于过程的输入输出信号一般总是可以测量的, 而且过程的动态特性必定表现在这些输入输出数据中 , 那么就可以利用输入输出数据所供应的信息来建立过程的数学模型. 这种建模方法就称为系统辨识. 把辨识建模称作 “黑箱建模” .系统辨识又分为参

4、数辨识和阶次辨识 , 在本文中只争辩参数辨识问题最小二乘递推算法所用的模型： Zk=B uk+vk最小二乘递推算法为 :可编辑资料 - - - 欢迎下载ukvkGz 1 +ykN z 1 ekzk可编辑资料 - - - 欢迎下载图 1 最小二乘递推算法辨识实例结构图可编辑资料 - - - 欢迎下载vk 是听从 N0,1分布的不相关随机噪声.可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载G z 1 B z 1 1, N z 1 D z 1 ,1可编辑资料 - - - 欢迎下载A zC z考虑如图 1 示仿真对象 , 系统的差分方程为zk=3.808*zk-1-5.434*zk-

5、2+3.445*zk-3-0.8187*zk-4+uk-3.935*uk-1+5.8可编辑资料 - - - 欢迎下载06*uk-2-3.807*uk-3+0.9362*uk-4+v k 3.69可编辑资料 - - - 欢迎下载仿真对象选择如下的模型结构可编辑资料 - - - 欢迎下载zka1 zk1a2 zk2a3zk3a4zk4可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载b1ukb2uk1b3uk2b4uk3b5uk43.68vk可编辑资料 - - - 欢迎下载其中, v k 是听从正态分布的白噪声N0,1 .输入信号接受4 位移位寄存器产生的M序列,可编辑资料 - -

6、- 欢迎下载幅度为 1.按式 3.69构造 h k .加权阵取单位阵 LI .利用式 3.61运算 K k ， . k可编辑资料 - - - 欢迎下载和 P k ,运算各次参数辨识的相对误差,精度中意要求式3.67后停机.最小二乘递推算法辨识的 Malab6.0 程序流程图：可编辑资料 - - - 欢迎下载最小二乘递推算法辨识程序clear % 清理工作间变量L=35; % M 序列的周期y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值for i=1:L;%开头循环,长度为Lx1=xory3,y4; % 第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”x2=

7、y1; % 其次个移位寄存器的输入是第一个移位寄存器的输出x3=y2; % 第三个移位寄存器的输入是其次个移位寄存器的输出x4=y3; % 第四个移位寄存器的输入是第三个移位寄存器的输出yi=y4; %取出第四个移位寄存器的幅值为"0" 和"1" 的输出信号,即 M序列if yi>0.5,ui=-1; %假如 M序列的值为 "1",辨识的输入信号取“ -1 ” else ui=1; %假如 M序列的值为 "0",辨识的输入信号取“ 1”end %小循环终止y1=x1;y2=x2;y3=x3;y4=x4; %为

8、下一次的输入信号做预备end %大循环终止,产生输入信号u figure1; %第一个图形stemu,grid on %显示出输入信号径线图并给图形加上网格z2=0;z1=0; z3=0;z4=0;%设 z 的前四个初始值为零for k=5:25; %循环变量从 5 到 25zk=3.808*zk-1-5.434*zk-2+3.445*zk-3-0.8187*zk-4+uk-3.935*uk-1+5.806*uk-2-3.807*uk-3+0.9362*uk-4; %输出采样信号end%RLS递推最小二乘辨识c0=0.001 0.001 0.001 0.001 0.001 0.001 0.00

9、1 0.001 0.001' %直接给出被辨识参数的初始值 , 即一个充分小的实向量p0=106*eye9,9; %直接给出初始状态P0,即一个充分大的实数单位矩阵E=0.000000005; % 取相对误差 E=0.000000005c=c0,zeros9,24; %被辨识参数矩阵的初始值及大小e=zeros9,25; %相对误差的初始值及大小for k=5:25; %开头求 Kh1=-zk-1,-zk-2,-zk-3,-zk-4,uk,uk-1,uk-2,uk-3,uk-4'x=h1'*p0*h1+1; x1=invx; %开头求 Kk k1=p0*h1*x1;%求

10、出 K 的值d1=zk-h1'*c0; c1=c0+k1*d1; %求被辨识参数 c e1=c1-c0; %求参数当前值与上一次的值的差值e2=e1./c0; %求参数的相对变化e:,k=e2; %把当前相对变化的列向量加入误差矩阵的最终一列c0=c1; %新获得的参数作为下一次递推的旧参数c:,k=c1; %把辨识参数 c 列向量加入辨识参数矩阵的最终一列可编辑资料 - - - 欢迎下载p1=p0-k1*k1'*h1'*p0*h1+1; %求出 pk的值p0=p1; %给下次用if e2<=E break; %假如参数收敛情形中意要求,终止运算end %小循环终

11、止end %大循环终止c,e%显示被辨识参数及其误差 收敛 情形%分别参数a1=c1,:;a2=c2,:;a3=c3,:;a4=c4,:;b1=c5,:;b2=c6,:;b3=c7,:; b4=c8,:; b5=c9,:;ea1=e1,:;ea2=e2,:;ea3=e3,:;ea4=e4,:;eb1=e5,:;eb2=e6,:; eb3=e7,:; eb4=e8,:; eb5=e9,:;figure2; %其次个图形i=1:25; %横坐标从 1 到 25 ploti,a1,'r',i,a2,':',i,a3,'r',i,a4,'k

12、9;,i,b1,'y',i,b2,':',i,b3,'m',i,b4,':',i,b5,'g' %画出 a1, a2, a3,a4,b1 , b2,b3,b4,b5的各次辨识结果title'参数变化曲线 ' % 图形标题figure3; %第三个图形i=1:25; %横坐标从 1 到 25 ploti,ea1,'r',i,ea2,'k:',i,ea3,'y',i,ea4,'m',i,eb1,'g',i,eb2,'

13、c:',i,eb3,'r',i,eb4,':',i,eb5,'g' %画出 a1, a2, a3,a4,b1 , b2,b3,b4,b5的各次辨识结果的收敛情形title'误差曲线'%图形标题可编辑资料 - - - 欢迎下载参数变化图相对误差图仿真结果说明, 大约递推到第十五步时, 参数辨识的结果基本达到稳固状态,a2=5.434,a 3=-3.445,a4= -0.8187,b1=1,b2=-3.935b 3=5.806 b 4=-3.807即 a1=-3.808,b 5=0.9362 .此时,参数的相对变化量E0.00

14、0000005 .从整个辨识过程来看,精度的要求直接影响辨识的速度. 虽然最终的精度可以达到很小,但开头阶段的相对误差会很大,从相对误差图形中可看出,参数的最大相对误差会达到三位数可编辑资料 - - - 欢迎下载增 广 最 小 二 乘 法 算 法 所 用 的 模 型 ：Zk=Buk+Dek增广最小二乘法算法为：模型结构选用如下形式：zka1zk1a2zk2a3zk3a4z k4b1ukb2uk1b3u k2b4uk3b5uk4d1* vk -1d2 * vk - 2d3 * vk - 3d4* vk - 4可编辑资料 - - - 欢迎下载增广最小二乘算法流程图如下页图所示可编辑资料 - - -

15、 欢迎下载增广最小二乘辨识程序clearL=55; %4 位移位寄存器产生的M序列的周期y1=1;y2=1;y3=1;y4=0; %4个移位寄存器的输出初始值for i=1:L;x1=xory3,y4; %第一个移位寄存器的输入信号x2=y1; %其次个移位寄存器的输入信号x3=y2; % 第三个移位寄存器的输入信号x4=y3; %第四个移位寄存器的输入信号yi=y4; %第四个移位寄存器的输出信号,M序列 ,幅值"0" 和"1" , if yi>0.5,ui=-1; %M序列的值为 "1" 时, 辨识的输入信号取“ -1 ”

16、else ui=1; %M序列的值为 "0" 时, 辨识的输入信号取“ 1”endy1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号预备endfigure1; %画第一个图形subplot2,1,1; %画第一个图形的第一个子图stemu,grid on %画出 M序列输入信号v=randn1,60; %产生一组 60 个正态分布的随机噪声subplot2,1,2; %画第一个图形的其次个子图plotv,grid on; %画出随机噪声信号u ,v % 显示输入信号和噪声信号z=zeros13,55; %输出采样矩阵z2=0; z1=0;z3=0; z4

17、=0;%输出采样的初值c0=0.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.0010.001' %直接给出被辨识参数的初始值, 即一个充分小的实向量p0=106*eye13,13; %直接给出初始状态P0,即一个充分大的实数单位矩阵E=0.00000000005; % 相对误差 E=0.000000005c=c0,zeros13,54; %被辨识参数矩阵的初始值及大小e=zeros13,55; %相对误差的初始值及大小for k=5:55; %开头求 Kzk=3.808*zk-1-5.434*zk-2+3.445*zk-

18、3-0.8187*zk-4+uk-3.935*uk-1+5.8 06*uk-2-3.807*uk-3+0.9362*uk-4+4.004e-010*vk-1+4.232e-009*vk-2+4.066e-009*vk-3+3.551e-010*vk-4; %系统在 M序列输入下的输出采样信号h1=-zk-1,-zk-2,-zk-3,-zk-4,uk,uk-1,uk-2,uk-3,uk-4,vk-1,vk-2,vk-3,vk-4' %为求 Kk 作预备x=h1'*p0*h1+1; x1=invx; k1=p0*h1*x1; %Kd1=zk-h1'*c0; c1=c0+k1

19、*d1; %辨识参数 c e1=c1-c0; e2=e1./c0; %求参数误差的相对变化e:,k=e2;可编辑资料 - - - 欢迎下载c0=c1; %给下一次用c:,k=c1; %把递推出的辨识参数c 的列向量加入辨识参数矩阵p1=p0-k1*k1'*h1'*p0*h1+1; %find pk p0=p1; %给下次用if e2<=E break; %如收敛情形中意要求,终止运算end %判定终止end %循环终止c, e, %显示被辨识参数及参数收敛情形%分别变量a1=c1,:;a2=c2,:;a3=c3,:;a4=c4,:;b1=c5,:;b2=c6,:;b3=c

20、7,:;b4=c8,:; b5=c9,:; %分别出 a1,a2, a3,a4,b1 , b2,b3,b4,b5d1=c10,:; d2=c11,:; d3=c12,:; d4=c13,:;%分别出 d1， d2 ， d3 d4 ea1=e1,:;ea2=e2,:;ea3=e3,:;ea4=e4,:;eb1=e5,:;eb2=e6,:;eb3=e7,:; eb4=e8,:; eb5=e9,:; %分别出 a1, a2, a3,a4,b1 , b2,b3,b4,b5的收敛情形ed1=e10,:; ed2=e11,:; ed3=e12,:;ed4=e13,:; %分别出 d1 ，d2 ，d3 d4

21、的收敛情形figure2; %画其次个图形i=1:55;ploti,a1,'r',i,a2,':',i,a3,'r',i,a4,'k',i,b1,'y',i,b2,':',i,b3,'m',i,b4,':',i,b5,'g',i,d1,'k',i,d2,'g:',i,d3,'m',i,d4,'r' %画出各个被辨识参数title'参数变化曲线 ' % 标题figure3;

22、i=1:55; %画出第三个图形ploti,ea1,'r',i,ea2,'k:',i,ea3,'y',i,ea4,'m',i,eb1,'g',i,eb2,'c:',i,eb3,'r',i,eb4,':',i,eb5,'g',i,ed1,'y',i,ed2,'g:',i,ed3,'k',i,ed4,'m' %画出各个参数收敛情形title'误差曲线 ' % 标题可编辑资料 -

23、 - - 欢迎下载参数变化曲线相对误差曲线可编辑资料 - - - 欢迎下载仿真结果说明, 递推到第十步时, 参数辨识的结果基本达到稳固状态,即 a1=-3.808,a2=5.434, a3=-3.445,a4= -0.8187,b1=1,b2=-3.935b 3 =5.806 b 4=-3.807b 5=0.9362 , d1=0.0000,d2=0.0000,d3=0.0000 , d4=0.0000 .此时,辨识参数的相对变化量E0.000000005 .与最小二乘递推算法相比, 增广最小二乘递推算法虽然考虑了噪声模型,同样具有速度快， 辨识结果精确的特点.广义最小二乘法所用的模型是：Zk

24、=Buk+ek可编辑资料 - - - 欢迎下载cz kz k1zk2b u k1b u k2e k可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载aa1212e kc e k1e k2k可编辑资料 - - - 欢迎下载12其中我们选定模型参数：a1=1.5,a2=-0.7,b1=1,b2=0.5,c1=0,c2=0广义最小二乘递推算法的运算步骤：可编辑资料 - - - 欢迎下载.0充分小的实向量 可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载1. 给定初始条件0Pfe. 0a I a为充分大的数 20可编辑资料 - - - 欢迎下载ezfP 0I可编辑资料 -

25、 - - 欢迎下载2. 利用式 kzuffkCzCz z k11 zk,运算 k 和 ufk .可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载1na1nb3. 利用式 a ,a ,b ,b ,构造h f k .可编辑资料 - - - 欢迎下载fhKffhf k zf k1, zf kna ,uf k1,uf knb 可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载.k .k1 k zkk. k1可编辑资料 - - - 欢迎下载ffKf4. 利用式kk1hkkP k1hk1递推运算. k .可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载Pfhf

26、fPKfff1f kIkhk P k1可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载5. 利用e. kzkh k. k 和可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载hkzk1,zkna , u k1,u knb 运算 e.k .可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载6. 依据h e ke. k1,e. knc ehe来构造h e k .可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载ee.k. k1k e.kk.k1可编辑资料 - - - 欢迎下载可编辑资料 - - - 欢迎下载7.利用K e kPe k1he

27、khe k Pe k1h e k1可编辑资料 - - - 欢迎下载PeKe1eKeekIkhkP k1可编辑资料 - - - 欢迎下载广义最小二乘法程序代码clear % 清理工作间变量L=55; % M 序列的周期y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值for i=1:L;%开头循环,长度为Lx1=xory3,y4; % 第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”x2=y1; % 其次个移位寄存器的输入是第一个移位寄存器的输出x3=y2; % 第三个移位寄存器的输入是其次个移位寄存器的输出x4=y3; % 第四个移位寄存器的输入是第三个移位寄

28、存器的输出yi=y4; %取出第四个移位寄存器的幅值为"0" 和"1" 的输出信号,即 M序列if yi>0.5,ui=-1; %假如 M序列的值为 "1",辨识的输入信号取“ -1 ” else ui=1; %假如 M序列的值为 "0",辨识的输入信号取“ 1”end %小循环终止y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做预备end %大循环终止,产生输入信号u z2=0;z1=0; %设 z 的前四个初始值为零for k=3:45; %循环变量从 5 到 45zk=1.5*z

29、k-1-0.7*zk-2+1*uk-1+0.5*uk-2; %输出采样信号end%RGLS广义最小二乘辨识c0=0.001 0.001 0.001 0.001' %直接给出被辨识参数的初始值, 即一个充分小的实向量pf0=106*eye4,4; %直接给出初始状态P0,即一个充分大的实数单位矩阵ce0=0.001 0.001'pe0=eye2,2;c=c0,zeros4,44; %被辨识参数矩阵的初始值及大小ce=ce0,zeros2,44;e=zeros4,45; %相对误差的初始值及大小ee=zeros2,45;for k=3:45; %开头求 K zfk=zk+0*zk-

30、1+0*zk-2;ufk=uk+0*uk-1+0*uk-2;hf1=-zfk-1,-zfk-2,ufk-1,ufk-2' x=hf1'*pf0*hf1+1; x1=invx; %开头求 Kk k1=pf0*hf1*x1;%求出 K 的值d1=zfk-hf1'*c0; c1=c0+k1*d1; %求被辨识参数 ce1=c1-c0; %求参数当前值与上一次的值的差值e2=e1./c0; %求参数的相对变化e:,k=e2; %把当前相对变化的列向量加入误差矩阵的最终一列c0=c1; %新获得的参数作为下一次递推的旧参数c:,k=c1; %把辨识参数 c 列向量加入辨识参数矩阵的最终一列可编辑资料 - - - 欢迎下载pf1=pf0-k1*hf1'*pf0;