系统辨识最小二乘法大作业

1、系统辨识大作业最小二乘法及其相关估值方法应用Ff学院：自动化学院专业：信息工程学号：2007302171姓名：马志强日期：2010.11.14基于最小二乘法的多种系统辨识方法研究1. 最小二乘法的引出在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。设单输入一单输出线性定长系统的差分方程为x(k) + ax(k 1) anQc n) = bou(lc) + + bnu(k n), k = k, 2, 3,(5.1.1)式中:u(k)为随机干扰；x(k)为理论上的输出值。x(A)只有通过观测才能得到，在观测过程 中往往附加有随机干扰。x(A)的观测值y(幻可表示为y(k)=x(k) + n(

2、k)(5.1.2)式中:”(k)为随机干扰。由式(5.1.2)得x(k)=y(k)-n(k)(5.1.3)将式(5.1.3)带入式(5.1.1)得y(k) + "(A 一 1) + + %y(k - )n=bou(k) + b1u(k 1) + bnu(k -n) + n(k) + > q (k - i)i=l(5.1.4)我们可能不知道n(A)的统计特性，在这种情况卜，往往把7i(k)看做均值为0的白噪声。设i=l(5.1.5)则式(5.1.4)可写成y(k) = -aiy(k - 1) - a2y(k - 2)Ony(k n) + 如"(幻 + biu(jc 1)

3、 + + bnu(k-n) +(5.1.6)在观测u(k)时也有测危误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定 f(A)不仅包含了x(A)的测量误差，而且包含了u(幻的测最误差和系统内部噪声。假定f(幻是 不相关随机序列(实际上f(A)是相关随机序列)。现分别测出7i + N个随机输入ttiy(l), y(2),，y(n 4- N), u(l) , rz(2),，u(n + N),则可写成N个方程，即y(n + 1) = -aiy(n) - azy(n - 1)any(l) + bQu(n + 1) + biu(n) + + bnu(l)+ f (" + 1)y(n+

4、2) = -aLy(n + 1) - a2y(n)any(2) + bou (n + 2) + bLu(n + 1) + + bnu(2)+ f(" + 2)yn + N) = aiy(n + N 1) a2y(n + N 2)Ony(N) + bou(n + N)+ bi“(n+ N - 1) + bnU(N) +f(n + N)上述TV个方程可写成向最矩阵形式ry(n + i)idr.m+l)ly =y(” + 2).,0 =Qn 如，f =f(” + 2) y(n + N). bn -k(”+N)Jy(n + l)y(n + 2).y(n + N).y(”)-y(”+1)-y(

5、n + N k)D2)Ao rk zfk ? -y-y-y(u 2+ + m(nm u u u阳+1)f(”+2).油+3).(5.1.7)一 y(“)-y(” +1)yy(2)u(n + 1)u(n + 2)u(2)-y(n + N- 1)-yWu(n + N)则式(5.1.7)可写为y = <I>6 + 5(5.1.8)式中：y为N维输出向后：毛为N维噪声向量；。为(2n + l)维参数向鸵<!>为N x (2n + 1)测 量矩阵。因此式(5.1.8)是一个含有(2n+l)个未知参数，由N个方程组成的联立方程组。如 果N v 2n + 1,方程数少于未知数数目，则

6、方程组的解是不定的，不能唯一地确定参数向量。如果N = 2n+1,方程组正好与未知数数目相等，当噪声(=0时，就能准确地解出0 = <x»-1y(5.1.9)如果噪声毛壬0则0 =(5.1.10)从上式可以看出噪声玉对参数估计是有影响的，为了尽量较小噪声玉对0估值的影响。在给定 输出向量y和测量矩阵切的条件下求系统参数。的估值，这就是系统辨识问题。可用最小二乘 法来求。的估值，以下讨论最小二乘法估计。2. 最小二乘法估计算法设6表示。的最优估值，0表示y的最优估值，则有y = <!>§(5.1.11)欧+1)y=如 了), §)(n + N).写

7、出式(5.1.11)的某一行，则有y(fc)= 一甫y(k - 1) - aiyk - 2)ay(k-n) + bou(k) + biu(k - 1) + n+ 鼠u(k-n) + f(A) = - Z 吊y(k 一 j) + 2- i)，ki=l1=0=n + 1 f n + 2, ,+'(5.1.12)设。0)表示y(k)与0(/c)之差，即e(k) = yW-y(k)=ny。) -yW = y(k) > 名 y« -o +i=li=0(1 +ajz-1 + +盆)y(k) (b° + bLzL + + 扁 )u(k)=a(z-1)y(k) b(z-1)u

