最小二乘参数辨识标准算法——第三讲

Harbin Institute of Technology– HIT系统辨识与自适应控制黄显林、班晓军 控制理论与制导技术研究中心 哈尔滨工业大学 banxiaojun@2010-3-15控制理论与制导技术研究中心第1页Harbin Institute of Technology– HIT第四讲 最小二乘参数辨识标准算法内容提要： 1. 最小二乘数学方法引例； 2. 最小二乘辨识方法的基本计算公式； 3. 算法演示与仿真分析； 4. 加权最小二乘法介绍。

2010-3-15控制理论与制导技术研究中心第2页Harbin Institute of Technology– HIT最小二乘方法的典故：1801年左右,德国数学家Gauss，在“星体轨道估计中”就发明了最小二 乘方法。

Gauss, K. F. (1809), Theoria Motus Corporum Celestium, English Translation: Theory of the Motion of Heavenly bodies. Dover(1963), New York.当时的天文界正在为火星和木星间庞大的间隙烦恼不已，认为火星和 木星间应该还有行星未被发现。

在1801年，意大利的天文学家 Piazzi， 发现在火星和木星间有一颗新星。

它被命名为「谷神星」(Cere)。

现在 我们知道它是火星和木星的小行星带中的一个，但当时天文学界争论 不休，有人说这是行星，有人说这是彗星。

必须继续观察才能判决， 但是 Piazzi只能观察到它 9 度的轨道，再来，它便隐身到太阳後面去 了。

因此无法知道它的轨道，也无法判定它是行星或彗星。

2010-3-15控制理论与制导技术研究中心第3页Harbin Institute of Technology– HIT高斯这时对这个问是产生兴趣，他决定解决这个捉摸不到的星体轨 迹的问题。

高斯自己独创了只要三次观察，就可以来计算星球轨道 的方法。

他可以极准确地预测行星的位置。

果然，谷神星准确无误 的在高斯预测的地方出现。

这个方法－－虽然他当时没有公布－－ 就是「最小二乘法」 (Method of Least Square)。

1802年，他又准确预测了小行星二号－－智神星 (Pallas) 的位置，这 时他的声名远播，荣誉滚滚而来，俄国圣彼得堡科学院选他为会 员，发现 Pallas 的天文学家 Olbers 请他当哥廷根天文台主任，他没 有立刻答应，到了1807年才前往哥廷根就任。

1809年他写了《天体运动理论》二册，第一册包含了微分方程、圆 椎截痕和椭圆轨道，第二册他展示了如何估计行星的轨道。

2010-3-15控制理论与制导技术研究中心第4页Harbin Institute of Technology– HIT一、引例 (参考《高等数学》 第四版 下册 同济大学数学教研室 高等教育出版社 PP. 78)问题：为了测定刀具的磨损速度，我们做这样的实验： 经过一定时间（每隔一个小时），测量一次刀具的厚 度，得到一组实验数据如下：2010-3-15控制理论与制导技术研究中心第5页Harbin Institute of Technology– HIT2010-3-15控制理论与制导技术研究中心第6页Harbin Institute of Technology– HIT27 26.5 26 25.5 25 24.5 0y (mm)123 4 Time (Hours)567图1. 引例中的数据图示 2010-3-15 控制理论与制导技术研究中心 第7页Harbin Institute of Technology– HIT7 7 7 ⎧ ⎧140a + 28b = 717 2 ⎨ ∑ t a =7 (∑ ⎪(⎩ 28ai +)8b + 208.5ti )b = ∑ ti yi i ,f ( = M⎪(a= 0⎧a t=−0.30360+−7 f (ti i)]02 = [i y 指标： 目的： ⎨ ⇒ b) =) ∑at= i b 结果： 推导过程： 0参见板书 7 ⎨ i= ⎪(∑ tb)= 27.125 = ∑ yi ⎩ i a + 8b ⎪ ⇒ 0y = f (t ) = −0.3036t + 27.125 i =0 ⎩ i=2010-3-15控制理论与制导技术研究中心第8页Harbin Institute of Technology– HIT27.5 27 26.5 y (mm) 26 25.5 25 24.5 0123 4 Time (Hours)567图2. 应用最小二乘法拟合的结果 2010-3-15 控制理论与制导技术研究中心 第9页Harbin Institute of Technology– HIT注：是以“偏差的平方和最小”为指标，在“线性 函数”这一类函数中找到的最优解。

