《系统辨识》实验手册哈尔滨工业大学控制与仿真中心2012年8月目录实验1白噪声和M序列的产生----------------------------------------------------------2实验2脉冲响应法的实现----------------------------------------------------------------5实验3最小二乘法的实现---------------------------------------------------------------9实验4递推最小二乘法的实现----------------------------------------------------------12附录实验报告模板----------------------------------------------------------------------16实验1白噪声和M序列的产生一、实验目的1、熟悉并掌握产生均匀分布随机序列方法以及进而产生高斯白噪声方法2、熟悉并掌握M序列生成原理及仿真生成方法二、实验原理1、混合同余法混合同余法是加同余法和乘同余法的混合形式,其迭代式如下:111(*)mod/nnnnxaxbMRxM式中a为乘子,0x为种子,b为常数,M为模。混合同余法是一种递归算法,即先提供一个种子0x,逐次递归即得到一个不超过模M的整数数列。2、正态分布随机数产生方法由独立同分布中心极限定理有:设随机变量12,,....,,...nXXX相互独立,服从同一分布,且具有数学期望和方差:2(),()0,(1,2,...)kkEXDXk则随机变量之和1nkiX的标准化变量:1111()()nnnkkkiiinkiXEXXnYnDX近似服从(0,1)N分布。如果nX服从[0,1]均匀分布,则上式中0.5,2112。即10.512nkiXnYn近似服从(0,1)N分布。3、M序列生成原理用移位寄存器产生M序列的简化框图如下图所示。该图表示一个由4个双稳态触发器顺序连接而成的4级移位寄存器,它带有一个反馈通道。当移位脉冲来到时,每级触发器的状态移到下一级触发器中,而反馈通道按模2加法规则反馈到第一级的输入端。三、实验内容1、生成均匀分布随机序列(1)利用混合同余法生成[0,1]区间上符合均匀分布的随机序列,并计算该序列的均值和方差,与理论值进行对比分析。要求序列长度为1200,推荐参数为a=65539,M=2147483647,0x0M。(2)将[0,1]区间分为不重叠的等长的10个子区间,绘制该随机序列落在每个子区间的频率曲线图,辅助验证该序列的均匀性。(3)对上述随机序列进行独立性检验。(该部分为选作内容)2、生成高斯白噪声利用上一步产生的均匀分布随机序列,令n=12,生成服从N(0,1)的白噪声,序列长度为100,并绘制曲线。3、生成M序列M序列的循环周期取为63126PN,时钟节拍Sec1t,幅度1a,逻辑“0”为a,逻辑“1”为-a,特征多项式65()Fsss。生成M序列的结构图如下所示。要求编写Matlab程序生成该M序列,绘制该信号曲线,并分析验证M序列的性质。四、实验步骤1.分别画出三部分实验内容的程序框图(流程图);2.编制MATLAB的M文件;3.运行编制的M文件;4.查看程序运行结果并进行分析;5.填写实验报告。五、实验报告格式参见附录一。C1C2C3C4C5C6CPM(6)M(5)+M(4)M(3)M(2)M(1)M(0)实验2脉冲响应法一、实验目的通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。二、实验原理一个单入单出线性定常系统的动态特性可用它的脉冲响应函数g(σ)来描述。这样,只要记录x(t)、y(t)的值,并计算它们的互相关函数,即可求得脉冲响应函数g(τ)。而在系统有正常输入的情形下,辨识脉冲响应的原理图如下图所示。0()()()ytgxtd则000()11lim()()(){lim()()}TTTTxtytxtdtgxtxtdtdTT上式两端同乘,进而取时间均值,有0()()()xyxRgRd则这就是著名的维纳霍夫积分方程。0()()(),()()()()()()()()xxxyxxyxtRkRkRgRdkgRgk如果输入是,这时的自相关函数为则根据维纳霍夫积分方程可得或者白噪声三、实验内容下图为本实验的原理框图。系统的传递函数为)(sG,其中Sec26TSec,3812021..,TK;)()(kzku和分别为系统的输入和输出变量;)(kv为测量白噪声,服从正态分布,均值为零,方差为2v,记作),(~)(20vNkv;)(kg0为系统的脉冲响应理论值,)(ˆkg为系统脉冲响应估计值,)(~kg为系统脉冲响应估计误差。系统的输入采用M序列(采用实验1中的M序列即可),输出受到白噪声)(kv的污染。根据过程的输入和输出数据)(),(kzku,利用相关分析法计算出系统的脉冲响应值)(ˆkg,并与系统的脉冲响应理论值)(kg0比较,得到系统脉冲响应估计误差值)(~kg,当k时,应该有0)(~kg。1、模拟过程传递函数)(sG,获得过程的输入和输出数据)(),(kzku(采样时间取1秒)。(1)惯性环节其中,T为惯性环节的时间常数,K为惯性环节的静态放大倍数。若采样时间记作0T,则惯性环节的输出可写成:相关分析法v(k)u(k)z(k)ckRtaNNkgMzPP)()()(ˆ2121210TtkTtkeeTTKkg//)()(ˆ)()(~kgkgkg0))(()(1121sTsTKsGTsK/1u(k)y(k)0011111000TkukuTeTTKkueTKkyekyTTTTTT)()())()()()()(///(2)传递函数)(sG仿真(串联)21211111TsTsTTKsG//)(令211TTKK=,则)(sG的表达框图为:2、互相关函数的计算PPNrNiPMzizkiurNkR)()()()(111其中,r为周期数,1PNi表示计算互相关函数所用的数据是从第二个周期开始的,目的是等过程仿真数据进入平稳状态。(可分别令r=1、3,对比仿真结果)3、c的补偿补偿量c应取)(1PMzNR-,不能取)(PMzNR-。因为)(kRMz是周期函数,则有)()(0MzPMzRNR=,故不能取)(PMzNR-。4、计算脉冲响应估计值●脉冲响应估计值ckRtaNNkgMzPP)()()(ˆ21●脉冲响应估计误差PPNkNkgkgkgkg120120)()(ˆ)(=四、实验步骤(1)掌握相关分析辨识方法的基本原理;(2)设计实验方案,画出程序框图;(3)编制实验程序;(4)调试并运行程序,记录数据;111TsK/u(k)x(k)211Ts/y(k)(5)分析实验结果,完成实验报告。五、实验报告格式参见附录一。实验3最小二乘法的实现一、实验目的理解并掌握系统辨识中的最小二乘法原理。二、实验原理给定系统12()(1)(2)()nykaykaykaykn01()(1)()()nbukbukbuknk(1)其中12,,,naaa,012,,,,nbbbb为待辨识的未知参数,()k是不相关随机序列。y为系统的输出,u为系统的输入。分别测出nN个输出、nN输入值(1),(2),(3),(),(1),(2),()yyyynNuuunN,则可写出N个方程,具体写成矩阵形式,有10(1)()(1)(1)(1)(1)(2)(1)(2)(2)(2)(2)()(1)()()()()nnaynynyununaynynyununbynNynNyNunNuNnNb(2)设10(1)(1)(2)(2),,()()nnaynnaynnybynNnNb,()(1)(1)(1)(1)(2)(2)(2)(1)()()()ynyunuynyunuynNyNunNuN则式(2)可写为y(3)式中:y为N维输出向量;为N维噪声向量;为21n维参数向量;为(21)Nn测量矩阵。为了尽量减小噪声对估值的影响,应取21Nn,即方程数目大于未知数数目。的最小二乘估计为1()TTy(4)三、实验内容对象的数学模型如下:()1.5(1)0.7(2)(1)0.5(2)()zkzkzkukukvk其中,)(kv是服从正态分布的白噪声N)1,0(。输入信号采用4阶M序列,幅度为1。选择如下形式的辨识模型:)()2()1()2()1()(2121kvkubkubkzakzakz设输入信号的取值是从k=1到k=16的M序列,则待辨识参数LSθˆ为LSθˆ=LτL1LτLzH)HH(。其中,被辨识参数LSθˆ、观测矩阵zL、HL的表达式为2121ˆbbaaLSθ,)16()4()3(zzzLz,)14()2()1()15()3()2()14()2()1()15()3()2(uuuuuuzzzzzzLH要求编制仿真程序,获取系统输入输出数据,并运用最小二乘法对这一系统的参数进行辨识,并将辨识结果与实际参数进行对比。四、实验步骤1.写出系统结构、实际参数、噪声源及输入信号等内容;2.画出程序框图(流程图);3.编制MATLAB的M文件;4.运行编制的M文件;5.查看程序运行结果并进行分析;6.填写实验报告。五、实验报告格式参见附录一。实验4递推最小二乘法的实现一、实验目的熟悉并掌握递推最小二乘法的算法原理。二、实验原理给定系统12()(1)(2)()nykaykaykaykn01()(1)()()nbukbukbuknk(1)其中12,,,naaa,012,,,,nbbbb为待辨识的未知参数,()k是不相关随机序列。y为系统的输出,u为系统的输入。分别测出nN个输出、nN输入值(1),(2),(3),(),(1),(2),()yyyynNuuunN,则可写出N个方程,具体写成矩阵形式,有10(1)()(1)(1)(1)(1)(2)(1)(2)(2)(2)(2)()(1)()()()()nnaynynyununaynynyununbynNynNyNunNuNnNb(2)设10(1)(1)(2)(2),,()()nnaynnaynnybynNnNb,()(1)(1)(1)(1)(2)(2)(2)(1)()()()ynyunuynyunuynNyNunNuN则式(2)可写为y(3)式中:y为N维输出向量;为N维噪声向量;为21n维参数向量;为(21)Nn测量矩阵。为了尽量减小噪声对估值的影响,应取21Nn,