8、(k)» k = n + l,n + 2, , n + N(5.1.13)式中=14- azl HF az-71B(z-i)=扃 + bz-1 + .+ 慕厂"e(k)成为残差。把Sn+L n + 2,，n+N分别代入式(5.1.13)可得残差e(n +1), e(n + 2),，e(?i + N)。设e = e(n + l)e(n + 2) e(n + N)T则有e(k) = y y = y - <i>§(5.1.14)最小二乘估计要求残差的平方和为最小，即按照指数函数J = eTe = (y 4>0)，(y <X>0)(5.1.1

9、5)为最小来确定估值6,求/对0的偏导数并令其等于0可得= 2<1»丁3 «f>q)= 0(5.1.16)4>t00 = <I>Ty(5.1.17)由式(5.1.17)可得。的最小二乘估计0 = (<I>T0)-1<I>Ty(5.1.18)3 .递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。设己获得的观测数据长度为M 将式(5.1.8)中的y和毛分别用Yn，"孔来代替，即Yn =+ &(5.3.1)用0JY表示。的最小二乘估计，则(53.2)Pn= (n %)(53.5

10、)于是(53.6)如果再获得1组新的观测值u(“ + N + 1)和+ N + 1),则又增加1个方程yN+i =寸万+衣+切+i(53.7)式中)如+1 = y( + n +1)，切+i = f (” + n +1)传+i=y(n + N)y(N + 1) u(n 4- N + 1) u(N + 1)将式(5.3.1)和式(5.3.7)合并，并写成分块矩阵形式，可得比膈此(53.8)根据上式可得到新的参数估值=Pn+1 (的+ 如+1VV+1)(53.9)式中T t(% 6n+Wn+*Wn+i)根据矩阵求逆引理可以求得递推最小二乘法辨识公式®N+1 = M + KN+l(yN+l

11、Wn+；°N)(53.19)Kn+i = P n+iWn+i(1 + Wn+n'n+!)t(53.20) Pjv+l = Pn _ PnWn+1(1 +IPN(53.21) 由于进行递推计算需要给出箱和初值Po和&o,通过计算证明，可以取初值：0o = O. Po = C2h c是充分大的常数，为(2“ + l)(2n + l)单位矩阵，则经过若干次递推之后能够 得到较好的参数估计。3. 辅助变量法辅助变最法是一种可克服最小二乘有偏估计的一种方法，对于原辨识方程y = M + £(5.4.1)当杠k)是不相关随机序列时，最小二乘法可以得到参数向虽0的一致无偏

12、估计。但是，在实 际睥用中0)往往是相关随机序列。假定存在着一个(2n + 1) XN的矩阵Z满足约束条件lim= 0Nlim= Qn-8 n(542)式中Q是非奇异的。用Z乘以式(5.4.1)等号两边得ZTy = Zr00 + Zr(5.4.3)由上式得0 = (Z0Tz7y - (ZTfP)ZT(5.4.4)如果取0IV = (Z4>)Tzy(5.4.5) 作为。估值，则称Gw为辅助变量估值，矩阵z成为辅助变量短阵，z中的元素称为辅助变量。 常用的辅助变最法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。4. 广义最小二乘法广义最小二乘法是能克服最小二乘法有偏估计的另-种方法，这种

13、方法计算比较复杂 但效果比较好。下面直接介绍广义最小二乘法的计算步骤：应用得到的输入和输出数据u(k)和y(k)(k= 1, 2, 3,，n + N),按模型 Q(zT)y(A)=b(zT)u(k)+f(k)求出。的最小二乘估计§(1)=e (k) = 8 (zT)y(k) % (zT)tz(k)(3) 用残差e(幻代替f(幻，计算笋)=(nN)尸(口)W)(4) 计算宁(幻和五(幻y(1)W = y(k) + 樗)y(k 一 I) + . +智)y(k - m)W】)(k) = u(k) + 方(k - 1) + +曹u(k - m)(5) 应用得到的顼(幻和5(幻按模型a(z-1

14、)y伙)=b(z-1)u(k) 4- E(k)用最小二乘法重新估计。，得到。的第2次估值6。)。然后按步骤(2)计算残差e(2)(k),按步骤(3) 重新估计八得到估值再按照步骤(4)计算外气幻和五(2)(k),按照步骤(5)求。的第3 次估值6(3)。重复上述循环，之道。的估值6()收敛为止。5. 一种交替的广义最小二乘法求解技术(夏式法)这种方法是夏天长提出来的，又称夏式法。以上讨论过的广义最小二乘法的特点在于系 统的输入和输出信号反复过滤。一卜介绍的夏式法是一种交替的广义最小.乘法求解技术, 它不需要数据反岌过滤，因而计算效率较高。这种方法可消去最小二乘估计中的偏差，而且由这种方法导出的

