系统建模与仿真-哈尔滨工业大学

《系统辨识》实验手册哈尔滨工业大学控制与仿真中心2018年5月目录实验1 白噪声和M序列的产生---------------------------------------------------------- 2 实验2 脉冲响应法的实现---------------------------------------------------------------- 5 实验3 递推最小二乘法的实现---------------------------------------------------------- 9 附录实验报告模板---------------------------------------------------------------------- 13实验1 白噪声、M 序列的产生一、实验目的1、熟悉并掌握产生均匀分布随机序列方法以及进而产生高斯白噪声方法2、熟悉并掌握M 序列生成原理及仿真生成方法二、实验原理1、混合同余法混合同余法是加同余法和乘同余法的混合形式，其迭代式如下：111(*)mod /n n n n x a x b MR x M +++=+⎧⎨=⎩ 式中a 为乘子，0x 为种子，b 为常数，M 为模。

混合同余法是一种递归算法，即先提供一个种子0x ，逐次递归即得到一个不超过模M 的整数数列。

2、正态分布随机数产生方法由独立同分布中心极限定理有：设随机变量12,,....,,...n X X X 相互独立，服从同一分布，且具有数学期望和方差：2(),()0,(1,2,...)k k E X D X k μσ==>=则随机变量之和1nk i X =∑的标准化变量:()nnnkk kXE X Xn Y μ--==∑∑∑近似服从(0,1)N 分布。

如果n X 服从[0, 1]均匀分布，则上式中0.5μ=，2112σ=。

即0.5nkXn Y -=∑近似服从(0,1)N 分布。

3、M 序列生成原理用移位寄存器产生M 序列的简化框图如下图所示。

该图表示一个由4个双稳态触发器顺序连接而成的4级移位寄存器，它带有一个反馈通道。

当移位脉冲来到时，每级触发器的状态移到下一级触发器中，而反馈通道按模2加法规则反馈到第一级的输入端。

三、实验内容1、生成均匀分布随机序列（1）利用混合同余法生成[0, 1]区间上符合均匀分布的随机序列，并计算该序列的均值和方差，与理论值进行对比分析。

要求序列长度为1200，推荐参数为a=65539，M=2147483647，0<x 0<M 。

（2）将[0, 1]区间分为不重叠的等长的10个子区间，绘制该随机序列落在每个子区间的频率曲线图，辅助验证该序列的均匀性。

（3）对上述随机序列进行独立性检验。

（该部分为选作内容）2、生成高斯白噪声利用上一步产生的均匀分布随机序列，令n=12，生成服从N(0,1)的白噪声，序列长度为100，并绘制曲线。

3、生成M 序列M 序列的循环周期取为63126=-=P N ，时钟节拍Sec 1=∆t ，幅度1=a ，逻辑“0”为a ，逻辑“1”为-a ，特征多项式65()F s s s =⊕。

生成M 序列的结构图如下所示。

要求编写Matlab程序生成该M序列，绘制该信号曲线，并分析验证M序列的性质。

四、实验步骤1．分别画出三部分实验内容的程序框图（流程图）；2．编制MATLAB的M文件；3．运行编制的M文件；4．查看程序运行结果并进行分析；5．填写实验报告。

五、实验报告格式参见附录一。

实验2 相关分析法辨识脉冲响应一、实验目的通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。

二、实验原理一个单入单出线性定常系统的动态特性可用它的脉冲响应函数g(σ)来描述。

这样，只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数g(τ)。

而在系统有正常输入的情形下，辨识脉冲响应的原理图如下图所示。

0 ()()()y t g x t d σσσ∞=-⎰则000()11lim ()()(){lim ()()}T TT T x t y t x t dt g x t x t dt d T T ττσστσ∞→∞→∞--=--⎰⎰⎰上式两端同乘，进而取时间均值，有0 ()()()xy x R g R d τστσσ∞=--⎰则这就是著名的维纳霍夫积分方程。

0() ()(), ()() ()()()()()()x x xy x xy x t R k R k R g R d kg R g kτδττσδτστστσστττ∞=-=--=-==⎰如果输入是，这时的自相关函数为则根据维纳霍夫积分方程可得或者 白噪声三、实验内容下图为本实验的原理框图。

