相关辨识法上机

1、上机实习IV-1：相关辨识法上机12自动化指导老师：王国利摘要：系统辨识、状态估计和控制理论是现代控制论中相互渗透的三个领域。控制理论的应 用离不开系统辨识技术，实际中，许多控制系统的模型在工作中是变化的，为了实现自适应 控制，需要系统辨识技术不断更新模型参数。通过学习使用MATLAB软件，初步体验系统 辨识方法。关键字：系统辨识，控制理论，MATLABoPractice 1： Practice the method ofSystem identificationAbsti-act:System identification, State estimation and Hie Principl

2、e of Automatic Control are tliiee different disciplines of the modem control theory, wliich are inteipenetrated with one another, hi practice, the model of system is changing all the time. To control adaptively, the system model should be update its parameters, by the method of System identification

3、. By learning the using of MATLAB, we are supposed to practice the method of system identification.Key Words: System identification, System identification, MATLAB目录 TOC o 1-5 h z 、弓I言3 HYPERLINK l bookmark11 o Current Document 1.1介绍3 HYPERLINK l bookmark14 o Current Document 1.2实验目的3 HYPERLINK l boo

4、kmark17 o Current Document 二、实验内容和方法3 HYPERLINK l bookmark20 o Current Document 2.1实验内容3 HYPERLINK l bookmark26 o Current Document 2.2实验原理3 HYPERLINK l bookmark32 o Current Document 一次完成算法4 HYPERLINK l bookmark44 o Current Document 2.2.2递推算法52.3实验步骤6一次完成算法6 HYPERLINK l bookmark57 o Current Document

5、2.3.2递推算法7 HYPERLINK l bookmark64 o Current Document 2.4实验程序和结果7 HYPERLINK l bookmark67 o Current Document 三、结果分析8引言1.1介绍在自然科学和社会科学的许多领域中，人们越来越重视对系统进行定量的系统分析、系 统综合、仿真、控制和预测。将被控对象模型化，是开展工作的前提和基础。基于传递函数的经典控制理论分析前需要确定模型的各个参数，忽略工作中的变化，适 用于精度要求不高的情况下。而实际中许多控制系统内部参数是未知的，而且工作中有可能 会变化，这就需要系统辨识来确定模型参数。通常为机理分

6、析法和测试法相结合来辨识，对 于系统机理己知的部分采用机理分析法，机理未知的部分采用测试法。1.2实验目的采用相关辨识方法，按要求完成脉冲响应函数的辨识。描述实验验证的数据准备、基本过程和实验结 果。使用上次实验产生的M伪随机信号序列，对日行设计的系统进行激励，产生输出信号，并采用一次完 成算法和递推算法进行系统辨识，并比较结果。二、实验内容和方法2.1实验内容采用相关辨识方法，按要求完成脉冲响应函数的辨识。描述实验验证的数据准备、基本过程和实验结 果。自己设定一个稳定系统，结合上机实习III中的M序列生成程序产生输入信号，获得系统输出值。利用题1中的数据，采用一次完成算法，估计脉冲响应序列。

7、利用题1中的数据，采用递推算法估计脉冲响应序列，并与题2的结果比较。2.2实验原理M序列是循环周期为Np*At,自相关函数近似于8函数的一种随机序列，其统计特性近似于 白噪声。M序列的自相关函数为：R(0)= 1 ； R(t)=-1/Np维纳-霍夫方程离散形式Ryu(k)=Zj=QtfX)g(/)Rull(k)At2.2.1 一次完成算法利用周期为Np的M序列d构造M(n)进而生成周期为Tp= NpAt连续的随机信号 u(t)9输入自相关函数(以下M替代u)R.(k)= a2(l +Np)/Npd(k)-a2/Np当Tp充分大时，近似表示成RyM0)=j=G,Np 套用M序列作输入信号的一次完

8、成算法，根据上面式子得到Np个线性联立方程：做(0)0(。)+砌(1)0 +.+曲(M+l)0(Np-l)=/S(O)W砌(1)(。)+做(0 娘+. +RmQNp+2)卧pl)=RyM(g Rmp-1)部)+Rmp-2)念 Q)+ . +Rm(。)食(Np-l)=R)，M0p-l)W_R展=玲，今自其中：Rm(O) &m(T)&m(-Np + 1)_Rm(1)Rm(O) (Np + 2)Rm-：: &M(-Np + 1) &M(-Np + 2) Rm(O)g=5（0） 0.RyM =RyM（O） .笊Np - 1）.RyM（Np - 1）.进而有:1M1M其中:M(0)M(l)MWM(0)M