15、计算方法也比较简单。基于以上的几种方法，有(y =中。+ (If = + e(5.7.1) 因而有y = <P0 + f2f + e = (P Q + e(5.7.2) 应用最小二乘法可得到参数估值(5.7.3) 可以推出0 =一 1 寺Ty (<j>T0)-10n/(5.7.11)上式中的第1项是。最小二乘估计岛s，第2项是偏差项6驴 所以必须准确计算为了准确计算0B，可采用迭代的方法。6. 专题解答设但输入-单输出系统的差分方程为y(A) = a、y(k 1) 。2贝* 2) + bu (k 1) + 2) + (k)f (A) = E(k) + QE(k 1) 4-。2

16、巴(* 2)取真实H0T = ai az b bz = 1.642 0.715 0.39 0.3引，输入数据如下所示k(k)ktl(k)ku(k)11.14711-0.958210.48520.201120.810221.6333-0.78713-0.044230.0434-1.159140.947241.3265-1.05215-1.474251.70660.86616-0.719260.34071.15217-0.086270.89081.57318-1.099281.14490.626191.450291.177100.433201.15130-0.390用。的真实值利用查分方程求出y(

17、k)作为测虽:值,£(k)为均值为0,方差为0.1, 0.5的不相关随机序列。(1) 用最小二乘法估计参数OT = a± az b± bzo(2) 用递推最小二乘法估计eT = ax a2如b2o(3) 用辅助变量法估计参数3T = aA a2知奶。(4) 设f(k) + 2f(k) = e(k)，用广义最小二乘法估计参数6T = ax也bi如。(5) 用夏式法估计参数OT = ax a2 妇。(6) 详细分析和比较所获得的参数辨识结果，并说明上述参数便是方法的优缺点。 根据题目要求的解法，利用Matlab编程实现系统辨识的估值利用最小二乘法估计的结果如卜.：最小

18、二乘法方差。20。20.00011.62801.63161.63541.63621.63601.62890.70280.70590.71200.70820.71650.70460.39710.39470.39180.39700.39060.39083.4490.34940.34630.35270.35340.34370.0011.15431.55771.60501.60601.61951.56700.67660.63710.68600.68160.70300.65720.40640.38680.37370.35830.39070.37520.33040.32490.32440.31670.33

19、660.31400.011.35380.89561.00081.34031.05741.12310.50100.16370.20360.47070.22890.29630.34860.42370.42360.38260.36820.35920.17090.06970.12680.26150.13170.15060.11.14241.02550.88960.81820.81000.77150.27100.17360.11090.11140.01530.13110.32840.37660.38930.42980.41220.47140.22160.18440.08130.09230.13520.0

20、6400.50.97510.89380.19270.55060.75600.94590.10170.07400.01970.03920.04940.13770.02710.34840.37620.65100.33720.38150.07920.16340.15210.05840.14370.1853部分程序运行结果9.W5 乜”“ 0.岫 0劣。A»lIt»t f.WJ O.bJbb 加。冲KbOW O.bM O.TTP iwn递推最小二乘法方差«2bi0.00011.67540.67870.52070.39241.30760.29000.00750.32841.

21、51460.69631.14010.16391.57330.77820.41490.71101.16020.47530.67360.31591.20910.31920.52770.02750.0011.47670.40400.26790.55121.62590.75940.32530.37801.53930.47570.12680.43461.15480.17000.19260.82510.88580.07600.33850.04061.41290.31270.09920.83800.011.34850.34450.31940.37101.16390.32960.78130.21621.994

22、61.23231.48520.03041.39240.35430.33190.45721.39820.36080.77730.31521.63460.72290.57800.39470.11.56240.71320.44220.41121.73350.71520.08440.63991.47630.53660.32550.31161.44770.34890.22180.22651.62160.70820.65950.42751.51050.40000.01130.22130.51.79270.94110.27300.34711.55560.88770.59720.12171.78681.253

23、81.12480.21001.57330.74340.35890.13871.31930.60841.29710.30291.59590.53860.01410.6947部分程序运行结果:* t% w< . wt flit ，t葬。仇“”小f.noi sowAh»A"辅助变量法方差«2bi。20.00011.77990.85880.41470.40461.30760.29000.00750.32841.67350.75780.40600.34841.58120.65460.37710.36111.66570.74690.37720.35611.52810.6

24、5090.36450.33020.0011.62950.67750.40820.35011.64250.73050.39370.35431.55950.60520.35630.33621.41450.49250.40210.28731.63710.72700.36810.34181.35390.47330.39060.24890.011.34510.47000.38220.26801.36570.48930.43490.24961.37020.50090.43880.26111.18840.37070.34570.15211.36360.53300.41350.20641.31580.4836

25、0.45620.21690.11.55450.61670.41040.35121.59000.66480.40980.36641.66100.70290.40010.34951.51040.60140.39680.32561.56200.64960.39500.30971.44180.57060.41830.27830.51.49520.57040.37690.34831.55920.65410.43870.33301.36370.53020.38650.19161.55430.63240.32080.29291.53850.59940.36610.35761.35110.48390.3412

26、0.2208部分程序运行结果:ge =曲*as MmL”。 O.g”03“A”广义最小二乘法方差。2如。20.00011.64510.71820.38950.35031.64680.72050.38660.34921.63080.70380.38520.34571.64670.71880.38980.35441.64110.71560.38980.34771.64410.71850.38960.34880.0011.70270.76780.40590.37611.64770.71660.39730.36431.63910.71130.38360.35041.65860.71660.37910.

27、36721.68310.74720.39940.37311.59640.68500.36550.32310.011.68260.74690.38590.36001.72450.76960.35650.36421.65770.71860.39100.38201.66560.72630.34250.34131.69030.73920.36330.38381.69420.74390.38340.38800.10.89960.15150.44350.06041.03800.30440.34100.04641.65710.65690.45890.42591.21800.45840.31520.10300

28、.19570.09740.21640.10351.68590.72320.41440.44040.51.40530.67030.36720.13721.57830.61480.35120.41731.02330.06580.20350.40601.59730.81030.35440.21311.61910.67090.31110.48661.08040.22100.75000.5185寸役*0.。*.州宙广义*4二豢用奈中，软角 ane -0.9309 0.07230.3039。.2 弛 2白0.。，、!剃国广小二星埔施喜於代' I anrt.3100.39S00.b2280S939&

29、quot;韦小0.000E：利可广二集冷: an$ -Ijno 0.39K0 0.622P 0.K939Ko.ooom.利句广二嗑鼻: am *赫，03 O.bbOb 。.67 奶 O.TT92 "部分程序运行结果:夏式法方差。20.00011.63670.709003910035111.64510.71840.38760.34721.36180.21590.82050.58511.62910.70750.39110.34641.61730.69290.39550.34391.63070.70660.39120.34490.0011.53140.64070.40160.30921.5

30、6320.65750.42040.33191.61720.69460.39250.32581.56700.64890.40160.32331.53260.62430.39950.30701.62510.69460.36210.33600.010.88150.10890.37340.04401.37630.49720.41210.26541.50270.59520.38540.32771.34210.42680.39780.27901.24540.30240.36690.25251.33610.47410.49990.29000.10.66390.10740.68470.08271.14060.

31、31540.42590.23670.17910.29160.94540.21261.20680.38270.74600.32200.64120.07040.58020.08771.03840.28780.74930.37930.50.93300.07520.97410.58820.62910.16860.46790.78690.59860.27920.26050.20080.44410.80330.74421.76670.62160.14400.58210.49060.65760.09991.00020.0263部分程序运行结果:. / 一A:，一 一tin T<we<。.颂 7

32、0.博25 O.g 0.99孙，o.nn o.m omor.t»W ».2U6 O.MK 9.b9Hnut.A «R«3930$勺仇网of二.网，i结论：通过编程计算，获得在噪声方差比较小的情况I、，各种方法所获得的估值比较理想, 但随着噪声方差的增大，估值的偏差随之增大，横向比较看来夏式法与广义最小二乘法能够 更好地还原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值。Mathb源程序：%最小二乘估计clearu= 1. 1470. 2010. 6260. 433 -0. 985-0.086 -1. 099 1.450 -0. 3400. 8

33、900. 144-0. 787 -1.589 -1. 052 0. 8661. 1521.5730. 810-0. 044 D. 947-1. 474 -0.7191. 151 0. 485 1. 6330. 0431. 3261. 7061. 177-0.390;n=normrnd(0, sqrt (0. 1), 1, 31); z=zeros (1, 30);for k=3:31 z(k)=-1.642*z(k-l)-0. 715*z(k-2)+0. 39*u(k-l)+0. 35*u(k-2)+n(k)+l. 642*n( k-l)+0.715*n(k-2);endh0=-z(2) -z

34、(l) u(2) u'HLT=hO, zeros (4, 28);for k=3:30hl=-z(k) -z(k-1) u(k) u(k-l)* ;HLT(:, k-D=hl;endHL=HLT'y=z(3) ;z(4) ;z(5) ;z(6) ;z(7) ;z(8) ;z(9) ;z(10) ;z(ll) ;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)

