1最小二乘法的理论基础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)()()()()()()()()1201121n n y k a y k a y k a y k n b u k b u k b u k n k ξ=-------++-++-+()()()()()()101122,,n n a y n n y n a n y b y n N n N b ξξθξξ⎡⎤⎢⎥++⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥++⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦()()()()()()()()()()()()()()()()()()10111121222112n n y n y n y u n u y n y n y u n u y n N y n N y N u n N u N a n a n b n N b ξξξ+--+⎡⎤⎡⎤⎢⎥⎢⎥+-+-+⎢⎥⎢⎥=⨯⎢⎥⎢⎥⎢⎥⎢⎥+-+--+⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎢⎥+⎡⎤⎢⎥⎢⎥⎢⎥+⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦则式(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ˆˆ,ˆˆˆn n ay n a y n y b y n N b θ⎡⎤⎢⎥+⎡⎤⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦设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θ=Φ()()ˆˆTT J e e y y θθ==-Φ-Φˆθ()ˆ20ˆT J y θθ∂=-Φ-Φ=∂ˆT T y θΦΦ=Φ()1ˆT T y θ-=ΦΦΦ220ˆTJ θ∂=ΦΦ>∂1n N yyyu T+-ΦΦ⎡⎤当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 是随机的，所以θ是随机的，但注意θ不是随机值。

如果{}{}ˆE E θθθ==，则称ˆθ是θ无偏估计 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且服从正态分布的白噪声向量，则最小二乘参数估计值ˆθ服从正态分布，即..11y uy W P T yu u R R R R R N ⎡⎤ΦΦ−−−→=⎢⎥⎣⎦121ˆT Var E NN σθ-⎧⎫⎪⎪⎛⎫=ΦΦ⎨⎬⎪⎝⎭⎪⎪⎩⎭21ˆlim lim 0,..1N N Var R w p Nσθ-→∞→∞==ˆθ(){}121ˆT Var E M θσ--=ΦΦ=()()ˆˆln /ln /ˆˆTp y p y M E θθθθ⎧⎫⎡⎤⎡⎤∂∂⎪⎪⎢⎥⎢⎥=⎨⎬⎢⎥⎢⎥∂∂⎪⎪⎣⎦⎣⎦⎩⎭12ˆT -1.2递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识 设已获得的观测数据长度为N ，则估计方差矩阵为 式中 于是如果再获得1组新的观测值，则又增加1个方程得新的参数估计值 式中应用矩阵求逆引理，可得1N P +和N P 的递推关系式矩阵求逆引理 设A 为n ×n 矩阵，B 和C 为n ×m 矩阵，并且A ，A ＋BCT 和I+CTA-1B 都是非奇异矩阵，则有矩阵恒等式（4.2.2）得到递推关系式由于1TN N N P ψψ+是标量，因而上式可以写成（4.2.3）最后，得最小二乘法辨识公式（4.2.4）有2种给出初值的办法（1）设N0（N0>n ）为N 的初始值，可算出（4.2.5）N N NY θξ=Φ+()1ˆNT T N N N NY θ-=ΦΦΦ()122T N NN Var P θσσ-=ΦΦ=()1T N N N P -=ΦΦˆNT N N N P Yθ=Φ111T N N N y ψθξ+++=+()1111ˆT N N N N N N P Y y θψ++++=Φ+()11111111TN N TT T N N N N N N P P ψψψψ---+++++⎧⎫ΦΦ⎡⎤⎡⎤⎪⎪==+⎨⎬⎢⎥⎢⎥⎣⎦⎣⎦⎪⎪⎩⎭()()111111T T T A BC A A B I C A B C A ------+=-+()11111T TN N N N N N N N NP P P I P P ψψψψ-++++=-+()111111T T N N N N N N N N NP P P P P ψψψψ-++++=-+()10000000ˆ,TT N N N N N N N P P Y θ-=ΦΦ=Φ()1111ˆˆˆT N N N N N NK y θθψθ++++=+-()11111T N N N N N N K P P ψψψ-+++=+()1111111T T N N N N N N N N NP P P P P ψψψψ-+++++=-+（4.2.1）（2）假定200ˆ0,,P c I θ==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.5000 0.7000 1.0000 0.5000 a1 = 1.5000 a2 = 0.7000 b1 =1.00002.1.4辨识数据比较2.1.5程序运行曲线-101输入u (k )0102030405060-505输出z (k)0102030405060-505离散输出z (k )2.2递推最小二乘法的实现2.2.1递推算法实现步骤1）输入M 序列的产生程序，可以参见3.2当中M 序列产生方法（具体程序见附录。

） 2）初始化。

一种简单的方法是选择ˆ0nθ=和2n P C I =，其中C 为尽量大的数，然后以它们为起始值进行计算(参照式2.2.3)。

可以证明，当C →∞时，按照递推公式算得的将与上面用批处理算法求得的结果相同，当N 很大时，C 的选择便无关紧要了。

3）递推算法的停机标准：()()()1ˆˆ1maxˆi i iik k k θθεθ-∀--<，ε是适当小的数。

4）最小二乘估计的递推算法系统参数的步骤如下：（1）产生系统输入信号M 序列，输入系统阶次，计算输入和输出数据u （i ），y （i ），i=1,2,…,n;（2）置N=n ，10ˆ0,10N NP I θ==（I …单位矩阵）； （3）形成行向量[]1()(1)()(1N y n y N n u N u N n ψ+=--+-+-（4）计算()()()111111;TTN N N K P N N P N ψψψ-+⎡⎤=++++⎣⎦（5）输入新测量数据y ( N+1 )和u （N+1）；（6）计算()11ˆˆˆ(1)1;N N N NK y N N θθψθ++⎡⎤=++-+⎣⎦（7）计算[]111;N N N N P I K P ψ+++=- （8）输出1ˆN θ+和1N P +； （9）N ←N+1,转（3）； 2.2.2程序编制思路：（1）产生M 序列，通过计算，依据算出输出值，设定递推初始值，采用4.2.4方法2设定辨识参数初值，0ˆθ为一个很小的量，P0为一个很大值的4×4矩阵，设定误差量，作为停机标准。

（2）依据设定，计算[]1(2)(1)(1)(2)Ty y u u ψ=--，1011T P ψψ⨯⨯+可以算出为一标量，依据式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递推最小二乘法程序框图如下所示：具体程序参见附录2.2.4程序运行曲线-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5图6M 序列信号的波形102030405060708090100-0.4-0.200.20.40.60.811.21.41.6Parameter Identification with Recursive Least Squares Method0102030405060708090100-400-2000200400600800100012001400Identification Precision2.2.5测试结果 见附录62.2.6地退数据表差的分布趋势。