完整word版系统辨识实验手册

1、系统辨识实验手册哈尔滨工业大学 控制与仿真中心2012 年 8 月实验 1白噪声和 M 序列的产生实验脉冲响应法的实现实验最小二乘法的实现实验递推最小二乘法的实现附录实验报告模板121614实验1白噪声和M序列的产生、实验目的1、熟悉并掌握产生均匀分布随机序列方法以及进而产生咼斯白噪声方法2、熟悉并掌握M序列生成原理及仿真生成方法 二、实验原理1、混合同余法混合同余法是加同余法和乘同余法的混合形式，其迭代式如下:xn 1(a* xn b) mod MRn 1 Xn 1/M式中a为乘子，X0为种子，b为常数，M为模。混合同余法是一种递归算法, 即先提供一个种子X0，逐次递归即得到一个不超过模 M

2、的整数数列。2、正态分布随机数产生方法由独立同分布中心极限定理有：设随机变量Xi,X2,，Xn,相互独立，服从同一分布，且具有数学期望和方差:E(Xk),D(Xk)n则随机变量之和Xk的标准化变量：i 1nXki 1近似服从N(0,1)分布。nnXk E( Xk)i 1i 1Jd/xq如果Xn服从0, 1均匀分布，则上式中0.5，2 。即12nXk 0.5 nY j 1 L 近似服从N(0,1)分布。3、M序列生成原理4个双用移位寄存器产生M序列的简化框图如下图所示。该图表示一个由稳态触发器顺序连接而成的4级移位寄存器，它带有一个反馈通道。当移位脉冲 来到时，每级触发器的状态移到下一级触发器中

3、， 而反馈通道按模2加法规则反 馈到第一级的输入端。4三、实验内容1、生成均匀分布随机序列（1）利用混合同余法生成0,1区间上符合均匀分布的随机序列，并计算该序列的均值和方差，与理论值进行对比分析。要求序列长度为1200,推荐参数为 a=65539, M=2147483647, 0X0M。（2）将0, 1区间分为不重叠的等长的10个子区间，绘制该随机序列落在 每个子区间的频率曲线图，辅助验证该序列的均匀性。（3）对上述随机序列进行独立性检验。（该部分为选作内容）2、生成高斯白噪声利用上一步产生的均匀分布随机序列，令n=12,生成服从N（0,1）的白噪声,序列长度为100,并绘制曲线。3、生成M

4、序列M序列的循环周期取为Np 26 1 63，时钟节拍t 1Sec,幅度a 1，逻辑“C为a,逻辑“ 1为-a，特征多项式F（s） s6 s5。4. 生成的M序列为：00000100001100010100111101000111001001011011101100110101011111生成M序列的结构图如下所示。要求编写Matlab程序生成该M序列，绘制该信号曲线，并分析验证M序列 的性质。四、实验步骤1分别画出三部分实验内容的程序框图（流程图）2. 编制 MATLAB的M文件;3. 运行编制的M文件;4. 查看程序运行结果并进行分析;5. 填写实验报告。五、实验报告格式参见附录一。实验2

