应用最小二乘一次完成法和递推最小二乘法算法的系统辨识

1、目录1 引言 21.1 概述 21.2 辨识的基本步骤 22 系统辨识输入信号的产生方法和理论依据 32.1 白噪声序列 32.1.1 白噪声序列的产生方法 32.2 M 序列的产生 42.2. 1 伪随机噪声 42.2.2 M 序列的产生方法 43 应用经典辨识方法的辨识方案。 63.1 经典辨识方法概述 63.2 经典辨识方法的实现 64 最小二乘法的理论基础 74.1 最小二乘法 74.1.1 最小二乘法估计中的输入信号 94.1.2 最小二乘估计的概率性质 94.2 递推最小二乘法 105 两种算法的实现方案 115.1 最小二乘法一次完成算法实现 115.1.1 最小二乘一次完成算法

2、程序框图 115.1.2 一次完成法程序 115.1.3 一次完成算法程序运行结果 115.1.4 辨识数据比较 125.1.5 程序运行曲线 125.2 递推最小二乘法的实现 125.2.1 递推算法实现步骤 125.2.2 程序编制思路： 135.2.3 递推最小二乘法程序框图 145.2.4 程序运行曲线 155.2.5 测试结果 165.2.6 递推数据表 166 结论 167 参考文献 178 附录 176应用最小二乘一次完成法和递推最小二乘法算法的系统辨识摘要 ：本题针对一个单输入单输出系统的便是问题，辨识的输入信号采用的是伪随机二位式序列（ M 序列），系统噪声为独立 同分布高斯

3、随机向量序列 （白噪声），辨识的算法是递推最小二乘法和广义最小二乘法， 本文简单描述应用经典辨识方法的辨识 方案，详细描述了输入信号、噪声的产生方法及 matlab 程序，阐述了用两种不同算法的辨识原理并对它们的推导过程及辨识程 序编制思路做了详细的描述。最后结合真值与估计值对不同辨识算法的优劣进行了比较。关键词 ：系统辨识 M 序列 最小二乘法1 引言1.1 概述系统辨识是现代控制理论中的一个分支， 它是根据系统的输入输出时间函数来确定描述系统行为的数学模型。 通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当 前可测量的系统的输入和输出预测系统输

4、出的未来演变，以及设计控制器。对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。对系统进行控制的主要问题 是根据系统的特性设计控制输入，使输出满足预先规定的要求。而系统辨识所研究的问题恰好是这些问题的 逆问题。通常，预先给定一个模型类卩= M（即给定一类已知结构的模型），一类输入信号u和等价准则J=L（y，yM）（ 般情况下，J是误差函数，是过程输出 y和模型输出yM的一个泛函）；然后选择使误差函数 J 达到最小的模型，作为辨识所要求的结果。系统辨识包括两个方面：结构辨识和参数估计。在实际的辨识过 程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是

5、可以交织在一起进行 的。1.2 辨识的基本步骤 先验知识和建模目的的依据。先验知识指关于系统运动规律、数据以及其他方面的已有知识。这些知识对 选择模型结构、设计实验和决定辨识方法等都有重要作用。用于不同目的的模型可能会有很大差别。 实验设计。辨识是从实验数据中提取有关系统信息的过程，设计实验的目标之一是要使所得到的数据能包 含系统更多的信息。主要包括输入信号设计，采样区间设计，预采样滤波器设计等。 结构辨识。即选择模型类中的数学模型M的具体表达形式。除线性系统的结构可通过输入输出数据进行辨识外 ，一般的模型结构主要通过先验知识获得。 参数估计。知道模型的结构后，用输入输出数据确定模型中的未知参

6、数。实际测量都是有误差的，所以参 数估计以统计方法为主。 模型适用性检验。造成模型不适用主要有三方面原因：模型结构选择不当；实验数据误差过大或数据代表性太差；辨识算法存在问题。检验方法主要有利用先验知识检验和利用数据检验两类。11.3设待辨识系统如图 1 所示。二阶系统为:参数真值为:1设 f (k)1待辨识系统图图a(z 1) b(z1) f(z1)1 21aza?z1 2b1zb2z1 2f1zf2zai1.5; a20.7;b|1；b20.50; fo1；f11；f20.2z1)1, (k)(k),(k)为白色噪声，简单描述应用经典辨识方法的辨识方案。12 (k)为有色噪声，f(z )f

