下载文档原格式
《系统辨识基础》第20讲要点第6章 最小二乘类参数辨识方法6.1 引言最小二乘法是一种最基本的辨识方法，但如果模型的噪声不是白噪声，最小二乘法不一定能给出无偏、一致估计。
以下着重讨论模型噪声是有色噪声时的各种最小二乘辨识方法。
6.2 增广最小二乘法 6.2.1 增广最小二乘原理 考虑如下模型A z z kB z u k N z v k ()()()()()()---=+111式中u (k )和z (k ) 分别为模型输入和输出变量；v (k ) 是均值为零、方差为σv 2的不相关随机噪声或称白噪声；N z ()-1为噪声模型；A z ()-1 和B z ()-1为迟延算子多项式，记作A z a z a z a zB z b z b z b z n n n n a ab b()()--------=++++=+++⎧⎨⎪⎩⎪11122111221 其中n a 和n b 为模型阶次。
为了运用最小二乘原理来辨识这种模型的参数，需要把模型(4-1)式写成最小二乘格式)()()(k v k k z +=θτh这样就必须把噪声模型的参数包含在参数向量θ 中，从而引出增广概念，用来构造模型的参数向量θ和数据向量h(k )，具体的构成形式会因噪声模型的结构不同而不同。
下面是三种不同噪声模型的向量构成方法：① 若N z D z d z d z d z n n dd()()==++++----111221 ，可按下式构成参数向量和数据向量：⎪⎪⎩⎪⎪⎨⎧--------+==--------=∑∑∑===)(ˆ)1(ˆ)()1(ˆ)()1(ˆ)()(ˆ],,,,,,,,[)](ˆ,),1(ˆ),(,),1(),(,),1([)(111111i k v k d i k u k b i k z k a k z k v d d b b a a n k v k v n k u k u n k z k z k db a d b a n i i n i i n i i n n n d b a ττθ h② 若N z C zc zc zc zn n c c()()==++++----11111122，参数向量和数据向量的构成形式为：⎪⎪⎩⎪⎪⎨⎧-----+==----------=∑∑==)()1(ˆ)()1(ˆ)()(ˆ],,,,,,,,[)](ˆ,),1(ˆ),(,),1(),(,),1([)(11111i k u k b i k z k a k z k e c c b b a a n k e k e n k u k u n k z k z k ba cb a n i i n i i n n nc ba ττθ h③ 若N z D z C zd z d z d z c zc zc zn n n n d dc c()()()==++++++++--------111122112211 ，参数向量和数据向量的构成形式为：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧-----+=-----+==------------=∑∑∑∑====)()1(ˆ)()1(ˆ)()(ˆ)(ˆ)1(ˆ)(ˆ)1(ˆ)(ˆ)(ˆ],,,,,,,,,,,[)](ˆ,),1(ˆ),(ˆ,),1(ˆ),(,),1(),(,),1([)(11111111i k u k b i k z k a k z k e i k v k d i k e k c k e k v d d c c b b a a n k v k v n k e k e n k u k u n k z k z k b a dc cc b a n i i n i i n i i n i i n n n nd c b a ττθ h以上这种构成参数向量和数据向量的思想就是所谓的增广原理，它是增广最小二乘法的根本。
——第三章辨识方法
2 最小二乘法
1
清华大学电机系
辨识技术
辨识技术
()
*
lim k k α
α
→∞
=辨识技术
远离最小值点
接近最小值点
最速下降法
高斯牛顿法
清华大学电机系
辨识技术
42
2
三种算法的比较
y
香蕉函数最小化问题：
f (x ) = 100 ( X 2–X 12)2+ ( 1 –X 1)2
高斯-牛顿法：11次迭代
最速下降法：60次迭代
阻尼最小二乘法：18次迭代
清华大学电机系
辨识技术
43
例：非线性函数的参数辨识
x
C e
C x f ⋅⋅=21)(高斯-牛顿法：80次迭代
最速下降法：401次依然未收敛
阻尼最小二乘法：10次迭代
真值C * = [5.420136187 -0.25436189]
当x m = [ 1 2 3 4 5 ]时
测量值：f m (x ) = [ 4.20 3.25 2.52 1.95 1.51 ]f *(x ) = [ 4.202834 3.258924 2.527006 1.959469 1.519394 ]
真值
清华大学电机系
辨识技术
清华大学电机系辨识技术
清华大学电机系辨识技术
辨识技术
稳态电路。
---------------------------------------------------------------最新资料推荐------------------------------------------------------系统辨识—最小二乘法最小二乘法参数辨识 1 引言系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。
现代控制理论中的一个分支。
通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。
对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。
对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。
而系统辨识所研究的问题恰好是这些问题的逆问题。
通常，预先给定一个模型类={M}（即给定一类已知结构的模型），一类输入信号 u 和等价准则 J=L(y，yM)(一般情况下，J 是误差函数，是过程输出 y 和模型输出 yM 的一个泛函)；然后选择使误差函数J 达到最小的模型，作为辨识所要求的结果。
系统辨识包括两个方面：结构辨识和参数估计。
在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。
2 系统辨识的目的在提出和解决一个辨识问题时，明确最终使1 / 17用模型的目的是至关重要的。
它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。
通过辨识建立数学模型通常有四个目的。
①估计具有特定物理意义的参数有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。
这就需要通过能观测到的输入输出数据，用辨识的方法去估计那些参数。
②仿真仿真的核心是要建立一个能模仿真实系统行为的模型。
用于系统分析的仿真模型要求能真实反映系统的特性。
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．噪声的方差对辨识精度影响很大，其影响程度对不同 的被控对象有所不同。