如果函数类 扩大到“非线性”方程，有可能找到一个使“偏差 的平方和”更小的方程。

2010-3-15控制理论与制导技术研究中心第10页Harbin Institute of Technology– HIT二、基本定义和基本计算公式1. 最小二乘模型类和标准格式y (k ) + a1 y (k − 1) + L + an y (k − n) = b1u (k − 1) + L + bnu (k − n) + e(k )其中，{u ( k ), y ( k )} 为实测输入输出序列；{e(k )} 为白噪声序列；n为模型的阶数。

2010-3-15 控制理论与制导技术研究中心 第11页Harbin Institute of Technology– HIT2. 最小二乘准则(最小二乘估计原理)和正规方程组准则：J = eT ( N ) e ( N )正规方程组： ΦT ( N ) y ( N ) − Φ T ( N )Φ ( N )θ ( N ) = 02010-3-15控制理论与制导技术研究中心第12页Harbin Institute of Technology– HIT3. 参数解的表达式θˆ( N ) = [Φ T ( N )Φ ( N )]−1 Φ T ( N ) y ( N )4. 对 [ΦT (N )Φ(N )] 非奇异的进一步说明 结论：从矩阵理论上讲，[ΦT (N )Φ(N )] 非奇异要求 Φ(N ) 的行数（N）至 少要大于等于 N ≥ 2n (若待辨识参数为 2n 时） ；从物理背景上看N > 2n 。

问题： N > 2n ， [ΦT (N )Φ(N )] 一定就非奇异了吗？ 答：不一定。

2010-3-15 控制理论与制导技术研究中心 第13页Harbin Institute of Technology– HIT5. 开环可辨识性条件：2010-3-15控制理论与制导技术研究中心第14页Harbin Institute of Technology– HIT⎧U L = [ FuL , F 2uL ,L , F 2 nu L ] ⎪ uL = [u (1), u (2),L , u ( L)]T ⎪ ⎪ ⎛0 ⎪ 0⎞ ⎜ ⎟ ⎨ F =⎜1 O ⎟ ⎪ ⎜ ⎪ 0 1 0 ⎟ L× L ⎝ ⎠ ⎪ ⎪ n = max(na , nb ) ⎩参见：Astrom, K. J. and T. Bohlin, 1965, Numerical Identification of Linear Dynamic Systems from Normal Operating Records, IFAC Symposium On Theory of self-adaptive systems, Teddington, England.2010-3-15控制理论与制导技术研究中心第15页Harbin Institute of Technology– HIT开环可辨识性条件对实践的指导： 辨识信号不能任意选择，否则会造成不可辨识。

常用的 信号有： (1). 白噪声序列； (2). M序列或逆M序列； (3). n种频率的正弦信号的组合成的信号，但是各个 频率不能满足整数倍关系。

注：在推导最小二乘法的结果时，并不需要考虑噪声的 统计特性。

但是在评价最小二乘估计的性质时，需要用 到噪声的统计特性。

2010-3-15 控制理论与制导技术研究中心 第16页Harbin Institute of Technology– HIT三、算法程序示例（M语言）…… nb = 2; % 分子参数个数； na = 2; % 分母参数个数； N = 2; % 观测矩阵的行数； Fai_N = zeros(N,na+nb); % 对观测矩阵初始化0； % 对观测矩阵赋值 for i = 1 : nb for j = 1 : N Fai_N(j,na+nb+1-i) = R_input(i-1+j); end end for i = 1 : na for j = 1 : N Fai_N(j,na+1-i) = -R_out(i-1+j); end end 2010-3-15 控制理论与制导技术研究中心 第17页Harbin Institute of Technology– HIT% 赋值结束 Y_N = zeros(N,1); for i = 1: N Y_N(i) = R_out(na+i); end % 一次性完成算法 Saita = inv(Fai_N'*Fai_N)*Fai_N'*Y_N ……2010-3-15控制理论与制导技术研究中心第18页Harbin Institute of Technology– HIT四、算法演示参见：Lecture4文件夹中的M文件与mdl文件。

例一、稳定被控对象G( z) =0.45 z + 0.23 z 2 − 1.81z + 0.82y ( k ) − 1.81 y ( k − 1) + 0.82 y ( k − 2) = 0.45u (k − 1) + 0.23u (k − 2)2010-3-15控制理论与制导技术研究中心第19页Harbin Institute of Technology– HIT1. 考虑未加噪声干扰的情况 (B_Msequenc_10_20080314.m, Example20080314.mdl, B_leastsquare.m)输入信号10阶 M序列，幅值 ±1 , ∆t = 0.005 秒。

2010-3-15控制理论与制导技术研究中心第20页Harbin Institute of Technology– HITInput-Signal2 1.5 1Amplitude0.5 0 -0.5 -1 -1.5 -2 0 2 4 6 Time (Seconds)图3. 输入信号8102010-3-15控制理论与制导技术研究中心第21页Harbin Institute of Technology– HITOutput-Signal50 40 30Amplitude20 10 0 -10 -20 0 2 4 6 Time (Seconds)图4. 输出信号8102010-3-15控制理论与制导技术研究中心第22页Harbin Institute of Technology– HITa. N = 1000 时： Saita = -1.810000000000200e+000 8.200000000002071e-001 4.499999999999995e-001 2.299999999999036e-0012010-3-15控制理论与制导技术研究中心第23页Harbin Institute of Technology– HITb. N = 10时：⎡ 0.4500 ⎢ 0.5945 ⎢ ⎢ 0.0270 ⎢ ⎢ -1.1185 ⎢ -1.8267 Φ (10) = ⎢ ⎢ -1.7092 ⎢ -0.9157 ⎢ ⎢ -0.4759 ⎢ -0.7905 ⎢ ⎢ -1.7205 ⎣ 0 0.4500 0.5945 0.0270 1.0000 1.0000 1.0000 -1.0000 -1.0000 ⎤ 1.0000 ⎥ ⎥ 1.0000 ⎥ ⎥ 1.0000 ⎥ -1.0000 ⎥ ⎥ -1.0000 ⎥ -1.0000 ⎥ ⎥ 1.0000 ⎥ 1.0000 ⎥ ⎥ 1.0000 ⎥ ⎦-1.1185 -1.0000 -1.8267 -1.0000 -1.7092 -0.9157 -0.4759 1.0000 1.0000 1.0000-0.7905 -1.0000Y (10) = [-0.5945 -0.0270 1.1185 1.8267 1.7092 0.9157 0.4759 0.7905 1.7205 2.2460]T2010-3-15控制理论与制导技术研究中心第24页Harbin Institute of Technology– HITSaita = -1.809999999999999e+000 8.199999999999982e-001 4.499999999999996e-001 2.300000000000010e-0012010-3-15控制理论与制导技术研究中心第25页Harbin Institute of Technology– HITc. N = 4时： Saita = -1.809999999999996e+000 8.200000000000198e-001 4.499999999999904e-001 2.299999999999858e-0012010-3-15控制理论与制导技术研究中心第26页Harbin Institute of Technology– HITd. N = 3时 Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.926938e-018.2010-3-15控制理论与制导技术研究中心第27页Harbin Institute of Technology– HIT2. 考虑在数据输出端加上高斯白噪声序列 (B_Msequenc_10_20080314.m, Example20080314_n.mdl, B_leastsquare.m)2010-3-15控制理论与制导技术研究中心第28页Harbin Institute of Technology– HITInput-Signal2 1.5 1Amplitude0.5 0 -0.5 -1 -1.5 -2 0 2 4 6 Time (Seconds)图5. 输入信号 控制理论与制导技术研究中心8102010-3-15第29页Harbin Institute of Technology– HITOutput-Signal50 40 30Amplitude20 10 0 -10 -20 -30 0 2 4 6 Time (Seconds)图6. 输出信号 控制理论与制导技术研究中心8102010-3-15第30页Harbin Institute of Technology– HITa. 方差为 1时 Saita = -9.077092776149172e-001 -7.797576161866845e-002 4.807762051378330e-001 5.890277898736898e-0012010-3-15控制理论与制导技术研究中心第31页Harbin Institute of Technology– HITb. 方差为0.1时 Saita = -1.585198048771056e+000 5.955460280952445e-001 4.447524994244899e-001 3.111582289597879e-0012010-3-15控制理论与制导技术研究中心第32页Harbin Institute of Technology– HITc. 方差为0.01时 Saita = -1.785264848760534e+000 7.952709611711922e-001 4.467479*********e-001 2.345254993357168e-0012010-3-15控制理论与制导技术研究中心第33页Harbin Institute of Technology– HIT0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 0 2 4图8. 噪声信号68102010-3-15控制理论与制导技术研究中心第34页Harbin Institute of Technology– HITOutput-Signal50 40 30Amplitude20 10 0 -10 -20 0 2 4 6 Time (Seconds)图9. 输出信号8102010-3-15控制理论与制导技术研究中心第35页Harbin Institute of Technology– HIT3．考虑在系统的输入端加上高斯白噪声序列 (B_Msequenc_10_20080314.m,Example20080314_n1.mdl, B_leastsquare.m)2010-3-15控制理论与制导技术研究中心第36页Harbin Institute of Technology– HITa. 方差为1时 Saita = -1.834534885486862e+000 8.462761013050796e-001 4.550514915338836e-001 2.496264624523160e-0012010-3-15控制理论与制导技术研究中心第37页Harbin Institute of Technology– HITb. 方差为5时 Saita = -1.845194188501218e+000 8.582413447176447e-001 4.611804937181372e-001 2.824455404636859e-0012010-3-15控制理论与制导技术研究中心第38页Harbin Institute of Technology– HITc. 方差为10时 Saita = -1.847290220363152e+000 8.606997922255449e-001 4.658284739127945e-001 3.096861249074584e-0012010-3-15控制理论与制导技术研究中心第39页Harbin Institute of Technology– HIT例2 对临界稳定对象 (B_Msequenc_10_20080314.m, Example_unstable.mdl, B_leastsquare.m)0.45 z + 0.23 G( z) = 2 z − 1.82 z + 0.82y ( k ) − 1.82 y ( k − 1) + 0.82 y ( k − 2) = 0.45u (k − 1) + 0.23u (k − 2)2010-3-15控制理论与制导技术研究中心第40页Harbin Institute of Technology– HITSaita = -1.819999999988698e+000 8.199999999887081e-001 4.500000000000221e-001 2.300000000052445e-0012010-3-15控制理论与制导技术研究中心第41页Harbin Institute of Technology– HITInput-Signal2 1.5 1Amplitude0.5 0 -0.5 -1 -1.5 -2 0 2 4 6 Time (Seconds)图11. 输入信号 控制理论与制导技术研究中心8102010-3-15第42页Harbin Institute of Technology– HITOutput-Signal500Amplitude-50-100-150 024 6 Time (Seconds)图12. 输出信号 控制理论与制导技术研究中心8102010-3-15第43页Harbin Institute of Technology– HIT例三、对于不稳定系统 (B_Msequenc_10_20080314.m, Example_unstable1.mdl, B_leastsquare.m)0.45 z + 0.23 G( z) = 2 z − 1.83 z + 0.82Root1: 1.046244047484068e+000 Root2: 7.837559525159326e-0012010-3-15控制理论与制导技术研究中心第44页Harbin Institute of Technology– HITInput-Signal2 1.5 1Amplitude0.5 0 -0.5 -1 -1.5 -2 0 0.1 0.2 0.3 Time (Seconds)图12. 输入信号 控制理论与制导技术研究中心0.40.52010-3-15第45页Harbin Institute of Technology– HITOutput-Signal50 0Amplitude-50 -100 -150 -200 00.10.2 0.3 Time (Seconds)图13. 输出信号 控制理论与制导技术研究中心0.40.52010-3-15第46页Harbin Institute of Technology– HITN = 20时， Saita = -1.829999999999989e+000 8.199999999999954e-001 4.499999999999979e-001 2.300000000000012e-0012010-3-15控制理论与制导技术研究中心第47页Harbin Institute of Technology– HIT2010-3-15控制理论与制导技术研究中心第48页Harbin Institute of Technology– HIT2010-3-15控制理论与制导技术研究中心第49页Harbin Institute of Technology– HIT仿真结论： 1．最小二乘算法既可以应用到“稳定对象”也可以应用 到“不稳定被控对象”；对于自身不稳定的对象进行辨识 时，特别对于实际的不稳定对象，测试时要注意不能使 系统太远离平衡位置，以免超出线性方程所能描述的范 围，而实际情况，这一点很难保证； 2．噪声的方差对辨识精度影响很大，其影响程度对不同 的被控对象有所不同。