最小二乘法

第3章 线性动态模型参数辨识－最小二乘法3.1 辨识方法分类根据不同的辨识原理，参数模型辨识方法可归纳成三类： ① 最小二乘类参数辨识方法，其基本思想是通过极小化如下准则函数来估计模型参数：min )()ˆ(ˆ==∑=θθLk k J 12ε 其中)(k ε代表模型输出与系统输出的偏差。

典型的方法有最小二乘法、增广最小二乘法、辅助变量法、广义最小二乘法等。

② 梯度校正参数辨识方法，其基本思想是沿着准则函数负梯度方向逐步修正模型参数，使准则函数达到最小，如随机逼近法。

③ 概率密度逼近参数辨识方法，其基本思想是使输出z 的条件概率密度)|(θz p 最大限度地逼近条件0θ下的概率密度)|(0θz p ，即)|()ˆ|(0m a x θθz p z p −−→−。

典型的方法是极大似然法。

3.2 最小二乘法的基本概念● 两种算法形式 ① 批处理算法：利用一批观测数据，一次计算或经反复迭代，以获得模型参数的估计值。

② 递推算法：在上次模型参数估计值)(ˆ1-k θ的基础上，根据当前获得的数据提出修正，进而获得本次模型参数估计值)(ˆk θ，广泛采用的递推算法形式为() ()()()~()θθk k k k d z k =-+-1K h其中)(ˆk θ表示k 时刻的模型参数估计值，K (k )为算法的增益，h (k -d ) 是由观测数据组成的输入数据向量，d 为整数，)(~k z 表示新息。

● 最小二乘原理定义：设一个随机序列)},,,(),({L k k z 21∈的均值是参数θ 的线性函数E{()}()T z k k θ=h其中h (k )是可测的数据向量，那么利用随机序列的一个实现，使准则函数21()[()()]LT k J z k k θθ==-∑h达到极小的参数估计值θˆ称作θ的最小二乘估计。

● 最小二乘原理表明，未知参数估计问题，就是求参数估计值θˆ，使序列的估计值尽可能地接近实际序列，两者的接近程度用实际序列与序列估计值之差的平方和来度量。

● 如果系统的输入输出关系可以描述成如下的最小二乘格式()()()T z k k e k θ=+h 式中z (k )为模型输出变量，h (k )为输入数据向量，θ为模型参数向量，e (k )为零均值随机噪声。

为了求此模型的参数估计值，可以利用上述最小二乘原理。

根据观测到的已知数据序列)}({k z 和)}({k h ，极小化下列准则函数21()[()()]LT k J z k k θθ==-∑h即可求得模型参数的最小二乘估计值θˆ。

● 最小二乘估计值应在观测值与估计值之累次误差的平方和达到最小值处，所得到的模型输出能最好地逼近实际系统的输出。

3.3 最小二乘问题的描述 (1) 考虑模型)()()()()(11k e k u z B k z z A +=--式中u (k )和z (k ) 分别为过程的输入和输出变量，e (k )是均值为零、方差为2nσ的随机噪声，)(1-z A 和)(1-z B 为迟延算子多项式，写成A z a z a z a zB z b z b z b z n nn n a a b b ()()--------=++++=+++⎧⎨⎪⎩⎪11122111221(2) 假定模型阶次n a 和n b 为已知，且有b a n n ≥，也可设n n n b a ==,并定义1212()[(1),,(),(1),,()][,,,,,,,]a b Ta b Tn n k z k z k n u k u k n a a a b b b θ⎧=------⎪⎨=⎪⎩h (3) 将模型写成最小二乘格式()()()T z k k e k θ=+h对于L k ,,, 21= (L 为数据长度)，可以构成如下线性方程组L L L θ=+e z H 式中T T T [(1),(2),,()](0)(1)(0)(1)(1)(1)(2)(1)(2)(2)(1)()(1)()()[(1),(2),,()]T L a b a b L a b TLz z z L z z n u u n z z n u u n z L z L n u L u L n L n n n L ⎧=⎪----⎡⎤⎡⎤⎪⎢⎥⎢⎥⎪----⎪⎢⎥⎢⎥==⎨⎢⎥⎢⎥⎪⎢⎥⎢⎥⎪------⎣⎦⎣⎦⎪=⎪⎩z h h H h e(4) 噪声的统计性质{}{}{}E 0,cov E L T L L L e ===∑e e e e (5) 噪声与输入不相关{}E ()()0,,e k e k l k l -=∀ (6) 数据长度L 充分大3.4 最小二乘问题的解 ● 考虑模型11()()()()()A z z k B z u k e k --=+● 准则函数取2()()[()()]LT k 1J k z k k θθ==Λ-∑h其中)(k Λ为加权因子，对所有的k ，)(k Λ都必须大于零。

● 准则函数又可写成()()()T L L L L L J θθθ=-Λ-z H z H式中L Λ为加权矩阵，它是正定的对角阵，由加权因子)(k Λ构成ΛL =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡)()()(L ΛΛΛ 0020001 ● 该准则函数)(θJ 可用以衡量模型输出与实际系统输出的接近情况．极小化这个准则函数，即可求得模型的参数估计值，使模型的输出能最好地预报系统的输出。

● 当TL L L ΛH H 是正则矩阵时，模型的加权最小二乘解为1WLSˆ()T T LLLLL Lθ-=ΛΛH H H z ● 通过极小化准则函数)(θJ 求得模型参数估计值WLS θˆ的方法称作加权最小二乘法，记作WLS (Weighted Least Squares algorithm)，对应的WLS θˆ称为加权最小二乘估计值。

● 如果加权矩阵取单位阵，即I =L Λ，则加权最小二乘解退化成普通最小二乘解1LS ˆ()T T L L L Lθ-=H H H z 这时的LS θˆ称之为最小二乘估计值，对应的估计方法称作最小二乘法，记作LS(Least Squares algorithm)。

最小二乘法是加权最小二乘法的一种特例。

3.5 最小二乘估计的可辨识性 ● 基于参数模型的辨识问题实际上可以归结为模型参数的最优化问题。

当给定输入输出数据时，对假定的模型结构是否能唯一地确定模型的参数，这就是可辨识问题。

● 可控性、可观性和可辨识性之间的关系● 可辨识性和输入信号的关系：过程的所有模态都必须被输入信号持续激励。

● 常用的输入信号：随机序列、伪随机序列、离散序列● 最小二乘估计的可辨识条件为矩阵L L L H H Λτ必须是非奇异的。

● 这一要求与数据集是“提供信息”的或辨识输入信号是“持续激励”的概念密切相关。

例3.3 考虑仿真对象)()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z +-+-=-+-- (3.41) 其中，)(k v 是服从正态分布的白噪声N )1,0(。

输入信号采用4阶M 序列，幅度为1。

选择如下形式的辨识模型)()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+ (3.42)设输入信号的取值是从k =1到k =16的M 序列，则待辨识参数LSθˆ为LS θˆ=L τL 1L τL z H )H H -(。

其中，被辨识参数LSθˆ、观测矩阵z L 、H L 的表达式为 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=2121ˆb b a a LSθ , ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)16()4()3(z z z L z ， ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------=)14()2()1()15()3()2()14()2()1()15()3()2(u u u u u u z z z z z z L H (3.43) 程序框图如图3.2所示。

Matlab6.0仿真程序如下：%二阶系统的最小二乘一次完成算法辨识程序u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列z=zeros(1,16); %定义输出观测值的长度for k=3:16z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值endsubplot(3,1,1) %画三行一列图形窗口中的第一个图形stem(u) %画输入信号u的径线图形subplot(3,1,2) %画三行一列图形窗口中的第二个图形i=1:1:16; %横坐标范围是1到16，步长为1plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线subplot(3,1,3) %画三行一列图形窗口中的第三个图形stem(z),grid on %画出输出观测值z的径线图形，并显示坐标网格u,z %显示输入信号和输出观测信号%L=14 %数据长度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) -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)] %给样本矩阵H L赋值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)] % 给样本矩阵z L赋值%Calculating Parametersc1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c=c2*c3 %计算并显示θˆLS%Display Parametersa1=c(1), a2=c(2), b1=c(3),b2=c(4) %从θˆ中分离出并显示a1、a2、b1、b2LS%End程序运行结果：>>u =[ -1，1，-1，1，1，1，1，-1，-1，-1，1，-1，-1，1，1]z =[ 0，0，0.5000，0.2500，0.5250，2.1125，4.3012，6.4731，6.1988，3.2670，-0.9386，-3.1949，-4.6352，6.2165，-5.5800，-2.5185]HL =0 0 1.0000 -1.0000-0.5000 0 -1.0000 1.0000-0.2500 -0.5000 1.0000 -1.0000-0.5250 -0.2500 1.0000 1.0000-2.1125 -0.5250 1.0000 1.0000-4.3012 -2.1125 1.0000 1.0000-6.4731 -4.3012 -1.0000 1.0000-6.1988 -6.4731 -1.0000 -1.0000-3.2670 -6.1988 -1.0000 -1.00000.9386 -3.2670 1.0000 -1.00003.1949 0.9386 -1.0000 1.00004.6352 3.1949 -1.0000 -1.00006.2165 4.6352 1.0000 -1.00005.58006.2165 1.0000 1.0000ZL =[ 0.5000，0.2500，0.5250，2.1125，4.3012，6.4731，6.1988，3.2670，-0.9386，-3.1949，-4.6352，-6.2165，-5.5800，-2.5185]Tc =[ -1.5000，0.7000，1.0000，0.5000]Ta1 = -1.5000a2 = 0.7000b1 = 1.0000b2 =0.5000>>-101-10010-10010何噪声成分，所以辨识结果也无任何误差。