7、1z 1f2z11.00; f20.2 ,(k)为独立同分布的高斯序列,0.4,2系统辨识输入信号的产生方法和理论依据2.1白噪声序列如果随机序列 w t 均值为0,并且两两不相的关的，对应自相关函数为Rw 12 l,l 0, 1, 2,ggg式中l0； 00则称这种随机序列为白噪声序列 2.1.1白噪声序列的产生方法F面主要介绍(0, 1)均匀分布和正态分布随机数的产生方法在计算机上产生(0, 1)均匀分布随机数的方法很多，其中最简单、最方便的是数学方法。产生伪随机 数的数学方法很多，其中最常用的是乘同余法和混合同余法。(1)乘同余法这种方法先用递推同余式产生正整数序列Xi=Axi-1(mo

8、dM),i=1,2,3式中：M为2的方幕，k为大于2的整数；A = 3(mod8)或A = 5(mod8),且A不能太小；初值 x0取正奇数，例如取 x0=1.再令 i 十 1,2,.则是伪随机序列，循环周期可达i(2)混合同余法混合同余法产生伪随机数的递推同余式为X Axi 1 c(mod M )式中 M 2k,K为大于2的整数；A三1(mod4),即A 2n 1，其中n为满足关系式2W nW 34的整数。初Xk值X0为非负整数。令 i 丄，则 i是循环周期为2的伪随机数序列M2.2 M序列的产生在进行系统辨识时，选用白噪声作为辨识输入信号可以保证获得较好的便是效果，但工程上难以实现。M序列

9、是一种很好的辨识输入信号，它具有近似白燥声的性质，不仅可以保证有较好的辨识效果，而且工程上 易于实现。M序列是伪随机二位式序列的一种形式。在介绍 M序列之前，先介绍一下伪随机噪声的概念。2.2. 1 伪随机噪声对白噪声的一个样本函数w(t)截取0,T时间内一段，对其它时间段T,2T,2T,3T,，以 周期T延拖下去，这样获得的函数 w(t)是周期T的函数，在0,T时间内是白噪声，在此时间之外是重复的白噪声，它的自相关 函数 RwE w t w t 的周期也是T。由于在0,T时间内自相关函数R,就是白噪声的自相关函数，它具有周期性，称为w(t)为伪随机噪声。2.2.2 M序列的产生方法M序列是一

10、种离散二位式：随机序列，所谓二位式”是指每个随机变量只有 2种状态。可用多级线性反馈 移位寄存器产生 M序列。M序列是最长线性反馈移存器序列的简称，是由带线性反馈的移存器产生的周期最长的一种序列。具有较强的抗干扰能力和较低的截获概率，而且长的M序列更容易在一定的强噪声中被提取，这样就能够充分保证数据的正常通信。通常产生伪随机序列的电路为反馈移存器.一般说来，由n级移位寄存器产生的周期为N=2?-1的M序列，在一个循环周期内，“0”出现的次数为 N 1/2, “1”出现的次数为 N 1/2。现在我们引入 M序列的本原多 项式的概念。若一个 n次多项式f (x)满足以下条件(1) f (X)为既约

11、的。(2) f (x)可整除(xm 1), m 2n 1。(3) f (x)除不尽(xq 1), q m则f (x)为本原多项式。一个4级M序列可以通过线性反馈移位寄存器产生，如下图所示：图2周期为15的伪随机序列产生器图每级移位寄存器由双稳态触发器和门电路组成，称为1位，分别以0和1来表示2种状态。当移位脉冲到来时，每位的内容移至下一位，最后1位移出的内容即为输出。为了保持连续工作，将最后2级寄存器的内容经过适当的逻辑运算后反馈到第1级寄存器作为输入。当一个移位脉冲到来后，第1级寄存器的内容送到第2级，第2级寄存器的内容送到第3级，第3级寄存器的内容送到第 4级，而第3级和第4级寄存器的内容