35、;%求出FAIcl=HL'*HL;c2=inv(cl);c3二HL'*y;c二c2*c3;最小二乘法估计辨识参数。如下：%display('方差=0. 1时，al=c(l);a2=c(2);bl=c(3);b2=c(4);clear%递推最小二乘法估计0. 201-0. 9851.4500. 144u=0. 626-0. 086-0. 3401. 1470.433-1. 0990. 890-0. 787 -1. 589 -1. 052 0. 8661. 1521. 5730. 810-0. 044 0. 947-1. 474 -0. 7191. 151 0. 485 1

36、. 6330. 0431. 3261. 7061. 177-0. 390;z(2)=0;z(l)=0;1, 31);n=normrnd (0, sqrt (0. 1), for k=3:31 z(k)=-l. 642*z(k-l)-0. 715*z(k-2)+0. 39*u(k-l)+0. 35*u(k-2)+n(k)+l. 642*n( k-l)+O.715*n(k-2);endcO=O. 001 0. 001 0. 001 0.001' %直接给出被辨识参数的初始值，即一个 充分小的实向量p0=10"6*eye(4, 4); %直接给出初始状态P0,即一个充分大的实数单位

37、矩阵 E=0. 000000005; %取相对误差E二0. 000000005c = c0, zeros (4, 30) ; %被辨U成矩阵的初始值及大小e=zeros(4, 30); %相对误差的初始值及大小for k=3:30; %开始求Khl = -z(k-l), -z(k-2), u(k-l), u(k-2)'x=hr *p0*hl+l; xl=inv(x); %开始求K(k)kl=p0*hl*xl；% 求出 K 的值dl=z(k)-hl'*c0; cl=c0+kl*dl; %求被辨识参数cel=cl-c0; %求参数当前值与上一次的值的差值e2=el. /c0; %求

38、参数的相对变化e(:, k)=e2; %把当前相对变化的列向带加入误差矩阵的最后一列c0=cl; %新获得的参数作为下一次递推的旧参数c(:,k)=cl； %把辨识参数c列向量加入辨识矩阵的最后一列pl=pO_kl*kl，*hr *pO*hl+l; %求出 p(k)的值p0=pl: %给下次用if e2<=Ebreak; %如果参数收敛情况满足要求，终止计算endend%display C方差为0. 0001递推最小二乘法辨识后的结果是：');al=c(l,：);a2=c(2,:);bl=c(3,:);b2=c(4,:);%displayC al, a2, bl, b2经过递推最

39、小二乘法辨识的结果是：');for i=3:31;if(c(l, i)=0)ql=c(l, i-1)；break;endendfor i=3:31;if(c(2, i)=0)q2=c(2, i-1);break;end endfor 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;endendal=ql;a2=q2;bl=q3 ；b2=q4;%clear%辅助变量递推最小二乘法估计na=2;nb=2;siitt=1.642 0.715 0. 39

40、si it 0=0. 001*eye (na+nb, 1);p=10 6*eye(na+nb);siit (:, 1) =siitO;y(2)=0;y(l)=0;x(l)=0;x(2)=0;j=0;u=0. 626-0. 086-0. 3401. 1470.433-1. 0990. 8900. 201-0. 9851.4500. 144-0. 7870. 8101. 151 0. 485 1. 6331. 177-0. 390;-1. 589 -1. 052 0. 8661. 1521.573-0. 044 0. 947-1. 474 -0.7190. 0431. 3261. 706n=nor

41、mrnd(0, sqrt (0. 01), 1, 31);for k=3：31;h=-y(k-1), -y(k-2), u(k-l), u(k-2)'y(k)二h *siitt+n(k)+l. 642*n(kl)+0. 715*n(k2);hx=-x(k-l), -x(k-2), u(kl), u(k2)J，;kk=p*hx/ (h' *p*hx+l);p= eye (na+nb) -kk*h' *p;siit (:, kl) =siitO+kk*y(k) -h' *siitO;x(k)二hx' *siit (:, kl);j=j+(y(k)-h'

