1、系统辨识大作业最小二乘法及其相关估值方法应用学号：2012302259姓名：王家琦基于最小二乘法的多种系统辨识方法研究1. 最小二乘法的引出在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。设单输入-单输出线性定长系统的差分方程为 ()+1(1)+()=0()+()， =1， 2， 3， (5.1.1)式中: 为随机干扰； 为理论上的输出值。 只有通过观测才能得到，在观测过程 () () ()中往往附加有随机干扰。 的观测值 可表示为() ()()=()+()(5.1.2)式中: 为随机干扰。由式(5.1.2)得 () ()=()()(5.1.3)将式(5.1.3) 带入式(5.1.1)
2、得()+1(1)+()=0()+1(1)+()+()+=1()(5.1.4)我们可能不知道 的统计特性，在这种情况下，往往把 看做均值为 0 的白噪声。() ()设(k)=n()+ =1()(5.1.5)则式(5.1.4) 可写成() =1(1)2(2)()+0()+1(1)+()+(k)(5.1.6)在观测 时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假()定 不仅包含了 的测量误差，而且包含了 的测量误差和系统内部噪声。假定(k) () ()是不相关随机序列(实际上 是相关随机序列)。(k) (k)现分别测出 个随机输入值+，则可写成 个方程，即(1)， (2)， ， (
3、+)， (1)， (2)， ， (+) N(+1) =1()2(1)(1)+0(+1)+1()+(1)+(+1)(+2) =1(+1)2()(2)+0(+2)+1(+1)+(2)+(+2)(+) =1(+1)2(+2)()+0(+)+1(+1)+()+(+)上述 个方程可写成向量-矩阵形式N(+1)(+2)(+)=()(+1) (1)(+1)(1)(2)(+2)(2) (+1)()(+)()10+(+1)(+2)(+3)(5.1.7)设=(+1)(+2)(+)， =10， =(+1)(+2)(+)= ()(+1) (1)(+1)(1)(2)(+2)(2) (+1)()(+)()则式(5.1.7
4、) 可写为 =+(5.1.8)式中： 为 维输出向量； 为 维噪声向量； 为 维参数向量； 为 测y N N (2n+1) N(2+1)量矩阵。因此式(5.1.8)是一个含有 个未知参数，由 个方程组成的联立方程组。如(2n+1) N果 ，方程数少于未知数数目，则方程组的解是不定的，不能唯一地确定参数向N2n+1量。如果 ，方程组正好与未知数数目相等，当噪声 时，就能准确地解出N=2n+1 =0=1(5.1.9)如果噪声 ，则0 =11(5.1.10)从上式可以看出噪声 对参数估计是有影响的，为了尽量较小噪声 对 估值的影响。在给定 输出向量 和测量矩阵 的条件下求系统参数 的估值，这就是系统
5、辨识问题。可用最小二 乘法来求 的估值，以下讨论最小二乘法估计。2. 最小二乘法估计算法设 表示 的最优估值， 表示 的最优估值，则有 =(5.1.11)=(+1)(+2)(+)， =10写出式(5.1.11)的某一行，则有() =1(1)2(2)()+0()+1(1)+()+(k)=1()+=0()， =+1， +2， ， +(5.1.12)设 表示 与 之差，即e(k)()()-e(k)= ()()=()-()=()=1()+=0()=(1+11+)()(0+11+)()=(1)()(1)()， =+1， +2， ， +(5.1.13)式中(1)=1+11+(1)=0+11+成为残差。把
6、分别代入式(5.1.13)可得残差e(k) =+1， +2， ， +。设e(+1)， e(+2)， ， e(+)=e(+1)e(+2)e(+)则有 e(k)= =(5.1.14)最小二乘估计要求残差的平方和为最小，即按照指数函数 = ee=()()(5.1.15)为最小来确定估值 。求 对 的偏导数并令其等于 0 可得 =2()=0(5.1.16)=(5.1.17)由式(5.1.17)可得 的最小二乘估计 =()1 (5.1.18)3.递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。设已获得的观测数据长度为 ，将式(5.1.8)中的 和 分别用 来代替，N ，
7、， 即 =+(5.3.1)用 的最小二乘估计，则表示 =()1 (5.3.2)设=（ ） 1(5.3.5)于是=(5.3.6)如果再获得 1 组新的观测值 和 ，则又增加 1 个方程(+1) (+1)+1=+1+1(5.3.7)式中 +1=(+1)， +1=(+1)+1=(+) (+1) (+1) (+1)将式(5.3.1) 和式(5.3.7) 合并，并写成分块矩阵形式，可得 +1= +1 + +1(5.3.8)根据上式可得到新的参数估值+1= +1 +11 +1 +1=+1 +1 +1=+1(+1+1)(5.3.9)式中+1= +1 +11=（ +1+1） 1根据矩阵求逆引理可以求得递推最小
8、二乘法辨识公式+1=+1(+1+1)(5.3.19)+1=+1+1(1+1+1)1(5.3.20)+1=+1(1+1+1)1+1(5.3.21)由于进行递推计算需要给出 和 初值 和 ，通过计算证明，可以取初值： ， 0 0 0=0，c 是充分大的常数， 为 单位矩阵，则经过若干次递推之后能够0=2 (2+1)(2+1)得到较好的参数估计。3. 辅助变量法辅助变量法是一种可克服最小二乘有偏估计的一种方法，对于原辨识方程 =+(5.4.1)当 是不相关随机序列时，最小二乘法可以得到参数向量 的一致无偏估计。但是，在实() 际应用中 往往是相关随机序列。()假定存在着一个 的矩阵 满足约束条件(2
9、n+1) Zlim1=0lim1=(5.4.2)式中 是非奇异的。用 乘以式(5.4.1) 等号两边得Q =+ (5.4.3)由上式得 =()1()1(5.4.4)如果取=()1(5.4.5)作为 估值，则称 为辅助变量估值，矩阵 成为辅助变量矩阵， 中的元素称为辅助变量。 常用的辅助变量法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。4. 广义最小二乘法广义最小二乘法是能克服最小二乘法有偏估计的另一种方法，这种方法计算比较复杂但效果比较好。下面直接介绍广义最小二乘法的计算步骤：(1)应用得到的输入和输出数据 和 ，按模型u(k)y(k)(=1， 2， 3， ， +)(1)()=(1)()
10、+()求出 的最小二乘估计 (1)=1(1)(1)0(1)(1)(2)计算残差 (1)() (1)()=(1)(1)()(1)(1)()(3)用残差 代替 ，计算(1)() () (1)=(1)(1)1(1)(1)(4)计算 和(1)()(1)()(1)()=()+(1)1(1)+(1)()(1)()=()+(1)1(1)+(1)()(5)应用得到的 和 按模型(1)()(1)()(1)(1)()=(1)(1)()+()用最小二乘法重新估计 ，得到 的第 2 次估值 。然后按步骤(2) 计算残差 ，按步 (2) (2)()骤(3)重新估计 ，得到估值 。再按照步骤(4)计算 和 ，按照步骤(5
11、)求 的第 (2) (2)()(2)() 3 次估值 。重复上述循环，之道 的估值 收敛为止。(3) ()5. 一种交替的广义最小二乘法求解技术(夏式法)这种方法是夏天长提出来的，又称夏式法。以上讨论过的广义最小二乘法的特点在于系统的输入和输出信号反复过滤。一下介绍的夏式法是一种交替的广义最小二乘法求解技术，它不需要数据反复过滤，因而计算效率较高。这种方法可消去最小二乘估计中的偏差，而且由这种方法导出的计算方法也比较简单。基于以上的几种方法，有=+=+(5.7.1)因而有=+=+(5.7.2)应用最小二乘法可得到参数估值=1(5.7.3)可以推出 =()1 ()1(5.7.11)上式中的第 1
12、 项是 最小二乘估计 ，第 2 项是偏差项 ，所以必须准确计算 。 为了准确计算 ，可采用迭代的方法。6. 专题解答设但输入-单输出系统的差分方程为 ()=1(1)2(2)+1(1)+2(2)+()()=()+1(1)+2(2)取真实值 ，输入数据如下所示=1 2 1 2=1.6420.7150.390.35k u(k) k u(k) k u(k)1 1.147 11 -0.958 21 0.4852 0.201 12 0.810 22 1.6333 -0.787 13 -0.044 23 0.0434 -1.159 14 0.947 24 1.3265 -1.052 15 -1.474 25
13、 1.7066 0.866 16 -0.719 26 -0.3407 1.152 17 -0.086 27 0.8908 1.573 18 -1.099 28 1.1449 0.626 19 1.450 29 1.17710 0.433 20 1.151 30 -0.390用 的真实值利用查分方程求出 作为测量值， 为均值为 0，方差为 0.1，0.5 的 () ()不相关随机序列。(1) 用最小二乘法估计参数 。=1 2 1 2(2) 用递推最小二乘法估计 。=1 2 1 2(3) 用辅助变量法估计参数 。=1 2 1 2(4) 设 ，用广义最小二乘法估计参数 。()+2()=() =1 2
14、 1 2(5) 用夏式法估计参数 。=1 2 1 2(6) 详细分析和比较所获得的参数辨识结果，并说明上述参数便是方法的优缺点。根据题目要求的解法，利用 Matlab 编程实现系统辨识的估值利用最小二乘法估计的结果如下：最小二乘法方差 1 2 b1 b20.0001 1.6280 0.7028 0.3971 3.44916316 0.7059 0.3947 0.34941.6354 0.7120 0.3918 0.34631.6362 0.7082 0.3970 0.35271.6360 0.7165 0.3906 0.35341.6289 0.7046 0.3908 0.34370.001
15、1.1543 0.6766 0.4064 0.33041.5577 0.6371 0.3868 0.32491.6050 0.6860 0.3737 0.32441.6060 0.6816 0.3583 0.31671.6195 0.7030 0.3907 0.33661.5670 0.6572 0.3752 0.31400.01 1.3538 0.5010 0.3486 0.17090.8956 0.1637 0.4237 0.06971.0008 0.2036 0.4236 0.12681.3403 0.4707 0.3826 0.26151.0574 0.2289 0.3682 0.13
16、171.1231 0.2963 0.3592 0.15060.1 1.1424 0.2710 0.3284 0.22161.0255 0.1736 0.3766 0.18440.8896 0.1109 0.3893 0.08130.8182 0.1114 0.4298 0.09230.8100 0.0153 0.4122 0.13520.7715 0.1311 0.4714 0.06400.5 0.9751 0.1017 0.0271 0.07920.8938 0.0740 0.3484 0.16340.1927 0.0197 0.3762 0.15210.5506 0.0392 0.6510
17、 0.05840.7560 0.0494 0.3372 0.14370.9459 0.1377 0.3815 0.1853部分程序运行结果递推最小二乘法方差 1 2 b1 b20.0001 1.6754 0.6787 0.5207 0.39241.3076 0.2900 0.0075 0.32841.5146 0.6963 1.1401 0.16391.5733 0.7782 0.4149 0.71101.1602 0.4753 0.6736 0.31591.2091 0.3192 0.5277 0.02750.001 1.4767 0.4040 0.2679 0.55121.6259 0.7
18、594 0.3253 0.37801.5393 0.4757 0.1268 0.43461.1548 0.1700 0.1926 0.82510.8858 0.0760 0.3385 0.04061.4129 0.3127 0.0992 0.83800.01 1.3485 0.3445 0.3194 0.37101.1639 0.3296 0.7813 0.21621.9946 1.2323 1.4852 0.03041.3924 0.3543 0.3319 0.45721.3982 0.3608 0.7773 0.31521.6346 0.7229 0.5780 0.39470.1 1.56
19、24 0.7132 0.4422 0.41121.7335 0.7152 0.0844 0.63991.4763 0.5366 0.3255 0.31161.4477 0.3489 0.2218 0.22651.6216 0.7082 0.6595 0.42751.5105 0.4000 0.0113 0.22130.5 1.7927 0.9411 0.2730 0.34711.5556 0.8877 0.5972 0.12171.7868 1.2538 1.1248 0.21001.5733 0.7434 0.3589 0.13871.3193 0.6084 1.2971 0.30291.5
20、959 0.5386 0.0141 0.6947部分程序运行结果：辅助变量法方差 1 2 b1 b20.0001 1.7799 0.8588 0.4147 0.40461.3076 0.2900 0.0075 0.32841.6735 0.7578 0.4060 0.34841.5812 0.6546 0.3771 0.36111.6657 0.7469 0.3772 0.35611.5281 0.6509 0.3645 0.33020.001 1.6295 0.6775 0.4082 0.35011.6425 0.7305 0.3937 0.35431.5595 0.6052 0.3563
21、0.33621.4145 0.4925 0.4021 0.28731.6371 0.7270 0.3681 0.34181.3539 0.4733 0.3906 0.24890.01 1.3451 0.4700 0.3822 0.26801.3657 0.4893 0.4349 0.24961.3702 0.5009 0.4388 0.26111.1884 0.3707 0.3457 0.15211.3636 0.5330 0.4135 0.20641.3158 0.4836 0.4562 0.21690.1 1.5545 0.6167 0.4104 0.35121.5900 0.6648 0
22、.4098 0.36641.6610 0.7029 0.4001 0.34951.5104 0.6014 0.3968 0.32561.5620 0.6496 0.3950 0.30971.4418 0.5706 0.4183 0.27830.5 1.4952 0.5704 0.3769 0.34831.5592 0.6541 0.4387 0.33301.3637 0.5302 0.3865 0.19161.5543 0.6324 0.3208 0.29291.5385 0.5994 0.3661 0.35761.3511 0.4839 0.3412 0.2208部分程序运行结果：广义最小二
23、乘法方差 1 2 b1 b20.0001 1.6451 0.7182 0.3895 0.35031.6468 0.7205 0.3866 0.34921.6308 0.7038 0.3852 0.34571.6467 0.7188 0.3898 0.35441.6411 0.7156 0.3898 0.34771.6441 0.7185 0.3896 0.34880.001 1.7027 0.7678 0.4059 0.37611.6477 0.7166 0.3973 0.36431.6391 0.7113 0.3836 0.35041.6586 0.7166 0.3791 0.36721.6
24、831 0.7472 0.3994 0.37311.5964 0.6850 0.3655 0.32310.01 1.6826 0.7469 0.3859 0.36001.7245 0.7696 0.3565 0.36421.6577 0.7186 0.3910 0.38201.6656 0.7263 0.3425 0.34131.6903 0.7392 0.3633 0.38381.6942 0.7439 0.3834 0.38800.1 0.8996 0.1515 0.4435 0.06041.0380 0.3044 0.3410 0.04641.6571 0.6569 0.4589 0.4
25、2591.2180 0.4584 0.3152 0.10300.1957 0.0974 0.2164 0.10351.6859 0.7232 0.4144 0.44040.5 1.4053 0.6703 0.3672 0.13721.5783 0.6148 0.3512 0.41731.0233 0.0658 0.2035 0.40601.5973 0.8103 0.3544 0.21311.6191 0.6709 0.3111 0.48661.0804 0.2210 0.7500 0.5185部分程序运行结果：夏式法方差 1 2 b1 b20.0001 1.6367 0.7090 0.391
26、0 0.35111.6451 0.7184 0.3876 0.34721.3618 0.2159 0.8205 0.58511.6291 0.7075 0.3911 0.34641.6173 0.6929 0.3955 0.34391.6307 0.7066 0.3912 0.34490.001 1.5314 0.6407 0.4016 0.30921.5632 0.6575 0.4204 0.33191.6172 0.6946 0.3925 0.32581.5670 0.6489 0.4016 0.32331.5326 0.6243 0.3995 0.30701.6251 0.6946 0.
27、3621 0.33600.01 0.8815 0.1089 0.3734 0.04401.3763 0.4972 0.4121 0.26541.5027 0.5952 0.3854 0.32771.3421 0.4268 0.3978 0.27901.2454 0.3024 0.3669 0.25251.3361 0.4741 0.4999 0.29000.1 0.6639 0.1074 0.6847 0.08271.1406 0.3154 0.4259 0.23670.1791 0.2916 0.9454 0.21261.2068 0.3827 0.7460 0.32200.6412 0.0
28、704 0.5802 0.08771.0384 0.2878 0.7493 0.37930.5 0.9330 0.0752 0.9741 0.58820.6291 0.1686 0.4679 0.78690.5986 0.2792 0.2605 0.20080.4441 0.8033 0.7442 1.76670.6216 0.1440 0.5821 0.49060.6576 0.0999 1.0002 0.0263部分程序运行结果：结论：通过编程计算，获得在噪声方差比较小的情况下，各种方法所获得的估值比较理想，但随着噪声方差的增大，估值的偏差随之增大，横向比较看来夏式法与广义最小二乘法能够更
29、好地还原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值。Matlab 源程序：%最小二乘估计%clearu= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177 -0.390;n=normrnd(0, sqrt(0.1), 1, 31);z=zeros(1,30);for
30、 k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);endh0=-z(2) -z(1) u(2) u(1);HLT=h0,zeros(4,28);for k=3:30h1=-z(k) -z(k-1) u(k) u(k-1);HLT(:,k-1)=h1;end HL=HLT;y=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(17);z(18);z(19
31、);z(20);z(21);z(22);z(23);z(24);z(25);z(26);z(27);z(28);z(29);z(30);z(31);%求出FAIc1=HL*HL;c2=inv(c1);c3=HL*y;c=c2*c3;%display(方差=0.1时，最小二乘法估计辨识参数如下：);a1=c(1);a2=c(2);b1=c(3);b2=c(4);clear%递推最小二乘法估计u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.71
32、9 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177 -0.390;z(2)=0;z(1)=0;n=normrnd(0, sqrt(0.1), 1, 31);for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715*n(k-2);endc0=0.001 0.001 0.001 0.001; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(4
33、,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=c0,zeros(4,30); %被辨识参数矩阵的初始值及大小e=zeros(4,30); %相对误差的初始值及大小for k=3:30; %开始求K h1=-z(k-1),-z(k-2),u(k-1),u(k-2); x=h1*p0*h1+1; x1=inv(x); %开始求K(k)k1=p0*h1*x1;%求出K的值d1=z(k)-h1*c0; c1=c0+k1*d1; %求被辨识参数ce1=c1-c0; %求参数当前值与上一次的值的差值e2=e1./c0;
34、%求参数的相对变化e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1*h1*p0*h1+1; %求出 p(k)的值p0=p1; %给下次用if e2=E break; %如果参数收敛情况满足要求，终止计算endend%display(方差为0.0001递推最小二乘法辨识后的结果是：);a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%display(a1,a2,b1,b2经过递推最小二乘法辨识的结果是
35、：);for i=3:31;if(c(1,i)=0)q1=c(1,i-1);break;endendfor i=3:31;if(c(2,i)=0)q2=c(2,i-1);break;endendfor i=3:31;if(c(3,i)=0)q3=c(3,i-1);break;endendfor i=3:31;if(c(4,i)=0)q4=c(4,i-1);break;endenda1=q1;a2=q2;b1=q3 ;b2=q4;%clear%辅助变量递推最小二乘法估计na=2; nb=2; siitt=1.642 0.715 0.39 0.35; siit0=0.001*eye(na+nb,1
36、); p=106*eye(na+nb); siit(:,1)=siit0; y(2)=0;y(1)=0; x(1)=0;x(2)=0; j=0; u= 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177 -0.390;n=normrnd(0, sqrt(0.01), 1, 31);f
37、or k=3:31;h=-y(k-1),-y(k-2),u(k-1),u(k-2); y(k)=h*siitt+n(k)+1.642*n(k-1)+0.715*n(k-2); hx=-x(k-1),-x(k-2),u(k-1),u(k-2); kk=p*hx/(h*p*hx+1); p=eye(na+nb)-kk*h*p; siit(:,k-1)=siit0+kk*y(k)-h*siit0;x(k)=hx*siit(:,k-1); j=j+(y(k)-h*1.642 0.715 0.39 0.35)2; e=max(abs(siit(:,k-1)-siit0)./siit0);ess(:,k-
38、2)=siit(:,k-1)-siitt; siit0=siit(:,k-1); enda1=siit0(1);a2=siit0(2);b1=siit0(3);b2=siit0(4);clear%广义最小二乘估计clear;nn = normrnd(0,sqrt(0.5),1,31);uk=1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1
39、.706 -0.340 0.890 1.144 1.177 -0.390;yk(1)=0;yk(2)=0;for i=1:29;yk(i+2)=-1.642*yk(i+1)-0.715*yk(i)+0.39*uk(i+1)+0.35*uk(i)+nn(i+2)+1.642*nn(i+1)+0.715*nn(i);end;for i=1:29;A(i,:)=-yk(i+1) -yk(i) uk(i+1) uk(i);endsiit=inv(A*A)*A*(yk(3:31)+nn(2:30);e(1)=yk(1);e(2)=yk(2)+siit(1)*yk(1)-siit(3)*uk(1);for
40、 i=3:31;e(i)=yk(i)+siit(1)*yk(i-1)+siit(2)*yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2);endfor i=1:29;fai(i,:)=-e(i+1) -e(i);endf=inv(fai*fai)*fai*e(3:31);for i=3:31;yk(i)=yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);endyk(2)=yk(2)+f(1)*yk(1);for i=3:30;uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);
41、for j=1:30for i=1:29;A(i,:)=-yk(i+1) -yk(i) uk(i+1) uk(i);endsiit=inv(A*A)*A*yk(3:31);e(1)=yk(1);e(2)=yk(2)+siit(1)*(yk(1)-siit(3)*uk(1);for i=3:31;e(i)=yk(i)+siit(1)*(yk(i-1)+siit(2)*(yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2);endfor i=1:29;fai(i,:)=-e(i+1) -e(i);endf=inv(fai*fai)*fai*e(3:31);k1(j)=f(1);k2(j