12、作模和2相加后再反馈到第 1级寄存器。产生伪随机序列时要求寄存器的初始状态不全为0,因为全0初始状态将导致各级寄存器输出永远是0。如果寄存器的初始内容都是1,第1个移位脉冲来到后，4级寄存器的内容变以 0111， 个周期的变化规律为：1111 0111 0011 0001 1000 0100 0010 1001 1100 01101011 0101 1010 1101 1110 1111 一个周期结束后，产生的15种不同的状态。计算机模拟产生 M序列非常方便，先定义输出序列长度和一个数组，数组个数等于移位寄存器的个数，通过 使用异或指令，再利用for循环指令，即可完成任意长度和级数M序列的产生

13、。图3M序列产生的流程图图4MATLAB 仿真M序列图形如图3是产生M级长度为L的M序列程序流程图，图 4所示是Matlab仿真4级移位寄存器，长度 15的M 序列输出图形。23应用经典辨识方法的辨识方案。3.1经典辨识方法概述用正弦信号、单位阶跃信号和M序列辨识系统的方法都称为经典辨识方法。在经典辨识方法中，用得最多的是脉冲响应。用 M序列作为输入信号，再用相关法处理测试结果，可以很方便地得到系统的脉冲响应。 用白噪声作为输入信号确定系统脉冲响应的方法，称为相关法。相关法，就是根据维纳-霍夫积分方程，禾U 用输入信号的自相关函数和输入与输出的互相关函数确定系统脉冲响应的方法。采用周期性的伪随

14、机信号作 为输入信号可以使计算变得简单，所以用M序列辨识系统脉冲响应是用得最广泛的一种方法。3.2经典辨识方法的实现由系统辨识图可得如下关系统：y(k)= x(k)+ (k) g(k) (k)(3.2.1)其中：（k）为输入白噪声信号，x（k）为实际输入信号，y（k）为实际输出信号，g（ k ）为输入信号为u（ k ）时的 脉冲响应。根据维纳-霍夫积分方程可得离散的维纳霍夫积分方程：N 1k)(322)g( k)Rx(k 0令式中Ruxg(k)Rx(k)(323)再设:RuxRux(0)Rux(1)Rux(2)Ru(0)Ru(1)Ru( 1)Ru(0)Ru( N 1)Ru( N 2)Rux(N

15、1)Ru(N 1)Ru(N 2)Ru(0)k则上式可写成以下形式:再根据上式（3.2.3）可得:RuxRgg 1R 1Rux(3.2.4)由二电平M序列的自相关函数的计算公式可以求得112 R aN1N11N1NR 1111NN根据递推公式求 m次观测的脉冲响应2 11N1 21(3.2.5)a N 11 12g(m)对系统进行m次观测，m由m次观测得到的Rux()用Rux( ,m )来表示，则有如下递推式:Rux(m,m)=m+ 1 k=oy(x) x(k- m)=Rux(m,m-1)+1m+ 1y(m) x(m- m) - Rux(m,m- 1)(326)综合以上公式就可以得到以下递推公式

16、：211x(m)1N121/ 、x(m- y(m)1)gm 1gm - gm-1 +.2m+ 1 a N+ 1 D112x(m-N+ 1)根据递推式，可从 gm 1及新的观测数据得到gm，随着观测数据的增加,(327)gm的精确度不断的提高。34最小二乘法的理论基础4.1最小二乘法设单输入单输出线性定长系统的差分方程表示为:y kak 1a2y k 2Lany k nb)u kb|Uk 1Lbnu knk其中nk n kk in(k)均值为0的白噪声，现分另iJ测出n+N个输出输入值y(1),y(2), ,y( n+N),u(1),u(2 ,iu(n+N),则可写出N个方程，写成向量一矩阵形式

17、y n 1y nLy 1u n 1Lu1y n 2y n 1Ly 2u n 2Lu2MMMMy n Nyn N 1 Ly N u n NLuNa1Mn 1(4.1.1)ann 2boMMn Nbn7yn 1Mn1yyn 2ann2M，b°JMyn NMnNbny nLy1un 1Lu 1y n 1Ly2un 2Lu 2MMMyn N 1LyNun NLu N则式(4.1.1)可写为y(4.1.2)式中：y为N维输出向量；E为N为维噪声向量；0为(2n+1)维参数向量；为N X(2n+1)测量矩阵。ai因此，式(4.1.1)是一个含有(2n+1)个未知参数，由 N个方程组成的联立方程组

18、。在给定输出向量y和测量矩阵的条件下求参数B的估计，这就是系统辨识问题。 设 ？表示的估计值，？表示y的最优估计，则有?式中：(4.1.3)? n? nMy nMb0Mbn设 e(k)=y(k)- ?(k), e(k)称为残差，则有e=y- ?=y-最小二乘估计要求残差的平方和最小，即按照指数函数Te eJ求j对 ？的偏导数并令其等于 0可得：(4.1.4)由式(4.1.5)可得的最小二乘估计:J为极小值的充分条件是:即矩阵T为正定矩阵，或者说是非奇异的。(4.1.5)(4.1.6)094.1.1最小二乘法估计中的输入信号当矩阵T的逆阵存在是，式(4.1.6)才有解。一般地，如果u(k)是随机

19、序列或伪随机二位式序列，则矩阵T是非奇异的，即( T)一1存在，式(4.1.6)有解。现在从T必须正定出发，讨论对u(k)的要求。n N 1Tyyyuk nuyuu(4.1.7)当N足够大时有1 TRW.P.1yRuyRNRyuRu(4.1.8)如果矩阵T正定，则Ru是是对称矩阵，并且各阶主子式的行列式为正。当N足够大时，矩阵 Ru才是是对称的。由此引出矩阵 T正定的必要条件是 u(k)为持续激励信号。如果序列u(k)的n+1阶方阵Ru是正定的，贝U称 u(k)的n+1阶持续激励信号。下列随机信号都能满足 Ru正定的要求1. 有色随机信号2. 伪随机信号3. 白噪声序列4.1.2最小二乘估计的

20、概率性质最小二乘估计的概率性质1) 无偏性由于输出值y是随机的，所以)是随机的，但注意B不是随机值。如果E ? E，则称？是 无偏估计2) 致性估计误差)的方差矩阵为2 彳1?1 TVarE(4.1.9)为有效值，参数(4.1.10)N N2limVar? lim R 10,w.p.1NN N上式表明当nis时，)以概率1趋近于e o当E (k)为不相关随机序列时，最小二乘估计具有无偏性和一致性3) 有效性如果式(4.1.2)中的E的均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值 估计偏差的方差达到 Cramer-Rao不等式的下界，即Var? 2E t 1 M 1式中M为Fishe

