matlab tskikk,最小二乘法汇总及matlab仿真

最小二乘算法总结1 支强 zhistar@ 西安交通大学自动化系 2009 年 9 1 本文仅限用于学习交流，转载请注明出处。文中如有不对之处，欢迎指正！ 西安交通大学科技报告 zhistar@ 1 目录 1. 一般最小二乘法 ...............................................................................1 1.1. 一次计算最小二乘算法 .............................................................1 1.2. 递推最小二乘算法 .....................................................................1 2. 遗忘因子最小二乘算法 ...................................................................4 2.1. 一次计算法 .................................................................................4 2.2. 递推算法......................................................................................4 3. 限定记忆最小二乘递推算法 ...........................................................7 4. 偏差补偿最小二乘法 .......................................................................9 5. 增广最小二乘法 .............................................................................11 6. 广义最小二乘法 .............................................................................13 7. 辅助变量法 .....................................................................................15 8. 二步法 .............................................................................................17 9. 多级最小二乘法 .............................................................................19 10. Yule-Walker 辨识算法....................................................................21 Matlab 程序附录......................................................................................22 附录 1、最小二乘一次计算法 ...............................................................22 附录 2、最小二乘递推算法 ...................................................................23 附录 3、遗忘因子最小二乘一次计算法 ...............................................24 附录 4、遗忘因子最小二乘递推算法 ...................................................25 附录 5、限定记忆最小二乘递推算法 ...................................................27 附录 6、偏差补偿最小二乘递推算法 ...................................................29 附录 7、增广最小二乘递推算法 ...........................................................30 西安交通大学科技报告 zhistar@ 2 附录 8、广义最小二乘递推算法 ...........................................................32 附录 9、辅助变量法 ...............................................................................34 附录 10、二步法......................................................................................36 附录 11、多级最小二乘法......................................................................37 附录 12、Yule-Walker 辨识算法............................................................40 图 1 一般最小二乘参数过渡过程 .....................................................2 图 2 一般最小二乘方差变化过程 ....................................................3 图 3 遗忘因子法参数过渡过程 ........................................................5 图 4 遗忘因子法方差变化过程 ........................................................6 图 5 限定记忆法参数过渡过程 ........................................................8 图 6 限定记忆法方差变化过程 ........................................................8 图 7 偏差补偿最小二乘参数过渡过程..........................................10 图 8 偏差补偿最小二乘方差变化过程..........................................10 图 9 增广最小二乘辨识模型 ..........................................................11 图 10 增广最小二乘参数过渡过程 ................................................12 图 11 广义最小二乘参数过渡过程 ................................................14 图 12 广义最小二乘方差变化过程 ................................................14 图 13 辅助变量法参数过渡过程 ....................................................16 图 14 辅助变量法方差变化过程 ....................................................16 图 15 二步法参数过渡过程 ............................................................18 图 16 二步法方差变化过程 ............................................................18 1. 一般最小二乘法一般最小二乘法 例 1 考虑如下仿真对象(2) 1.5 (1)0.7 ( )(1)0.5 ( )( )z kz kz ku ku kv k+−++=+++ 其中，( )v k为服从(0,1)N分布的白噪声。输入信号( )u k采用 M 序列，幅度为 1。M 序列 由 9 级移位寄存器产生， 49iii xxx −− =⊕。 选择如下的辨识模型 1212 (2)(1)( )(1)( )( )z ka z ka z kbu kb u kv k+= −+−++++ 观测数据长度取400L =。加权阵取IΛ =。 1.1. 一次计算最小二乘算法一次计算最小二乘算法 ^ 1 ^ ^ 2 1 ^ 1 ^ 2 -1.4916 0.7005 () 1.0364 0.4268 TT LS LLLL a a H HH Z b b θ − ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ === ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (1.1) 其中， (3) (4) ... (402) L Z Z Z Z ⎛⎞ ⎜⎟ ⎜⎟ = ⎜⎟ ⎜⎟ ⎝⎠ ， (3)(2)(1)(2)(1) (3)(2)(3)(2)(4) ............... (401)(400)(401)(400) (402) T T L T hZZuu ZZuuh H ZZuu h ⎛⎞ −−⎛⎞ ⎜⎟ ⎜⎟ −− ⎜⎟ ⎜⎟ == ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ −− ⎝⎠ ⎝⎠ Matlab 程序见附录 1。 1.2. 递推最小二乘算法递推最小二乘算法 递推最小二乘算法公式： ^^^ 1 ( )(1)( )[ ( )( ) (1)] 1 ( )(1) ( )[ ( ) (1) ( )] ( ) ( )[( ) ( )] (1) kkK kz kh kk K kP kh k h k P kh k k P kIK k h k P k θθθ − =−+−− =−−+ Λ =−− (1.2) 西安交通大学科技报告 zhistar@ 2 初始条件 ^ 1 ^ ^ 2 ^ 1 ^ 2 3 3 (0) 3 3 a a b b θ ⎛⎞ ⎜⎟ ⎛ ⎞ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ == ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎝ ⎠ ⎜⎟ ⎝⎠ , 4 4 (0)100*PI =。 经过编程计算，各个参数的估计值为 ^ 1 ^ ^ 2 ^ 1 ^ 2 -1.4976 0.6802 1.0284 0.3341 LS a a b b θ ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ == ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (1.3) Matlab 程序见附录 2。 050100150200250300350400450 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 待估参数过渡过程 b1 a2 b2 a1 图 1 一般最小二乘参数过渡过程 西安交通大学科技报告 zhistar@ 3 050100150200250300350400450 0 10 20 30 40 50 60 70 80 90 100 估计方差变化过程 图 2 一般最小二乘方差变化过程 西安交通大学科技报告 zhistar@ 4 2. 遗忘因子最小二乘算法遗忘因子最小二乘算法 采用的辨识模型与例 1 相同。 2.1. 一次计算法一次计算法 ^ 1 ^ ^ 2 **1** ^ 1 ^ 2 -1.4990 0.7060 () 0.8260 0.4919 TT LS LLLL a a HHHZ b b θ − ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ === ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (2.1) 其中， 1 2 * * (3) * (4) ... (402) L L L Z Z Z Z β β − − ⎛⎞ ⎜⎟ ⎜⎟ = ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ， 1 1111 2 2222 *(3)* (2)* (1)* (2)* (1) *(4)* (3)* (2)* (3)* (2) ............... (401)(400)(401)(400)(402) LT LLLL LT LLLL L T hZZuu hZZuu H ZZuuh βββββ βββββ − −−−− − −−−− ⎛⎞⎛⎞ −− ⎜⎟⎜⎟ −− ⎜⎟⎜⎟ == ⎜⎟⎜⎟ ⎜⎟⎜⎟ ⎜⎟⎜⎟ −− ⎝⎠⎝⎠ 衰减因子0.98β=，数据长度402L =。 Matlab 程序见附录 3。 2.2. 递推算法递推算法 遗忘因子递推最小二乘算法公式： ^^^ 1 ( )(1)( )[ ( )( ) (1)] ( )(1) ( )[ ( ) (1) ( )] 1 ( )[( ) ( )] (1) kkK kz kh kk K kP kh k h k P kh k P kIK k h k P k θθθ − =−+−− =−−+ =−− (2.2) 西安交通大学科技报告 zhistar@ 5 其中，01≤≤为遗忘因子， 此处取0.98。 数据长度L=402， 初始条件 ^ 1 ^ ^ 2 ^ 1 ^ 2 0.001 0.001 (0) 0.001 0.001 a a b b θ ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ == ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ , 4 4 (0)10*PI =。 经过编程计算，各个参数的估计值为 ^ 1 ^ ^ 2 ^ 1 ^ 2 -1.4852 0.6720 1.0734 0.4387 LS a a b b θ ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ == ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (2.3) Matlab 程序见附录 4。 050100150200250300350400450 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 待估参数过渡过程 b1 a2 b2 a1 图 3 遗忘因子法参数过渡过程 西安交通大学科技报告 zhistar@ 6 050100150200250300350400450 0 1 2 3 4 5 6 7 8 9 10 估计方差变化过程 图 4 遗忘因子法方差变化过程 西安交通大学科技报告 zhistar@ 7 3. 限定记忆最小二乘递推算法限定记忆最小二乘递推算法 辨识模型与例 1 相同。 限定记忆最小二乘算法： ^^^ (1,)( ,)(1,)[ (1)(1) ( ,)]kkLk kLK kkL z kLh kLk kLθθθ++=+−++++−+++ 1 (1,)( ,) (1)[1(1) ( ,) (1)]K kkLP k kL h kLh kLP k kL h kL − ++=+++++++++ (1,)[(1,) (1)] ( ,)P kkLIK kkL h kLP k kL++=++++++ ^^^ ( ,)( ,1)( ,)[ ()() ( ,1)]k kLk kLK k kL z kLh kLk kLθθθ+=+−+++−++− (3.1) 1 ( ,)( ,1) ()[1() ( ,1) ()]K k kLP k kLh kLh kL P k kLh kL − +=+−++++−+ ( ,)[( ,) ()] ( ,1)P k kLIK k kL h kL P k kL+=−+++− 初始条件 ^ 1 ^ ^ 2 ^ 1 ^ 2 3 3 (0) 3 3 a a b b θ ⎛⎞ ⎜⎟ ⎛ ⎞ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ == ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎝ ⎠ ⎜⎟ ⎝⎠ , 4 4 (0)100*PI =。 经过编程计算，各个参数的估计值为 ^ 1 ^ ^ 2 ^ 1 ^ 2 -1.4858 0.6788 0.9431 0.6777 LS a a b b θ ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ == ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (3.2) Matlab 程序见附录 5。 西安交通大学科技报告 zhistar@ 8 050100150200250300350400 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 待估参数过渡过程 b1 a2 b2 a1 图 5 限定记忆法参数过渡过程 050100150200250300350400 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 估计方差变化过程 图 6 限定记忆法方差变化过程 西安交通大学科技报告 zhistar@ 9 4. 偏差补偿最小二乘法偏差补偿最小二乘法 辨识模型与例 1 相同。 偏差补偿最小二乘递推算法如下： ^^^ ( )(1)( )[ ( )( )(1)] LSLSLSkkK kz kh kkθθθ=−+−− 1 ( )(1) ( )[1( ) (1) ( )]K kP kh kh k P kh k − =−+− ( )[( ) ( )] (1)P kIK k h k P k=−− ^ 2 [ ( )( )(1)] ( )(1) 1( ) (1) ( ) LSz kh kk J kJ k h k P kh k θ−− =−+ +− (4.1) 2 ^ ^^ ( ) ( ) [1(1)( )] w cLS J k k kkDk σ θθ = +− 0 00 a b n n I D ⎛⎞ = ⎜ ⎟ ⎜⎟ ⎝⎠ 2 ^^^^ ( )( )( ) ( )(1) cLSwckkkk P k Dkθθσθ=+− 数据长度 L=402，初始条件 ^ 1 ^ ^ 2 ^ 1 ^ 2 3 3 (0) 3 3 a a b b θ ⎛⎞ ⎜⎟ ⎛ ⎞ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ == ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎝ ⎠ ⎜⎟ ⎝⎠ , 4 4 (0)100*PI =，偏差补偿初始值为 ^ 2 3 (0) 1 3.5 cθ ⎛⎞ ⎜⎟ ⎜⎟ = ⎜⎟ ⎜⎟ ⎝⎠ 。 经过编程计算，各个参数的估计值为 ^ 1 ^ ^ 2 ^ 1 ^ 2 -1.5066 0.7113 0.9332 0.5312 LS a a b b θ ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ == ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (4.2) Matlab 程序见附录 6。 西安交通大学科技报告 zhistar@ 10 050100150200250300350400450 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 待估参数过渡过程 b1 a2 b2 a1 图 7 偏差补偿最小二乘参数过渡过程 050100150200250300350400450 0 10 20 30 40 50 60 70 80 90 100 估计方差变化过程 图 8 偏差补偿最小二乘方差变化过程 西安交通大学科技报告 zhistar@ 11 5. 增广最小二乘法增广最小二乘法 例 2 考虑如下仿真对象 (2) 1.5 (1)0.7 ( )(1)0.5 ( )(1)0.2 ( )z kz kz ku ku kv kv k+−++=++−++ 其中，( )v k为服从(0,1)N分布的白噪声。输入信号( )u k采用 M 序列，幅度为 1。M 序列 由 9 级移位寄存器产生， 49iii xxx −− =⊕。 12 12 1.00.5 1 1.50.7 zz zz −− −− + −+ 12 12 0.2 1 1.50.7 zz zz −− −− −+ −+ ( )v k ( )M k ( )z k + + 图 9 增广最小二乘辨识模型 选择如下的辨识模型 121212 (2)(1)( )(1)( )(1)( )z ka z ka z kbu kb u kd v kd v k+= −+−++++++ 观测数据长度取402L =。加权阵取IΛ =。 解答： 增广最小二乘递推算法： ^^^ ( )(1)( )[ ( )( ) (1)]kkK kz kh kkθθθ=−+−− 1 ( )(1) ( )[ ( ) (1) ( ) 1]K kP kh k h k P kh k − =−−+ (5.1) ( )[( ) ( )] (1)P kIK k h k P k=−− 其中，() ^ 121212 T aabbddθ=， ()( )(1)(2)(1)(2)(1)(2) T h kz kz ku ku kv kv k= −−−−−−−−。 仿真初始条件： () ^^^^^^^ 121212(0)3 3 3 3 3 3 T T aabbddθ ⎛⎞ == ⎜⎟ ⎝⎠ , 6 6 (0)100*PI =。 经过计算，得到各个参数的估计值： () ^^^^^^ 121212-1.4999 0.7000 1.0001 0.5002 -0.9999 0.2000 T T aabbdd ⎛⎞ = ⎜⎟ ⎝⎠ 西安交通大学科技报告 zhistar@ 12 Matlab 程序见附录 7。 050100150200250300350400450 -2 -1 0 1 2 3 4 5 待估参数过渡过程 d1 a1 d2 b2 a2 b1 图 10 增广最小二乘参数过渡过程 西安交通大学科技报告 zhistar@ 13 6. 广义最小二乘法广义最小二乘法 例 模型结构选用 1212 12 ( )(1)(2)(1)(2)( ) ( )(1)(2)( ) z ka z ka z kbu kb u ke k e kc e kc e kv k +−+−=−+−+⎧ ⎨ +−+−= ⎩ 其中，各个参数的真值为： 1 1.5a = −， 1 0.7a =， 1 1b =， 2 0.5b =， 1 0c =， 2 0c =。 解答：广义最小二乘算法： ^^^ ( )(1)( )[( )( ) (1)] fff kkKkzkhkkθθθ=−+−− 1 ( )(1)( )[1( )(1)( )] ffffff KkP khkhk P khk − =−+− ( )[( )( )](1) ffff P kIKk hk P k=−− (6.1) ^^^^ ( )(1)( )[ ( )( )(1)] eee ee kkK k e kh kkθθθ=−+−− 1 ( )(1)( )[1( )(1)( )] eeeeee K kP kh kh k P kh k − =−+− ( )[( )( )](1) eeee P kIK k h k P k=−− 仿真初始条件： 初始条件 ^ 1 ^ ^ 2 ^ 1 ^ 2 3 3 (0) 3 3 a a b b θ ⎛⎞ ⎜⎟ ⎛ ⎞ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ == ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎝ ⎠ ⎜⎟ ⎝⎠ , 4 4 (0)100*PI =， ^ ^ 1 ^ 2 0.5 (0) 0.3 e c c θ ⎛⎞ ⎛⎞ ⎜⎟ ==⎜ ⎟ ⎜⎟ ⎝⎠ ⎝⎠ , 2 2 (0)10* e PI =。 经过计算，各个参数的估计值为 ^ 1 ^ 2 ^ 1 ^ 2 ^ 1 ^ 2 -1.5120 0.7218 1.0319 0.5374 0.0097 0.0197 a a b b c c ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ =⎜ ⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎜⎟ ⎝⎠ (6.2) Matlab 程序见附录 8。 西安交通大学科技报告 zhistar@ 14 050100150200250300350400 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 待估参数过渡过程 b1 a2 b2 a1 图 11 广义最小二乘参数过渡过程 050100150200250300350400 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 噪声传递系数估计结果 c1 c2 图 12 广义最小二乘方差变化过程 西安交通大学科技报告 zhistar@ 15 7. 辅助变量法辅助变量法 辨识模型与例 1 相同，只不过在此处噪声为有色噪声，其产生过程为 ( )( )0.5* (1)0.2* (2)e kv kv kv k=+−+−，( )v k为零均值的不相关随机噪声。 参照 Tally 法选取辅助变量 ( )() d x kz kn=−， d n为误差传递函数的阶数，此处为 2. 则有 *( ) [(2 1)(22)(1)(2)]h kz kz ku ku k= −−−−−−−−。 辅助变量法的递推公式可写成 ^^^ **1 ( )(1)( )[ ( )( ) (1)] ( )(1)( )[ ( ) (1)( ) 1] ( )[( ) ( )] (1) kkK kz kh kk K kP kh k h k P kh k P kIK k h k P k θθθ − =−+−− =−−+ =−− (7.1) 初始条件 ^ 1 ^ ^ 2 ^ 1 ^ 2 3 3 (0) 3 3 a a b b θ ⎛⎞ ⎜⎟ ⎛ ⎞ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ == ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎝ ⎠ ⎜⎟ ⎝⎠ , 4 4 (0)100*PI =。 经过编程计算，各个参数的估计值为 ^ 1 ^ ^ 2 ^ 1 ^ 2 -1.3809 0.6050 1.1231 0.5115 a a b b θ ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ == ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (7.2) Matlab 程序见附录 9。 西安交通大学科技报告 zhistar@ 16 050100150200250300350400 -20 -15 -10 -5 0 5 10 15 待估参数过渡过程 a1 b1 a2 b2 图 13 辅助变量法参数过渡过程 050100150200250300350400 -100 -50 0 50 100 150 200 估计方差变化过程 图 14 辅助变量法方差变化过程 西安交通大学科技报告 zhistar@ 17 8. 二步法二步法 辨 识 模 型 与 例 1 相 同 ，(2) 1.5 (1)0.7 ( )(1)0.5 ( )( )z kz kz ku ku ke k+−++=+++， ( )( )0.5* (1)0.2* (2)e kv kv kv k=+−+−，( )v k为零均值的不相关随机噪声。 辅助变量法的递推公式可写成 ^^^ **1 ( )(1)( )[ ( )( ) (1)] ( )(1)( )[ ( ) (1)( ) 1] ( )[( ) ( )] (1) kkK kz kh kk K kP kh k h k P kh k P kIK k h k P k θθθ − =−+−− =−−+ =−− (8.1) 其中， *( ) [(1)(2)(3)(4)]h kM kM kM kM k= −−−−−−。( )M k为输入 M 序 列。 初始条件 ^ 1 ^ ^ 2 ^ 1 ^ 2 3 3 (0) 3 3 a a b b θ ⎛⎞ ⎜⎟ ⎛ ⎞ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ == ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎝ ⎠ ⎜⎟ ⎝⎠ , 4 4 (0)100*PI =。 经过编程计算，各个参数的估计值为 ^ 1 ^ ^ 2 ^ 1 ^ 2 -1.3809 0.6050 1.1231 0.5115 a a b b θ ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ == ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (8.2) Matlab 程序见附录 10。 西安交通大学科技报告 zhistar@ 18 050100150200250300350400 -250 -200 -150 -100 -50 0 50 100 150 200 待估参数过渡过程 图 15 二步法参数过渡过程 050100150200250300350400 -400 -300 -200 -100 0 100 200 300 400 500 估计方差变化过程 图 16 二步法方差变化过程 西安交通大学科技报告 zhistar@ 19 9. 多级最小二乘法多级最小二乘法 [Hsia, 1977] 考虑如下仿真对象 ( )0.9 (1)0.15 (2)0.02 (3)0.7 (1) 1.5 (2)( ) ( ) 1.0 (1)0.41 (2)( ) z kz kz kz ku ku ke k e ke ke kv kλ +−+−+−=−−−+⎧ ⎨ +−+−= ⎩ 其中，( )u k是输入变量，此处为 M 序列；( )v k是零均值、方差为 1 的不相关随机噪声，通 过控制λ的大小来控制信噪比。 解答：辨识模型结构选用 11 1 1 () ( )() ( )( ) () A zz kB zu kv k C z −− − =+ (9.1) 其中， 1123 123 112 12 112 12 ()1 () ()1 A za za za z B zb zb z C zc zc z −−−− −−− −−− ⎧= +++ ⎪ =+ ⎨ ⎪ = ++ ⎩ (9.2) 辨识过程如下： 第一级，辅助模型参数辨识 原模型可写为 1111 () () ( )() () ( )( )A zC zz kB zC zu kv k −−−− =+ (9.3) 令 11112345 12345 1111234 1234 ()() ()1 ()() () E zA zC ze ze ze ze ze z F zB zC zf zf zf zf z −−−−−−−− −−−−−−− ⎧=+++++ ⎪ ⎨ =+++ ⎪⎩ ? ? (9.4) 利用最小二乘法可获得辅助模型的参数无偏一致估计值 ^ 1 1 111 1 ()H HH zθ − = (9.5) 其中， 1 [ (1)(2)...( )]Tzzzz L=， 1 1 1 (1) ... ( ) h H h L ⎛⎞ ⎜⎟ =⎜ ⎟ ⎜⎟ ⎝⎠ ，数据长度 L=400， 1123451234 1 [] ( )[(1)(2)(3)(4)(5)(1)(2) (3)(4)] T T eeeeeffff h kz kz kz kz kz ku ku k u ku k θ⎧= ⎪ = −−−−−−−−−−−−−− ⎨ ⎪ −−−− ⎩ (9.6) 第二级，过程模型参数辨识 1111 () ()() ()A zF zB zE z −−−− = (9.7) 西安交通大学科技报告 zhistar@ 20 根据最小二乘算法可以获得过程模型的参数估计值为 ^ 1 2 2222 ()H HH zθ − = (9.8) 其中， 21234 [000]Tzffff=， 212312 []Taaabbθ=， 11 2121 232132 43243 4354 45 00010 001 0 0 000 fe ffee Hfffee fffee ffee fe ⎛⎞ ⎜⎟ − ⎜⎟ ⎜⎟−− ⎜⎟ = −−− ⎜⎟ ⎜⎟ −−− ⎜⎟ −− ⎜⎟ ⎜⎟ − ⎝⎠ (9.9) 第三级，噪声模型参数辨识 111 111 ()() () ()() () E zA zC z F zB zC z −−− −−− = = (9.10) 根据最小二乘算法可以获得过程模型的参数估计值为 ^ 1 3 3333 ()H HH zθ − = 其中， 3112233452234 []Tzeaeaeaeefbff=−−−−， 312 []Tccθ=， 1 21 32 3 3 1 21 2 10 1 0 0 0 a aa aa H a b bb b ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ = ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ 。 经过计算，可以得到各级、各个参数的估计值如下： ^^^^^^^^^ 12345 1234 [] [1.85061.36760.52000.11810.04010.7616-0.8593-1.2119-0.5786] T T eeeeeffff= (9.11) ^^^^^ 12312[] [0.90580.13850.05990.7598-1.5767] T T aaabb= (9.12) 西安交通大学科技报告 zhistar@ 21 ^^ 12[][0.94630.3681] TT cc= (9.13) Matlab 程序见附录 11。 10. Yule-Walker 辨识算法辨识算法 考虑如下仿真对象( )0.9 (1)0.36 (2)0.054 (3)( )z kz kz kz kv k+−+−+−=，其中， ( )z k是可观测变量；( )v k是均值为零，方差为 1 的不相关随机噪声；数据长度取 L=1024。 相关函数按下式计算 1 1 ( )( ) () L z kj Rjz k z kj L = + =− ∑ (10.1) 参数的估计算法按下式计算 1 ^ 1 2 3 zz a aR r a θ − ⎛⎞ ⎜⎟ == − ⎜⎟ ⎜⎟ ⎝⎠ (10.2) ^^ (0) zz dRrθ=+ (10.3) 其中， [(1)(2)(3)] zzzz rRRR=， (0)(1)(2) (1)(0)(1) (2)(1)(0) zzz zzzz zzz RRR RRRR RRR ⎛⎞ ⎜⎟ = ⎜⎟ ⎜⎟ ⎝⎠ 。 经过计算，各个参数的估计值如下 ^ 1 ^ 2 ^ 3 ^ 0.8988 0.3368 0.0868 1.0143 a a a d ⎛⎞ ⎜⎟ ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ = ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎜⎟ ⎝⎠ (10.4) Matlab 程序见附录 12。 西安交通大学科技报告 zhistar@ 22 Matlab 程序附录程序附录 附录附录 1、最小二乘一次计算法、最小二乘一次计算法 Matlab 程序： clear clc %========================================== %最小二乘法辨识 % Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k) % %==========产生 M 序列作为输入=============== x=[0 1 0 1 1 0 1 1 1]; %initial value n=403; %n 为脉冲数目 M=[]; %存放 M 序列 for i=1:n temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end %产生均值为 0，方差为 1 的高斯白噪声 v=randn(1,400); z=[]; z(1)=-1; z(2)=0; for i=3:402 z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); end H=zeros(400,4); for i=1:400 H(i,1)=-z(i+1); H(i,2)=-z(i); H(i,3)=M(i+1); 西安交通大学科技报告 zhistar@ 23 H(i,4)=M(i); end Estimate=inv(H*H)*H*(z(3:402)) 附录附录 2、最小二乘递推算法、最小二乘递推算法 Matlab 程序： %最小二乘的递推算法 %Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k) %======================================== clear clc %==========400 个产生 M 序列作为输入=============== x=[0 1 0 1 1 0 1 1 1]; %initial value n=403; %n 为脉冲数目 M=[]; %存放 M 序列 for i=1:n temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end %===========产生均值为 0，方差为 1 的高斯白噪声========= v=randn(1,400); %==============产生观测序列 z================= z=zeros(402,1); z(1)=-1; z(2)=0; for i=3:402 z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); end %递推求解 P=100*eye(4); %估计方差 西安交通大学科技报告 zhistar@ 24 Pstore=zeros(4,401); Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; Theta=zeros(4,401); %参数的估计值，存放中间过程估值 Theta(:,1)=[3;3;3;3]; % K=zeros(4,400); %增益矩阵 K=[10;10;10;10]; for i=3:402 h=[-z(i-1);-z(i-2);M(i-1);M(i-2)]; K=P*h*inv(h*P*h+1); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h*Theta(:,i-2)); P=(eye(4)-K*h)*P; Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; end i=1:401; figure(1) plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title(待估参数过渡过程) figure(2) plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title(估计方差变化过程) 附录附录 3、遗忘因子最小二乘一次计算法、遗忘因子最小二乘一次计算法 clear clc %========================================== %最小二乘法辨识 % Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k) % %==========产生M序列作为输入=============== x=[0 1 0 1 1 0 1 1 1]; %initial value n=403; %n为脉冲数目 M=[]; %存放M序列 for i=1:n temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end 西安交通大学科技报告 zhistar@ 25 x(1)=temp; end %产生高斯白噪声 v=randn(1,400); z=[]; z(1)=-1; z(2)=0; u=0.98;% 遗忘因子 L=400; for i=3:402 z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); zstar(i)=z(i)*u^(L-i+2); end H=zeros(400,4); for i=1:400 H(i,1)=-z(i+1)*u^(L-i); H(i,2)=-z(i)*u^(L-i); H(i,3)=M(i+1)*u^(L-i); H(i,4)=M(i)*u^(L-i); end Estimate=inv(H*H)*H*(zstar(3:402)) 附录附录 4、遗忘因子最小二乘递推算法、遗忘因子最小二乘递推算法 %最小二乘的递推算法 %Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k) %======================================== clear clc %==========400 个产生 M 序列作为输入=============== x=[0 1 0 1 1 0 1 1 1]; %initial value n=403; %n 为脉冲数目 M=[]; %存放 M 序列 for i=1:n 西安交通大学科技报告 zhistar@ 26 temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end %===========产生均值为 0，方差为 1 的高斯白噪声========= v=randn(1,400); %==============产生观测序列 z================= z=zeros(402,1); z(1)=-1; z(2)=0; for i=3:402 z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); end %递推求解 P=10*eye(4); %估计方差 Pstore=zeros(4,401); Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; Theta=zeros(4,401); %参数的估计值，存放中间过程估值 Theta(:,1)=[0.001;0.001;0.001;0.001]; % K=zeros(4,400); %增益矩阵 K=[10;10;10;10]; u=0.98; %遗忘因子 for i=3:402 h=[-z(i-1);-z(i-2);M(i-1);M(i-2)]; K=P*h*inv(h*P*h+u); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h*Theta(:,i-2)); P=(eye(4)-K*h)*P/u; Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; end %==========================输出结果及作图============================= disp(参数 a1 a2 b1 b2 的估计值:) Theta(:,401) i=1:401; figure(1) plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) 西安交通大学科技报告 zhistar@ 27 title(待估参数过渡过程) figure(2) plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title(估计方差变化过程) 附录附录 5、限定记忆最小二乘递推算法、限定记忆最小二乘递推算法 Matlab 程序如下： %限定记忆最小二乘的递推算法 %Z(k+2)=1.5*Z(k+1)-0.7*Z(k)+u(k+1)+0.5*u(k)+v(k) %======================================== clear clc %==========400 个产生 M 序列作为输入=============== x=[0 1 0 1 1 0 1 1 1]; %initial value n=403; %n 为脉冲数目 M=[]; %存放 M 序列 for i=1:n temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end

展开阅读全文