系统的传递函数为)(s G ,其中S e c 26T S e c ,3812021..,===T K ；)()(k z k u 和分别为系统的输入和输出变量；)(k v 为测量白噪声，服从正态分布，均值为零，方差为2v σ，记作),(~)(20v N k v σ；)(k g 0为系统的脉冲响应理论值，)(ˆk g 为系统脉冲响应估计值，)(~k g 为系统脉冲响应估计误差。

系统的输入采用M 序列（采用实验1中的M 序列即可），输出受到白噪声)(k v 的污染。

根据过程的输入和输出数据{})(),(k z k u ，利用相关分析法计算出系统的脉冲响应值)(ˆk g ，并与系统的脉冲响应理论值)(k g 0比较，得到系统脉冲响应估计误差值)(~k g，当∞→k 时，应该有0→)(~k g 。

1、模拟过程传递函数)(s G ，获得过程的输入和输出数据})(),(k z k u （采样时间取1秒）。

(1) 惯性环节其中，T 为惯性环节的时间常数，K 为惯性环节的静态放大倍数。

若采样时间记作0T ，则惯性环节的输出可写成：k g =)(ˆ]2T t k /∆[]011111000T k u k u T e T TK k u e TK k y e k y TT T T T T )()())()()()()(///--+-+--+-=--- (2) 传递函数)(s G 仿真（串联） 21211111T s T s T T K s G //)(++=令211T T KK ＝，则)(s G 的表达框图为：2、互相关函数的计算∑++=-=PP N r N i PMz i z k i u rN k R )()()()(111其中，r 为周期数，1+=P N i 表示计算互相关函数所用的数据是从第二个周期开始的，目的是等过程仿真数据进入平稳状态。

（可分别令r =1、3，对比仿真结果） 3、c 的补偿补偿量c 应取)(1-P Mz N R －，不能取)(P Mz N R －。

因为)(k R Mz 是周期函数，则有)()(0Mz P Mz R N R ＝，故不能取)(P Mz N R －。

4、计算脉冲响应估计值● 脉冲响应估计值 []c k R ta N N k g Mz P P+∆+=)()()(ˆ21 ● 脉冲响应估计误差 ()∑∑==⎪⎭⎫ ⎝⎛-PPN k N k g k gk g k g 12120)()(ˆ)(＝δ四、实验步骤(1) 掌握相关分析辨识方法的基本原理； (2) 设计实验方案，画出程序框图； (3) 编制实验程序；(4) 调试并运行程序，记录数据；(5) 分析实验结果，完成实验报告。

五、实验报告格式参见附录一。

实验3 递推最小二乘法一、实验目的熟悉并掌握基本最小二乘法和递推最小二乘法的算法原理。

二、实验原理给定系统12()(1)(2)()n y k a y k a y k a y k n =-------+01()(1)()()n b u k b u k b u k n k ξ+-++-+ （1）其中12,,,n a a a ，012,,,,n b b b b 为待辨识的未知参数，()k ξ是不相关随机序列。

y 为系统的输出，u 为系统的输入。

分别测出n N +个输出、n N +输入值(1),(2),(3),(),(1),(2),(y y y y n N u u u n N ++，则可写出N 个方程，具体写成矩阵形式，有10(1)()(1)(1)(1)(1)(2)(1)(2)(2)(2)(2)()(1)()()()()n n a y n y n y u n u n a y n y n y u n u n b y n N y n N y N u n N u N n N b ξξξ⎡⎤⎢⎥+--++⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥+-+-++⎢⎥⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥+-+--++⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦（2） 设10(1)(1)(2)(2),,()()n n a y n n a y n n y b y n N n N b ξξθξξ⎡⎤⎢⎥++⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥++⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦，()(1)(1)(1)(1)(2)(2)(2)(1)()()()y n y u n u y n y u n u y n N y N u n N u N --+⎡⎤⎢⎥-+-+⎢⎥Φ=⎢⎥⎢⎥-+--+⎣⎦则式（2）可写为y θξ=Φ+ （3）式中：y 为N 维输出向量；ξ为N 维噪声向量；θ为21n +维参数向量；Φ为(21)N n ⨯+测量矩阵。

为了尽量减小噪声ξ对θ估值的影响，应取21N n >+，即方程数目大于未知数数目。

θ的最小二乘估计为1()T T y θ-=ΦΦΦ （4）为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。

设已获得的观测数据长度为N ，将式（3）中的y 、Φ和ξ分别用,,N N N Y ξΦ来代替，即N N N Y θξ=Φ+ （5）用N θ表示θ的最小二乘估计，则()1T TN N N N N Y θ-=ΦΦΦ （6）令()1T N N N P -=ΦΦ，则T N N N N P Y θ=Φ （7）如果再获得一组新的观测值(1)u n N ++和(1)y n N ++，则又增加一个方程111T N N N y ψθξ+++=+ （8）式中11(1),(1)N N y y n N n N ξξ++=++=++[]1()(1)(1)(1)TN y n N y N u n N u N ψ+=-+-++++将式（5）和式（8）合并，并写成分块矩阵形式，可得T 111N N N N N Y y ξθξψ+++⎡⎤Φ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦（9） 于是，类似地可得到新的参数估值1T T1T T T 1111N NN N N N N N N Y y θψψψ-+++++⎧⎫⎡⎤⎡⎤⎡⎤ΦΦΦ⎡⎤⎪⎪⎢⎥⎢⎥⎢⎥⎪⎪⎢⎥=⎨⎬⎢⎥⎢⎥⎢⎥⎢⎥⎪⎪⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦⎪⎪⎩⎭ ()T1T 11111N N N N N T N N N N N Y P y P Y y ψψ++++++⎡⎤Φ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦=Φ+ （10）式中()1T 1T T 111T11N NN N N T N N N N P ψψψψ-+++-++⎧⎫⎡⎤⎡⎤ΦΦ⎪⎪⎢⎥⎢⎥⎪⎪=⎨⎬⎢⎥⎢⎥⎪⎪⎢⎥⎢⎥⎣⎦⎣⎦⎪⎪⎩⎭=ΦΦ+ ()11T 11N N N P ψψ--++=+ （11）应用矩阵求逆引理，从求得1N P +与N P 的递推关系式出发，经过一系列的推导，最终可求得递推最小二乘法辨识公式：()T1111N N N N N N K y θθψθ++++=+- （12） ()1T 11111N N N N N N K P P ψψψ-++++=+ （13）()1TT111111N N N N N NN N N P P P P P ψψψψ-+++++=-+ （14）为了进行递推计算，需要给出N P 和N θ的初值0P 和0θ。