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

11最小二乘法的理论基础1.1最小二乘法设单输入单输出线性定长系统的差分方程表示为：其中δ(k)为服从N(0,1)的随机噪声,现分别测出n+N个输出输入值y(1),y(2),…,y(n+N),u(1),u(2),…,u(n+N),则可写出N个方程，写成向量－矩阵形式(4.1.1)1201121nnykaykaykayknbukbukbuknk101122,,11112221nnaynnynanybynNnNbynyunuynyunuynNyNunNuN10111121222112nnynynyunuynynyunuynNynNyNunNuNananbnNb2则式(1.1.1)可写为（4.1.2）式中：y为N维输出向量；ξ为N为维噪声向量；θ为（2n+1）维参数向量；Φ为N×（2n+1）测量矩阵。因此，式（4.1.1）是一个含有（2n+1）个未知参数，由N个方程组成的联立方程组。11y在给定输出向量y和测量矩阵Φ的条件下求参数θ的估计，这就是系统辨识问题。设表示的估计值，ŷ表示y的最优估计，则有（4.1.3）式中：10ˆˆ1ˆˆ2ˆˆ,ˆˆˆnnaynaynybynNb设e(k)=y(k)-ŷ(k),e(k)称为残差，则有e=y-ŷ=y-Φθ最小二乘估计要求残差的平方和最小，即按照指数函数（4.1.4）求Ｊ对的偏导数并令其等于0可得：（4.1.5）由式（4.1.5）可得的最小二乘估计：（4.1.6）J为极小值的充分条件是：即矩阵ΦTΦ为正定矩阵，或者说是非奇异的。1.1.1最小二乘法估计中的输入信号当矩阵ΦTΦ的逆阵存在是，式（1.1.6）才有解。一般地，如果u(k)是随机序列或伪随机二位式序列，则矩阵ΦTΦ是非奇异的，即（ΦTΦ）－1存在，式（1.1.6）有解。现在从ΦTΦ必须正定出发，讨论对u(k)的要求。yˆˆˆyˆˆTTJeeyyˆˆ20ˆTJyˆTTy1ˆTTy220ˆTJ1nNyyyuTknuyuu3（4.1.7）当N足够大时有（4.1.8）如果矩阵ΦTΦ正定，则Ru是是对称矩阵，并且各阶主子式的行列式为正。当N足够大时，矩阵Ru才是是对称的。由此引出矩阵ΦTΦ正定的必要条件是u(k)为持续激励信号。如果序列{u(k)}的n+1阶方阵Ru是正定的，则称{u(k)}的n+1阶持续激励信号。下列随机信号都能满足Ru正定的要求1.有色随机信号2.伪随机信号3.白噪声序列1.1.2最小二乘估计的概率性质最小二乘估计的概率性质1）无偏性由于输出值y是随机的，所以是随机的，但注意θ不是随机值。如果ˆEE，则称ˆ是无偏估计2）一致性估计误差的方差矩阵为（4.1.9）上式表明当N→∞时，以概率1趋近于θ。当ξ(k)为不相关随机序列时，最小二乘估计具有无偏性和一致性3)有效性如果式（1.1.2）中的ξ的均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值为有效值，参数估计偏差的方差达到Cramer-Rao不等式的下界，即（4.1.10）式中M为Fisher矩阵，且（4.1.11）4）渐近正态性如果式（4.1.2）中的ξ是均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值ˆ服从正态分布，即..11yuyWPTyuuRRRRRN121ˆTVarENN21ˆlimlim0,..1NNVarRwpNˆ121ˆTVarEMˆˆln/ln/ˆˆTpypyME12ˆ,TNE4（4.1.2）1.2递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识设已获得的观测数据长度为N，则估计方差矩阵为式中于是如果再获得1组新的观测值，则又增加1个方程得新的参数估计值式中应用矩阵求逆引理，可得1NP和NP的递推关系式矩阵求逆引理设A为n×n矩阵，B和C为n×m矩阵，并且A，A＋BCT和I+CTA-1B都是非奇异矩阵，则有矩阵恒等式（4.2.2）得到递推关系式由于1TNNNP是标量，因而上式可以写成（4.2.3）最后，得最小二乘法辨识公式（4.2.4）有2种给出初值的办法（1）设N0（N0n）为N的初始值，可算出（4.2.5）NNNY1ˆNTTNNNNY122TNNNVarP1TNNNPˆNTNNNPY111TNNNy1111ˆTNNNNNNPYy11111111TNNTTTNNNNNNPP111111TTTABCAABICABCA11111TTNNNNNNNNNPPPIPP111111TTNNNNNNNNNPPPPP10000000ˆ,TTNNNNNNNPPY1111ˆˆˆTNNNNNNKy11111TNNNNNNKPP1111111TTNNNNNNNNNPPPPP（4.2.1）5（2）假定200ˆ0,,PcIC是充分大的常数，I为(2n+1)×(2n+1)单位矩阵，则经过若干次递推之后能得到较好的参数估计。[1][2][4]2两种算法的实现方案2.1最小二乘法一次完成算法实现如果把式（1.1.6）中的ˆ取好足够的输入—输出数据以后一次计算出来，那么这种算法就式最小二乘法的一次完成法。2.1.1最小二乘一次完成算法程序框图2.1.2一次完成法程序具体程序参见附录42.1.3一次完成算法程序运行结果ans=1.50000.70001.00000.5000a1=1.5000a2=0.7000b1=1.0000赋输入信号初值u定义输出观测值的长度并计算系统的输出值画出输入和输出观测值的图形给样本矩阵HL和zL赋值计算参数ˆ从ˆ中分离出并显示出被辨识参数a1,a2,b1,b2停机图3.2最小二乘一次完成算法程序框图62.1.4辨识数据比较ˆ1ˆa2ˆa1ˆb2ˆb1.50000.70001.0000.5000参数真值1.5000.7001.000.50误差00002.1.5程序运行曲线05101520253035404550-101输入u(k)0102030405060-505输出z(k)0102030405060-505离散输出z(k)2.2递推最小二乘法的实现2.2.1递推算法实现步骤1）输入M序列的产生程序，可以参见3.2当中M序列产生方法（具体程序见附录。）2）初始化。一种简单的方法是选择ˆ0n和2nPCI，其中C为尽量大的数，然后以它们为起始值进行计算(参照式2.2.3)。可以证明，当C→∞时，按照递推公式算得的将与上面用批处理算法求得的结果相同，当N很大时，C的选择便无关紧要了。3）递推算法的停机标准：1ˆˆ1maxˆiiiikkk，ε是适当小的数。74）最小二乘估计的递推算法系统参数的步骤如下：（1）产生系统输入信号M序列，输入系统阶次，计算输入和输出数据u（i），y（i），i=1,2,…,n;（2）置N=n，10ˆ0,10NNPI（I…单位矩阵）；（3）形成行向量1()(1)()(1NynyNnuNuNn（4）计算111111;TTNNNKPNNPN（5）输入新测量数据y(N+1)和u（N+1）；（6）计算11ˆˆˆ(1)1;NNNNKyNN（7）计算111;NNNNPIKP（8）输出1ˆN和1NP；（9）N←N+1,转（3）；2.2.2程序编制思路：（1）产生M序列，通过计算，依据算出输出值，设定递推初始值，采用4.2.4方法2设定辨识参数初值，0ˆ为一个很小的量，P0为一个很大值的4×4矩阵，设定误差量，作为停机标准。（2）依据设定，计算1(2)(1)(1)(2)Tyyuu，1011TP可以算出为一标量，依据式4.2.5算出K1,K1为4×1矩阵。（3）依据式4.2.5计算出1ˆ，和1P，加入估计参数矩阵；（4）计算误差大小，求出相对变化率，加入误差矩阵第一列。（5）再以1ˆ和1P为已知参数，重复（2）－（3）计算出参数2ˆ，并加入估计参数矩阵。（6）如此循环，计算出参数估计ˆn。（7）判断误差变化率满足要求否，如果满足，则退出循坏，不再进行计算。（8）分离估计参数，绘出图形。2.2.3递推最小二乘法程序框图如下所示：具体程序参见附录8Y工作间清零产生输出采样信号给被辨识参数和P赋初值按照式(4.2.4)的第三式计算P(k)计算被辨识参数的相对变化量参数收敛满足要求？停机按照式(4.2.4)的第二式计算K(k)按照式(4.2.4)的第一式计算(k)第四个移位寄存器的输出取反，并将幅值变为0.03得到辨识系统的输入信号样本值给M序列的长度L和移位寄存器的输入赋初始值画出被辨识参数的各次递推估计值图形分离参数画出被辨识参数的相对误差的图形画出辨识的输入信号径线图形图5最小二乘递推算法辨识的Malab程序流程图92.2.4程序运行曲线0102030405060708090100-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5图6M序列信号的波形0102030405060708090100-0.4-0.200.20.40.60.811.21.41.6ParameterIdentificationwithRecursiveLeastSquaresMethod100102030405060708090100-400-2000200400600800100012001400IdentificationPrecision2.2.5测试结果见附录62.2.6地退数据表递推次数1a1ˆa误差11ˆaa2a2ˆa误差22ˆaa1b2ˆb误差11ˆbb2b2ˆb误差22ˆbb11.500.0011.4990.70.0010.6991.00.0010.9990.50.0010.499921.501.50.700.71.001.00.500.531.50.0011.4990.70.0010.6991.00.750.250.50.750.2541.5-1.500.70.0010.6991.00.750.250.50.750.2551.5-1.500.70.701.00.750.250.50.750.2561.5-1.500.70.701.01.000.50.5071.51.500.70.701.01.000.50.50从以上数据可以看出，随着递推次数的增加，参数估计的误差逐渐减少