9、(Np - 1)&M(Np 2)Y= y(o) yW . . . .M(Np + 1).M(Np + 2) M(0)yP-1).2.2.2递推算法第1时刻的自相关函数的估计值伙)为RyM(k)=j=QJp-y(f)M(j-k)/Np为了在线辨识系统的脉冲响应，将上式写成递推的形式。=g=o,a(/W (j-kyy(i)M O-/r)/(z+l)讷尸)(g(，W(的/(i+1)=RyM(f(k)i/Q+l)+yQ)MQk)/Q+l)则有欧=Rm=RiRy %Q+)+M(Oy(i)/(i+l)曷=g(/1)i/(z+l)+RA/1M(Oj0/(i+l)A其中M(i)=M(i) M(M) . M(i

10、-Np+1)2.3实验步骤2.3.1 一次完成算法1、按照上一次实验的结论生成一段长度为63位的M序列：Xl=l;X2=0;X3=l;X4=0; X5=0;X6=0;m=63;for i=l:mY6=X6;Y5=X5;Y4=X4; Y3=X3; Y2=X2;Y1=X1;X6=Y5;X5=Y4;X4=Y3; X3=Y2;X2=Y1;Xl=xor(Y6,Y5);ifY6=0U(i)=.l;elseU(i)=Y6;endendM=U;2、自己设计一个稳定的系统： num=2 1,den=l 5 4;G=tf(num,den)Transfer function2s+ 1 sA2 + 5 s + 43、

11、由于在实验步骤中要涉及到M序列的移位，可以用最普通的循环的办法来进行移位，但 是工作量大而旦容易弄错，这里我将M序列进行了拓展，拓展成三个周期，同时将激励响 应拓展成三个周期。w=MNM；Y1=Y,Y,Y;4、由于数据必须尽可能多地包含系统的动态特性信息，这和输入信号的选择有很大关系。 因此，实验之前要精心地确定心和fhiaxo通过对自行设定的传递函数的波特图和阶跃响应 分析，得出血ax约等于10hz, Ts约等于4s,所以At约等于0.033.dt=0.033;t=0.033:dt:63*dtr;Y=lsim(GMt);5、M矩阵的生成，选取第二个周期。for k=l:63Mm(k,: )=

12、w(l 26 -k: 18 8 -k);end6、得出结果：Gl=(l/(63+l)*l*0.033)*Q*Mm*Y;23-2递推算法23-2递推算法岸1、采用一次完成算法中的系统和M序列。2、由于在求G(2)的时候需要Ry (1)的值，所以需要用公式提前求得(选用第二个周期): for k=l:63Ry(k)=0.5*(w(126k 尸 Yl(126)+w(127k)*Yl(127); end3、由于在求G的初始值的时候，涉及到Ry (1-1)项，所以我从G(2)开始求起。 根据公式，同样选取第二个周期进行辨识，当K=2时：G2=(2/3*63*64)*Q*Ry/0.033+(l/3*63*

13、64)*Q*w(188:-l:126),*Yl(128)/0.033;其中：Q=eye(63)+ones(63,63);4、进入迭代：fort=3:63M_i=w(l,t+126:-1: tM4) G2=t/(t+l 尸 G2+1/0.033/(t+l) * Y1 (t+126 尸 Q* (M_i，尸63/64; end2.4实验程序和结果程序如下：Xl=l;X2=0;X3=l;X4=0; X5=0;X6=0;m=63;for i=l:mY6=X6;Y5=X5;Y4=X4; Y3=X3; Y2=X2;Y1=X1;X6=Y5;X5=Y4;X4=Y3; X3=Y2;X2=Y1;Xl=xor(Y6,

14、Y5);if Y6=0U=.1;elseU(i)=Y6;endendM=U;num=l 2;den=l 5 4;G=tf(num,den);Q=eye(63)+ones(63,63);w=MNM；for k=l:63Mm(k,:)=w(127 -k: 189-k);enddt=O.O15;t=0.015:dt:63*dtr;Y=lsim(GMt);Gl=(l/(63+l)*l*0.015)*Q*Mm*Y;%Y1=Y,Y,Y;for k=l:63Ry(k)=0.5*(w(188-k尸 Yl(126)+w(189k)*YlQ27);endG2=(2/3*63*64)*Q*Ry/0.033+(l/3*63*64)*Q*w(188:-l:126),*Yl(128)/0.033;fori=3:63M_i=w(i+188:l:i+126);G2=i