21、r矩阵，且114)渐近正态性In p y/M EIn p y /如果式(4.1.2)中的E是均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值2E T 14.2递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识设已获得的观测数据长度为N，则Yn(4.1.11) 服从正态分布，(4.1.2)(4.2.1)估计方差矩阵为式中?NVarPnPnNyn如果再获得1组新的观测值，则又增加 1个方程yN得新的参数估计值PN 1式中应用矩阵求逆引理，可得PN 1和PN的递推关系式2Pn|Yn1yN矩阵求逆引理设A为nx n矩阵，B和C为nx m矩阵,则有矩阵恒等式bct得到

22、递推关系式PnPn由于N PNN 1是标量，因而上式可以写成PN 1PnPn最后，得最小二乘法辨识公式?NPnPN 1PnPn1并且A + BCTI+CTA-1B1Bcta1B1cta 1KnPnyNN PNN PNn Pn1 Rl1 Fn1 N 1PNN 1 PN都是非奇异矩阵，(4.2.2)(4.2.3)1113(424)35有2种给出初值的办法(1) 设NO ( N0>n)为N的初始值，可算出PN0(2) 假定? O,P0 C2|,c是充分大的常数,PN 0noYN 0(425)I为(2n+1) x (2n+1)单位矩阵，则经过若干次递推之后能得到较好的参数估计。 5两种算法的实现

23、方案5.1最小二乘法一次完成算法实现如果把式(4.1.6)中的'取好足够的输入一输出数据以后一次计算出来，那么这种算法就式最小二乘法的一次完成法。5.1.1最小二乘一次完成算法程序框图图3.2最小二乘一次完成算法程序框图5.1.2 一次完成法程序具体程序参见附录45.1.3 一次完成算法程序运行结果ans =-1.50000.70001.00000.5000a1 =-1.50000.7000bl =1.00005.1.4辨识数据比较?at?t?1.50000.70001.0000.5000参数真值1.5000.7001.000.50误差00005.1.5程序运行曲线o-1o o oo

