系统辨识.docx
《系统辨识.docx》由会员分享，可在线阅读，更多相关《系统辨识.docx（24页珍藏版）》请在冰点文库上搜索。
系统辨识
5.3基于最小二乘法的多种系统辨识方法研究
设但输入-单输出系统的差分方程为
取真实值
，
取真实值，输入数据如下所示
k
u（k）
k
u（k）
k
u（k）
1
1.147
11
-0.958
21
0.485
2
0.201
12
0.810
22
1.633
3
-0.787
13
-0.044
23
0.043
4
-1.159
14
0.947
24
1.326
5
-1.052
15
-1.474
25
1.706
6
0.866
16
-0.719
26
-0.340
7
1.152
17
-0.086
27
0.890
8
1.573
18
-1.099
28
1.144
9
0.626
19
1.450
29
1.177
10
0.433
20
1.151
30
-0.390
用的真实值利用查分方程求出作为测量值，为均值为0，方差为0.1，0.5的不相关随机序列。
∙用最小二乘法估计参数。
∙用递推最小二乘法估计。
∙用辅助变量法估计参数。
∙设，用广义最小二乘法估计参数。
∙用夏式法估计参数。
∙详细分析和比较所获得的参数辨识结果，并说明上述参数便是方法的优缺点。
根据题目要求的解法，利用Matlab编程实现系统辨识的估值
利用最小二乘法估计的结果如下：
最小二乘法
方差
0.0001
1.6280
0.7028
0.3971
3.449
1..6316
0.7059
0.3947
0.3494
1.6354
0.7120
0.3918
0.3463
1.6362
0.7082
0.3970
0.3527
1.6360
0.7165
0.3906
0.3534
1.6289
0.7046
0.3908
0.3437
0.001
1.1543
0.6766
0.4064
0.3304
1.5577
0.6371
0.3868
0.3249
1.6050
0.6860
0.3737
0.3244
1.6060
0.6816
0.3583
0.3167
1.6195
0.7030
0.3907
0.3366
1.5670
0.6572
0.3752
0.3140
0.01
1.3538
0.5010
0.3486
0.1709
0.8956
0.1637
0.4237
0.0697
1.0008
0.2036
0.4236
0.1268
1.3403
0.4707
0.3826
0.2615
1.0574
0.2289
0.3682
0.1317
1.1231
0.2963
0.3592
0.1506
0.1
1.1424
0.2710
0.3284
0.2216
1.0255
0.1736
0.3766
0.1844
0.8896
0.1109
0.3893
0.0813
0.8182
0.1114
0.4298
0.0923
0.8100
0.0153
0.4122
0.1352
0.7715
0.1311
0.4714
0.0640
0.5
0.9751
0.1017
0.0271
0.0792
0.8938
0.0740
0.3484
0.1634
0.1927
0.0197
0.3762
0.1521
0.5506
0.0392
0.6510
0.0584
0.7560
0.0494
0.3372
0.1437
0.9459
0.1377
0.3815
0.1853
部分程序运行结果
递推最小二乘法
方差
0.0001
1.6754
0.6787
0.5207
0.3924
1.3076
0.2900
0.0075
0.3284
1.5146
0.6963
1.1401
0.1639
1.5733
0.7782
0.4149
0.7110
1.1602
0.4753
0.6736
0.3159
1.2091
0.3192
0.5277
0.0275
0.001
1.4767
0.4040
0.2679
0.5512
1.6259
0.7594
0.3253
0.3780
1.5393
0.4757
0.1268
0.4346
1.1548
0.1700
0.1926
0.8251
0.8858
0.0760
0.3385
0.0406
1.4129
0.3127
0.0992
0.8380
0.01
1.3485
0.3445
0.3194
0.3710
1.1639
0.3296
0.7813
0.2162
1.9946
1.2323
1.4852
0.0304
1.3924
0.3543
0.3319
0.4572
1.3982
0.3608
0.7773
0.3152
1.6346
0.7229
0.5780
0.3947
0.1
1.5624
0.7132
0.4422
0.4112
1.7335
0.7152
0.0844
0.6399
1.4763
0.5366
0.3255
0.3116
1.4477
0.3489
0.2218
0.2265
1.6216
0.7082
0.6595
0.4275
1.5105
0.4000
0.0113
0.2213
0.5
1.7927
0.9411
0.2730
0.3471
1.5556
0.8877
0.5972
0.1217
1.7868
1.2538
1.1248
0.2100
1.5733
0.7434
0.3589
0.1387
1.3193
0.6084
1.2971
0.3029
1.5959
0.5386
0.0141
0.6947
部分程序运行结果：
辅助变量法
方差
0.0001
1.7799
0.8588
0.4147
0.4046
1.3076
0.2900
0.0075
0.3284
1.6735
0.7578
0.4060
0.3484
1.5812
0.6546
0.3771
0.3611
1.6657
0.7469
0.3772
0.3561
1.5281
0.6509
0.3645
0.3302
0.001
1.6295
0.6775
0.4082
0.3501
1.6425
0.7305
0.3937
0.3543
1.5595
0.6052
0.3563
0.3362
1.4145
0.4925
0.4021
0.2873
1.6371
0.7270
0.3681
0.3418
1.3539
0.4733
0.3906
0.2489
0.01
1.3451
0.4700
0.3822
0.2680
1.3657
0.4893
0.4349
0.2496
1.3702
0.5009
0.4388
0.2611
1.1884
0.3707
0.3457
0.1521
1.3636
0.5330
0.4135
0.2064
1.3158
0.4836
0.4562
0.2169
0.1
1.5545
0.6167
0.4104
0.3512
1.5900
0.6648
0.4098
0.3664
1.6610
0.7029
0.4001
0.3495
1.5104
0.6014
0.3968
0.3256
1.5620
0.6496
0.3950
0.3097
1.4418
0.5706
0.4183
0.2783
0.5
1.4952
0.5704
0.3769
0.3483
1.5592
0.6541
0.4387
0.3330
1.3637
0.5302
0.3865
0.1916
1.5543
0.6324
0.3208
0.2929
1.5385
0.5994
0.3661
0.3576
1.3511
0.4839
0.3412
0.2208
部分程序运行结果：
广义最小二乘法
方差
0.0001
1.6451
0.7182
0.3895
0.3503
1.6468
0.7205
0.3866
0.3492
1.6308
0.7038
0.3852
0.3457
1.6467
0.7188
0.3898
0.3544
1.6411
0.7156
0.3898
0.3477
1.6441
0.7185
0.3896
0.3488
0.001
1.7027
0.7678
0.4059
0.3761
1.6477
0.7166
0.3973
0.3643
1.6391
0.7113
0.3836
0.3504
1.6586
0.7166
0.3791
0.3672
1.6831
0.7472
0.3994
0.3731
1.5964
0.6850
0.3655
0.3231
0.01
1.6826
0.7469
0.3859
0.3600
1.7245
0.7696
0.3565
0.3642
1.6577
0.7186
0.3910
0.3820
1.6656
0.7263
0.3425
0.3413
1.6903
0.7392
0.3633
0.3838
1.6942
0.7439
0.3834
0.3880
0.1
0.8996
0.1515
0.4435
0.0604
1.0380
0.3044
0.3410
0.0464
1.6571
0.6569
0.4589
0.4259
1.2180
0.4584
0.3152
0.1030
0.1957
0.0974
0.2164
0.1035
1.6859
0.7232
0.4144
0.4404
0.5
1.4053
0.6703
0.3672
0.1372
1.5783
0.6148
0.3512
0.4173
1.0233
0.0658
0.2035
0.4060
1.5973
0.8103
0.3544
0.2131
1.6191
0.6709
0.3111
0.4866
1.0804
0.2210
0.7500
0.5185
部分程序运行结果：
夏式法
方差
0.0001
1.6367
0.7090
0.3910
0.3511
1.6451
0.7184
0.3876
0.3472
1.3618
0.2159
0.8205
0.5851
1.6291
0.7075
0.3911
0.3464
1.6173
0.6929
0.3955
0.3439
1.6307
0.7066
0.3912
0.3449
0.001
1.5314
0.6407
0.4016
0.3092
1.5632
0.6575
0.4204
0.3319
1.6172
0.6946
0.3925
0.3258
1.5670
0.6489
0.4016
0.3233
1.5326
0.6243
0.3995
0.3070
1.6251
0.6946
0.3621
0.3360
0.01
0.8815
0.1089
0.3734
0.0440
1.3763
0.4972
0.4121
0.2654
1.5027
0.5952
0.3854
0.3277
1.3421
0.4268
0.3978
0.2790
1.2454
0.3024
0.3669
0.2525
1.3361
0.4741
0.4999
0.2900
0.1
0.6639
0.1074
0.6847
0.0827
1.1406
0.3154
0.4259
0.2367
0.1791
0.2916
0.9454
0.2126
1.2068
0.3827
0.7460
0.3220
0.6412
0.0704
0.5802
0.0877
1.0384
0.2878
0.7493
0.3793
0.5
0.9330
0.0752
0.9741
0.5882
0.6291
0.1686
0.4679
0.7869
0.5986
0.2792
0.2605
0.2008
0.4441
0.8033
0.7442
1.7667
0.6216
0.1440
0.5821
0.4906
0.6576
0.0999
1.0002
0.0263
部分程序运行结果：
结论：
通过编程计算，获得在噪声方差比较小的情况下，各种方法所获得的估值比较理想，但随着噪声方差的增大，估值的偏差随之增大，横向比较看来夏式法与广义最小二乘法能够更好地还原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值。
Matlab源程序：
%最小二乘估计%
clear
u=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];
n=normrnd（0,sqrt（0.1）,1,31）;
z=zeros（1,30）;
fork=3:
31
z（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）;
end
h0=[-z
（2）-z
（1）u
（2）u
（1）]';
HLT=[h0,zeros（4,28）];
fork=3:
30
h1=[-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）;z（20）;z（21）;z（22）;z（23）;z（24）;z（25）;z（26）;z（27）;z（28）;z（29）;z（30）;z（31）];%求出FAI
c1=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.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];
z
（2）=0;
z
（1）=0;
n=normrnd（0,sqrt（0.1）,1,31）;
fork=3:
31
z（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）;
end
c0=[0.0010.0010.0010.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量
p0=10^6*eye（4,4）;%直接给出初始状态P0，即一个充分大的实数单位矩阵
E=0.000000005;%取相对误差E=0.000000005
c=[c0,zeros（4,30）];%被辨识参数矩阵的初始值及大小
e=zeros（4,30）;%相对误差的初始值及大小
fork=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;%求被辨识参数c
e1=c1-c0;%求参数当前值与上一次的值的差值
e2=e1./c0;%求参数的相对变化
e（:
k）=e2;%把当前相对变化的列向量加入误差矩阵的最后一列
c0=c1;%新获得的参数作为下一次递推的旧参数
c（:
k）=c1;%把辨识参数c列向量加入辨识参数矩阵的最后一列
p1=p0-k1*k1'*[h1'*p0*h1+1];%求出p（k）的值
p0=p1;%给下次用
ife2<=E
break;%如果参数收敛情况满足要求，终止计算
end
end
%display（'方差为0.0001递推最小二乘法辨识后的结果是：
'）;
a1=c（1,:
）;
a2=c（2,:
）;
b1=c（3,:
）;
b2=c（4,:
）;
%display（'a1,a2,b1,b2经过递推最小二乘法辨识的结果是：
'）;
fori=3:
31;
if（c（1,i）==0）
q1=c（1,i-1）;
break;
end
end
fori=3:
31;
if（c（2,i）==0）
q2=c（2,i-1）;
break;
end
end
fori=3:
31;
if（c（3,i）==0）
q3=c（3,i-1）;
break;
end
end
fori=3:
31;
if（c（4,i）==0）
q4=c（4,i-1）;
break;
end
end
a1=q1;
a2=q2;
b1=q3;
b2=q4;
%clear
%辅助变量递推最小二乘法估计
na=2;
nb=2;
siitt=[1.6420.7150.390.35]';
siit0=0.001*eye（na+nb,1）;
p=10^6*eye（na+nb）;
siit（:
1）=siit0;
y
（2）=0;y
（1）=0;
x
（1）=0;x
（2）=0;
j=0;
u=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9850.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8900.1441.177-0.390];
n=normrnd（0,sqrt（0.01）,1,31）;
fork=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.6420.7150.390.35]'）^2;
e=max（abs（（siit（:
k-1）-siit0）./siit0））;
ess（:
k-2）=siit（:
k-1）-siitt;
siit0=siit（:
k-1）;
end
a1=siit0
（1）;
a2=siit0
（2）;
b1=siit0（3）;
b2=siit0（4）;
clear
%广义最小二乘估计
clear;
nn=normrnd（0,sqrt（0.5）,1,31）';
uk=[1.1470.201-0.787-1.589-1.0520.8661.1521.5730.6260.433-0.9580.810-0.0440.947-1.474-0.719-0.086-1.0991.4501.1510.4851.6330.0431.3261.706-0.3400.8901.1441.177-0.390];
yk
（1）=0;
yk
（2）=0;
fori=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;
fori=1:
29;
A（i,:
）=[-yk（i+1）-yk（i）uk（i+1）uk（i）];
end
siit=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）;
fori=3:
31;
e（i）=yk（i）+siit
（1）*yk（i-1）+siit
（2）*yk（i-2）-siit（3）*uk（i-
升级会员 