42、;*l. 642 0.715 0.39 0.35')2e=max(abs(siit (:, k-1)-siitO). /siitO);ess(:, k2) =siit (:, kl) -siitt;siitO=siit (:, kl);endal=siitO(l);a2=siitO(2);bl=siit0(3);b2=siit0(4);clear%广义最小二乘估计clear;nn = normrnd(0, sqrt (0. 5), 1, 31)'uk=l. 147 0. 201 -0. 787 -1. 589 -1. 052 0. 866 1. 152 1. 573 0. 62

43、6 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. 706 -0. 340 0. 890 1. 144 1. 177 -0. 390;yk(l)=0；yk(2)=0;for i=l:29;yk(i+2)=-1. 642*yk(i+l)-0. 715*yk(i)+0. 39*uk(i+l)+0. 35*uk(i)+nn(i+2)+l.642*nn(i+l)+0. 715*nn(i);end;for i=l:29;A(i, :

44、) = -yk(i+l) -yk(i) uk(i+l) uk(i);endsiit=inv(A，*A)*A' *(yk(3:31)+nn(2:30)')'e(l)=yk(l);e(2) =yk(2)+siit (l)*yk(l)-siit (3)*uk(l);for i=3:31;e(i) =yk(i)+siit(l)*yk(il)+siit(2)*yk(i2) -siit(3)*uk(il)-siit(4)*uk (i-2);endfor i=l:29;fai (i, :) = -e(i+l) -e(i);endf=inv(fai，*fai)*fai* *e (3:3