24、o OJK( 岀输 )岀输散离5.2递推最小二乘法的实现5.2.1递推算法实现步骤1）输入M序列的产生程序，可以参见3.2当中M序列产生方法（具体程序见附录。2）初始化。一种简单的方法是选择2（参照式Cl，其中C为尽量大的数，然后以它们为起始值进行计算4.2.3）。可以证明，当Cis时，按照递推公式算得的将与上面用批处理算法求得的结果相同，当N很大时，C的选择便无关紧要了。3）递推算法的停机标准:max?i k ?! k 1? k&是适当小的数。4）最小二乘估计的递推算法系统参数的步骤如下:u (i), y (i), i=1,2,n;1)产生系统输入信号 M 序列，输入系统阶次，计算输

25、入和输出数据2)置 N=n,?N0,PN 1010II 单位矩阵)3)形成行向量1 y(n)L y(Nn)u(N) Lu(N1n4)计算 K N 1PNN 1 PN5)输入新测量数据y ( N+1 )和u( N+1 )；6)计算?N?1NKN 1 y(N 1) N?N7)计算PN1IKN 1 N 1 PN ;8)输出?N1和 PN9)N N+1,转(3);5.2.2 程序编制思路 ：(1)产生M序列，通过计算，依据算出输出值，设定递推初始值，采用424方法2设定辨识参数初值，?为一个很小的量，P0为一个很大值的4X 4矩阵，设定误差量，作为停机标准。( 2)依据设定,计算 1y(2)y(1)

26、u(1) u(2) T , T1 P01 1可以算出为一标量,依据式425算出K1,K1为4X 1矩阵。( 3)依据式 4.2.5 计算出 ?1,和 P1 ,加入估计参数矩阵；( 4)计算误差大小,求出相对变化率,加入误差矩阵第一列。(5)再以?和R为已知参数，重复(2)-( 3)计算出参数 ?，并加入估计参数矩阵。6) 如此循环,计算出参数估计?n 。7) 判断误差变化率满足要求否,如果满足,则退出循坏,不再进行计算。8) 分离估计参数,绘出图形。523递推最小二乘法程序框图如下所示:图5最小二乘递推算法辨识的Malab程序流程图具体程序参见附录524程序运行曲线图6M序列信号的波形Para

27、meter Ide ntificati on with Recursive Least Squares Method-5000-l'l11t111-II-rrrr1IrIdentification Precision3002001000-100-200-300-4001020304050607080901005.2.5测试结果见附录65.2.6地退数据表递推次数a1召1误差a1 ?a2召2误;差a2?b误差b, bb2b2误;差d b1-1.500.0011.4990.70.0010.6991.00.0010.9990.50.0010.49992-1.501.50.700.71.001

28、.00.500.53-1.50.0011.4990.70.0010.6991.00.750.250.50.750.254-1.5-1.500.70.0010.6991.00.750.250.50.750.255-1.5-1.500.70.701.00.750.250.50.750.256-1.5-1.500.70.701.01.000.50.507-1.51.500.70.701.01.000.50.50从以上数据可以看出，随着递推次数的增加，参数估计的误差逐渐减少，最终趋于稳定。下面是辨识误差的分布趋势。（误差=新估计旧估计值古计值）56结论最小二乘法是应用广泛的辨识方法之一。可用于动态系统

29、，也可用于静态系统；可用于线性系统，也可用 于非线性系统。通过本课程设计，了解和掌握了基于最小二乘法原理的一次性完成算法和递推算法的实现方 法，完成了 Matlab 下的仿真计算，能够精确的估计辨识参数。最小二乘一次性完成算法是离线算法，需要采 集大量数据，一次性完成计算，因此，数据计算量大，当数据量很大时，数据输入不方便，但在本课程设计 过程当中，考虑到了此问题，运用相应的办法，解决了矩阵输入的问题。递推算法适合于在线算法，利用原 有参数估计进行下一步估计，可以做到运算量小，实时进行估计，根据仿真结果图示，可以看出计算一定量 的数据后，参数估计趋于稳定，只要满足误差要求，不必计算完所有数据。