5、脉冲响应法、实验目的通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。、实验原理一个单入单出线性定常系统的动态特性可用它的脉冲响应函数g( R)描述。则 y(t) 0 g( )x(t)d上式两端同乘x(t1 T Tim - 0 y(t)x(t则)，)dt进而取时间均值，有1 T0 g( )ITm T 0 x(t )x(t )dtd)Rx()d这就是著名的维纳霍夫积分方程。Ry()0 g(X(t)L线性系紡-g何如果输入是 白噪声，这时x(t)的自相关函数为Rx( ) k ( ),Rx(则根据维纳霍夫积分方程可得Rxy()0g( )Rx()d kg()kx(t)、y(t)的值，并计算它们的

6、互相关函数，即可求得脉冲或者这样，只要记录响应函数g( To)而在系统有正常输入的情形下，辨识脉冲响应的原理图如下图所示。线性系统乘法赛一积分踣SB三、实验内容下图为本实验的原理框图。系统的传递函数为G(s),其中K 120,Ti 8.3Sec,T2 6.2Sec ; u(k)和z(k)分别为系统的输入和输出变量；v(k)为测量白噪声，服从正态分布，均值为零，方差为；，记作v(k) N(0, 2) ； g0(k)9为系统的脉冲响应理论值，g(k)为系统脉冲响应估计值，g(k)为系统脉冲响应 估计误差。系统的输入采用M序列(采用实验1中的M序列即可)，输出受到白噪声v(k)的污染。根据过程的输入

7、和输出数据u(k),z(k)，利用相关分析法计算出系9统的脉冲响应值g(k)，并与系统的脉冲响应理论值g0(k)比较，得到系统脉冲响 应估计误差值(k)，当k 时，应该有(k)0。1模拟过程传递函数G(s)，获得过程的输入和输出数据u(k), z(k)(采样? g(k)k t/T?时间取1 秒)。(1)惯性环节u(k)Ky(k)s 1/T其中，T为惯性环节的时间常数，K为惯性环节的静态放大倍数。若采样时间记 作To，则惯性环节的输出可写成：y(k) eT0/Ty(k 1)TKT(eT0/TTK(1 e T0/T)u(k 1)1) T0)u(k) u(k 1)To传递函数G(s)仿真(串联)G(

8、s)T1T2 s 1/s 1/T2令 K1= KT1T2则G (s)的表达框图为:u(k)K1x(k)11 y(k).s 1/T1s 1/T2I2、互相关函数的计算1 (r 1)NpRMz(k) u(i k)z(i)rNp i Np 1其中，r为周期数，i Np 1表示计算互相关函数所用的数据是从第二个周期开始的，目的是等过程仿真数据进入平稳状态。(可分别令r =1、3,对比仿真结果)3、c的补偿补偿量c应取Rmz(N P 1)，不能取Rmz( Np)。因为RMz(k)是周期函数,则有 Rmz( Np)= Rmz(0)，故不能取一Rmz(N p)。4、计算脉冲响应估计值?脉冲响应估计值g(k)

9、Rmz ( k) CNp脉冲响应估计误差 g四、实验步骤(1)掌握相关分析辨识方法的基本原理;设计实验方案，画出程序框图;编制实验程序;调试并运行程序，记录数据;(5) 分析实验结果，完成实验报告。五、实验报告格式参见附录一。实验 3 最小二乘法的实现、实验目的理解并掌握系统辨识中的最小二乘法原理。、实验原理给定系统y(k)a1y(k 1) a2y(k 2) Lany(k n)(k)b0u(k) b1u(k 1) L bnu(k n)其中 a1,a2,L,an, b0,b,b2,L ,bn为待辨识的未知参数,(k) 是不相关随机序列。y为系统的输出，U为系统的输入。分别测出n N个输出、 n

10、N 输入值y(1),y(2), y(3),L y(n N),u(1),u(2),Lu(n N) ，则可写出N 个方程,具体写成矩阵形式,有y(n y(n My(n1)2)N)a1Manb0Mbn(n(nM(n1)2)N)a1y(n1)y(n)Ly(1)U(n1)LU(1)M(n1)y(n2)y(n 1)Ly(2)U(n2)LU(2)an(n2)MMMMMMb0My(nN)y(n N 1)Ly(N)U(nN)LU(N)M(nN)bn则式( 2)可写为y(n)Ly(1)U(n1)LU(1)y(n 1)Ly(2)U(n2)LU(2)MMMMMy(n N 1)Ly(N)U(nN)LU(N)式中：y为N

11、维输出向量； 为N维噪声向量；为2n 1维参数向量;N (2n 1) 测量矩阵。为了尽量减小噪声对 估值的影响，应取 N 2n1，即方程数目大于未知数数目。的最小二乘估计为三、实验内容对象的数学模型如下：z(k) 1.5z(k 1) 0.7z(k 2)u(k 1) 0.5u(k 2) v(k)其中， v(k) 是服从正态分布的白噪声N(0,1)。输入信号采用4阶M序列,幅度为 1。选择如下形式的辨识模型：z(k) a1z(k 1) a2z(k 2)b1u(k 1) b2u(k 2) v(k)设输入信号的取值是从k =1 到 k=16的M序列，则待辨识参数?Ls为?Ls=(HTHl) 1H Tz

12、l。其中,被辨识参数?LS、观测矩阵zL、H L 的表达式为a1?Ls：2b1b2zLz(3)z(4)z(16)HLz(2)z(3)z(1)z(2)u(2)u(3)u(1)u(2)z(15)z(14)u(15)u(14)要求编制仿真程序，获取系统输入输出数据，并运用最小二乘法对这一系统 的参数进行辨识，并将辨识结果与实际参数进行对比。四、实验步骤1 写出系统结构、实际参数、噪声源及输入信号等内容；2画出程序框图(流程图) ；3编制 MATLAB 的 M 文件；4运行编制的 M 文件；5查看程序运行结果并进行分析；6填写实验报告。五、实验报告格式参见附录一。实验 4 递推最小乘法的实现、实验目的

13、熟悉并掌握递推最小二乘法的算法原理。、实验原理给定系统y(k) a1y(k 1) a2y(k2) L any(k n)bnu(k n) (k)( 1)(k) 是不相关随机序列。 n N 个输出、 n N 输入值 则可写出 N 个方程， 具体写成矩则式( 2)可写为式中： y 为 N 维输出向量；维噪声向量； 为 2n1维参数向量；3)为b0u(k) b1u(k 1) L其中ai,a2丄,an, boQdL ,bn为待辨识的未知参数, y 为系统的输出， u 为系统的输入。分别测出 y(1),y(2), y(3),L y(n N),u(1),u(2),L u(n N) ， 阵形式，有a1y(n1

14、)y(n)Ly(1)u(n1)Lu(1)M(n1)y(n2)y(n 1)Ly(2)u(n2)Lu(2)an(n2)MMMMMMb0My(nN)y(n N 1)Ly(N)u(nN)Lu(N)M(nN)bnyMJb0 ,My(n N)M(nN)bny(n)Ly(1)u(n 1)Lu(1)y(n 1)Ly(2)u(n 2)Lu(2)MMMMMy(n N 1)Ly(N)u(n N)Lu(N)1)2)1)2)y(n y(n(n (na1ManN （2n 1）测量矩阵。为了尽量减小噪声 方程数目大于未知数数目。的最小二乘估计为 ）（T ）1对估值的影响，应取N 2n 1,即(4)为了实现实时控制，必须采用

15、递推算法，这种辨识方法主要用于在线辨识。 设已获得的观测数据长度为 N,将式（3）中的y、 替，即和分别用Yn，N， N来代Yn(5)用）N表示的最小二乘估计，则）NnYn(6)1N ，则Pn(7)如果再获得一组新的观测值u(n1)和 y(nN 1），则又增加一个方程式中yN 1yN 1 y(n1),1 (nN 1)y(n N) Ly(N 1)u(n N 1) L u(N1)将式（5）和式（8）合并，并写成分块矩阵形式，可得YnL LyN 1NL LTN 1L LN 1类似地可得到新的参数估值T1)N 1NNNL LLLL LTN 1TN 1TN 1TNYnPN 1L LL LTN 1yN1P

16、n 1N YnN 1y 1T于是,YnL LyN 1(9)(10)16式中PN 1NLLTN1PN1T N 1 N11NLLTN1T1N1111)应用矩阵求逆引理，从求得 PN1与Pn的递推关系式出发，经过一系列的推导,最终可求得 递推最小二乘法辨识公式 ：)N1KN 1 yN)1N12)KN 1PN1PNN113)N1 1Pn和)N的初值Po和推荐取值方法为：假定 )0 0,P0 c2I ， c 是充分大的常数，(2n 1) (2n 1)单位矩阵，则经过若干次递推之后能得到较好的参数估计。PN 1PNPNN 1PN为了进行递推计算，需要给出N1)0N 1PN14)三、实验内容给定系统y(k) a1y(k 1) a2y(k 2) 即 n 2。假设实际系统的参数为 但是不已知，即不可测。取 (k)b0u(k) b1u(k 1) b2u(k 2) (k)15)a1 2 ，a2 1.3 ，b0 0.4 ， b1 0.88 ， b2 0.1,0.1 的零均值随机序列。输入信号取为2.2，231245u(k) 1.5sin 0.2k(16)要求编制 MATLAB 程序，运用递推最小二乘法对