45、1)?;for i=3:31;yk(i) =yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);endyk(2) =yk +f(l)*yk ;for i=3:30;uk=uk(i)+f(l)*uk(i-l)+f(2)*uk(i-2);enduk(2)=uk(2)+f(l)*uk(l);for j=l:30 for i=l:29;A(i, :) = -yk(i+l) -yk(i) uk(i+l) uk(i);endsiit=inv(A'*A)*A'*yk(3:31)'e(l)=yk(l);e(2) =yk(2)+siit (l)*(yk(l) -siit (3

46、)*uk(l); for i=3:31;e(i) =yk(i)+siit (l)*(yk(iT) )+siit (2)*(yk(i2) -siit (3)*uk(iT)-siit (4 )*uk(i-2);endfor i=l：29;fai(i, ：) = -e(i+l) -e(i);endf=inv(fai'*fai)*fai'*e(3:31)'kl(j)=f(l);k2(j)=f(2);for i=3:31;yk(i) =yk(i)+f *(yk(iT)+f (2)*(yk(i-2):endyk(2)=yk(2)+f(l)*yk(l);for i=3:30uk(i) =uk(i)+f (l)*uk(i.T)+f (2)*uk(i-2);enduk(2)=uk(2)+f(l)*uk(l);endsiit'clear;u= 1. 1470. 2010. 6260. 433 -0. 985-0.086 -1. 099 1.450 -0. 3400. 8900. 144%用夏氏偏差修正法估计参数-0. 787 -1. 589 -1. 052 0. 8661. 1521. 5730. 810-0. 04