30、两种算法不足之处，在考虑噪声 干扰时，看成了白噪声，具有不相关性。如果噪声引起的方程误差是相关的，而不是白噪声，可以通过下面 两种方法进行估计先估计未受噪音污染的系统输出，然后再估计模型的参数，就是辅助变量法；使用一 种滤波器，将相关的方程的误差转换为白噪声再进行估计，就是增广矩阵法。7 参考文献1 李言俊 张科,系统辨识理论及应用 ,国防工业出版社 ,2006.7 .2 刘宏才 ,系统辨识与参数估计 ,北京,冶金工业出版社 ,1999.3 黄道平，MATLAB与控制系统的数字仿真及 CAD化学工业出版社,2004.104 侯媛彬汪梅王立奇系统辨识及其MATLAB仿真科学出版社 2004.25

31、 蔡季冰，系统辨识，北京，北京理工大学出版社，1991.8 附录附录 1：用乘同余法产生白噪声 编程如下：A=6; x0=1; M=255; f=2;N=100 ；%初始化；x0=1; M=255;for k=1: N x2=A*x0;x1=mod （x2，M）; v1=x1/256;%乘同余法递推 100次；%分别用x2和x0表示xi+i和xi-1 ；%取x2存储器的数除以M的余数放x1(xi)中；%将x1存储器中的数除以256得到小于1的随机数放v1中；v（:，k）=（v1-0.5 ）*f;%将v1中的数(Q减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v(:,k)表示行

32、不变、列随递推循环次数变化；x0=x1;v0=v1;% xi-1= xi；end v2=v%递推100次结束；%该语句后无；'实现矩阵存储器 v中随机数放在v2中，且可直接显示在MATLAB 的 window 中；k1=k; %grapher k=1:k1;%以下是绘图程序；plot(k，v，k，v， 'r ');xlabel(k')， ylabel( v');tktle( (-1，+1)均匀分布的白噪声 ') 1.2 程序运行结果如图 6所示。图6采用MATLAB产生的（-1,+1）均匀分布的白噪声序列附录2产生的（-1 , 1）均匀分布的白噪

33、声序列在程序运行结束后，产生的（-1, 1）均匀分布的白噪声序列，直接从 MATLAB的window界面中copy出来如下（v2中每行存6个随机数）：v2 =-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.718

34、80.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-

35、0.9531-0.71880.6875-0.8359-0.01560.92190.57030.4531-0.2500-0.48440.1016-0.36720.8047-0.13280.21880.3359-0.9531-0.71880.6875-0.8359附录3 M序列产生程序function Out=Mge nerator(le ngth,ple ngth,a)%plength=4;length=15; %M序列周期=2yiength-1 ,% 置M序列总长度和级数，a是幅度for n=1:plength%移位寄存器输入 X(i)初T态(1111 )X( n)=1;endfor i=1:

36、le ngthfor j=ple ngth:-1:2X(j)=X(j-1);endX(1)=xor(X(ple ngth-1),X(ple ngth);%异或运算if X(ple ngth)=OOut(i)=-1;elseOut(i)=X(ple ngth);endend%save ('mdata','Out');%绘图i1=ik=1:1:le ngth;plot(k,Out,k,Out,'rx')xlabel('k')ylabel('M 序列')title('移位寄存器产生的M序列')end附录4

37、程序运行结果如下图所示移位寄存器产生的ri序列图7移位寄存器产生的M序列附录4最小二乘一次完成法%Oncefinish %model is y(k)-1.5y(k-1)+0.7(k-2)=u(k-1)+0.5(k-2)+v(k); L=50;A=1;n=2; %output number; u=mseries(A,L);%u=Mserise(L,M,A); %系统辨识的输入信号为一个周期的M序列，信号长度L ;幅度A;z=zeros(1,16);%定义输出观测值的长度 1*16dd=zeros(1,4); gg=zeros(14,1);for k=n+1:L+1z(k)=1.5*z(k-1)-

38、0.7*z(k-2)+u(k-1)+0.5*u(k-2); % 用理想输出值作为观测值end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画出输入信号 u 的经线图形 ylabel(' 输入 u(k)')subplot(3,1,2) i=1:1:L+1;plot(i,z)ylabel(' 输出 z(k)')subplot(3,1,3)%画三行一列图形窗口中的第三个图形stem(z),grid on%画出输出观测值 z 的经线图形，并显示坐标网格ylabel(' 离散输出 z(k)')%u,z;%显示输入信号

39、和输出观测信号%L=L-1% 数据长度for i=n+1:L+1h=-z(i-1) -z(i-2) u(i-1) u(i-2);dd(i-2,:)=h;%dd=(L-1)*2nend%HL=-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11

40、) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14) % 给样本矩阵 HL 赋值 for j=n+1:L+1zz=z(j);gg(j-2,:)=zz;end%c=inv(dd'*dd)*dd'*gg;给样本矩阵 zL 赋值%ZL=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)%ca

41、lculating parameters% 计算参数c1=dd'*dd; c2=inv(c1); c3=dd'*gg; c=c2*c3;c'%c1=dd'*dd; c2=inv(c1); c3=dd'*ZL; c=c2*c3 % 计算并显示 %DISPLAY PARAMETERS a1=c(1), a2=c(2), b1=c(3), b2=c(4); % 从 中分离出并显示 a1 、 a2、 b1、 b2 %End 附录 5 递推最小二乘算法程序 clear;clc;% 清理工作间变量L=100;% M 序列的周期 y1=1;y2=1;y3=1;y4=0

42、;% 四个移位积存器的输出初始值for i=1:L;% 开始循环，长度为 Lx1=xor(y3,y4);% 第一个移位积存器的输入是第 3 个与第 4 个移位积存器的输出的“或”3 个移位积存器的输出2 个移位积存器的输出3 个移位积存器的输出x2=y1;% 第二个移位积存器的输入是第x3=y2;% 第三个移位积存器的输入是第x4=y3;% 第四个移位积存器的输入是第y(i)=y4;%取出第四个移位积存器幅值为"0" 和 "1" 的输出信号if y(i)>0.5,u(i)=-0.5;%如果M序列的值为"1"时，辨识的输入信号取“

43、 -0.03else u(i)=0.5;% 当M序列的值为"0"时,辨识的输入信号取“ 0.03 ” end% 小循环结束y1=x1;y2=x2;y3=x3;y4=x4;% 为下一次的输入信号做准备en d%大循环结束，产生输入信号 u figure(1);% 第 1 个图形stem(u),grid on% 以径的形式显示出输入信号并给图形加上网格 z(2)=0;z(1)=0;% 取 z 的前两个初始值为零 for k=3:L;% 循环变量从 3 到 25z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%给出理想的辨识输出采样信号e

44、nd %RLS递推最小二乘辨识c0=0.001 0.001 0.001 0.001'% 直接给出被辨识参数的初始值 ,即一个充分小的实向量p0=10A6*eye(4,4);% 直接给出初始状态P0,即一个充分大的实数单位矩阵E=0.000000005;% 相对误差 E=0.000000005c=c0,zeros(4,99);% 被辨识参数矩阵的初始值及大小e=zeros(4,100);% 相对误差的初始值及大小for k=3:L; % 开始求 Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2)' x=h1'*p0*h1+1; x1=inv(x); %开始

45、求 K(k)k仁pO*h1*x1;%求出K的值d1=z(k)-h1'*c0; c1=c0+k1*d1;%求被辨识参数 ce1=c1-c0;% 求参数当前值与上一次的值的差值e2=e1./c0;% 求参数的相对变化e(:,k)=e2; % 把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;% 新获得的参数作为下一次递推的旧参数c(:,k)=c1;% 把辨识参数 c 列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1'*h1'*p0*h1+1;% 求出 p(k) 的值p0=p1;% 给下次用if e2<=E break;% 若参数收敛满足要求，终止计算end

46、% 小循环结束end%大循环结束c%显示被辨识参数e%显示辨识结果的收敛情况%分离参数a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); figure(2);% 第 2 个图形i=1:L;% 横坐标从 1 到 15 plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出 a1， a2， b1， b2 的各次辨识结果title('Parameter Identification with Recursive Least Squares Method')%图形标题figure(3); % 第 3 个图形 i=1:L; % 横坐标从 1 到 15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出 a1，a2，b1，b2 的各次辨识结果的收敛情况title('Identification Precision') %图形标题附录 6 递推最小二乘法程序测试结果c =Columns 1 t