小二乘法在系统辨识中的应用包含相关的三种算法

1、2008级硕士研究生系统建模理论试卷已知一个三阶线性离散系统的输入.输出数据，共有40个采样值，试分别 用：最小二乘法（LS）、递推最小二乘法（RLS）、广义最小二乘法（GLS）进行 参数估计，给出相应的数学模型，并阐述相应的辨识原理。k 12345678910uO) 0. 82510. 09880. 4628-0.9168 2.2325 0.0777 2.36540. 34761. 1473-1. 9035y(力 1.5333-1. 06801. 0666-0. 5284 -0. 5835 3.1471 -3 71856. 2149-6. 30267. 2705k 1112131415161

2、7181920u(A)0.92291. 64000.84100. 7599-0. 4739 -0.17847760-1. 67221. 2959-0.0591 y -9.05521. 07868.1735-5. 90043. 98702. 2486 0. 9525 0. 5325-1. 52270. 4200k 21222324252627282930u -1 0576 一1 00711.1342-0. 07400. 6759 0. 5221 0. 99540. 5271 -1. 76560. 4936y -1. 5579 0. 6640 -1. 4222 2. 6444 -2. 9572 3

3、. 6340 -3.1281 3. 8334 -3. 25421.1568402.11455. 1507k 313233343536373839u() 1-4810 0 9591 -3.1293 -0.3604 -0.4251 0.4185 -0. 6728 -0. 0027 1.1157y(A) 0. 0615 0. 9120 -0. 0692 -3. 2731 3 7486 -4 3194 4. 7230 -5. 2781 -2. 7235一-最小二乘法（LS）1、数学模型设时不变SIS0动态过程的数学模型为A（z）z 伙）=伙）+ “伙）（1.1.1）其中，“伙）为过程的输入量，Z伙）为

4、过程的输岀量，n伙）是噪声，多项式A（z-*）和B（z）为：J A(z J) = 1 + z + (“z J + + % 广B(zl) = btzl + b2z2 + + bbznb在本题中,na=nb=3.即A(zi) = + alzl +a2z2 + Bb+b+b.z3将此模型写成最小二乘格式z伙)= h伙)8 +”伙)(1.1.2)其中，讹)是过程的输岀量：h伙)是可观测的数据向量；”是均值为零的随机噪声。Jh伙)=_祕_1),一2伙_2)厂2伙_3)川(_1),伙_2),“伙_3) 式中S8=4，2皿3，勺厶对于k = 2 L，方程式(1.1.2)构成一个线性方程组，可以把它写成勺=z

5、(l),z(2),z(/)r q =(1)/(2)丿(/)h产_hr(l)hr(2)_ - z(0)71)0-z(0)00“(0)H(l)0“(0)0 _0-2(/-2)-z(/-3)临一 2)w(/-3)_(1.1.3)利用数据序列滋小和小伙),极小化准则函数丿(8) = z伙)-h(切OF(1.1.4)使J(6) = min的6估计值记作6,称为参数的最小二乘估计值。通过极小化(1.1.4)式 来讣算6的方法称作最小二乘法，未知模型参数最可能的值是在实际观测值与计算值之累次误差的平方和达到最小处所得 到的，这种模型输岀能最好地接近实际过程的输出。2、辨识原理考虑模型(1.1.2)式的辨识问

6、题，其中択)和h伙)都是可观测的数据，是待估计参数, 准则函数取(1.1.4)根据(1.1.3)的泄义，准则函数丿(0)可写成二次型的形式J(e)= (z/-H/(9)r(z/-H,0)(1.2.1)显然上式中的Hf代表模型的输出，或者说是过程的输出预报值。因此丿(0)可以看作来衡量模型输出与实际过程输岀的接近情况。极小化7(0),求得参数。的估计值e = (H；HZ 尸 H；Z/(1.2.2)将使模型的输出最好的预报过程的输岀。3、辨识结果6 = q a2 a3 /?, b2 仇=卜 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000二、递推最小二乘法(RLS

7、)1、数学模型在第一部分中建立了最小二乘法，并用一次完成算法进行了计算。但是由于具体使用时占用 内存量大，而且不能用于在线辨识。所以引入了最小二乘参数估计的递推算法。递推算法的基本思想可以概括成：新的估计值&伙)二老的估计值0(k -1) +修正项(2.1.1)在此算法中，时不变SIS0动态过程的数学模型仍与最小二乘法的一样。模型的最小二乘格 式也相同，只是计算方法不同，具体计算方法，在递推最小二乘法的辨识原理中描述。2、辨识原理首先将1. 2. 2式一次完成算法写成6 = (H；HH；Z 2P(/)H；石=h(i)hp)“f h(讥(2 2 1定义h 伙)= h(讪七)2H：H匕(2.2.2

8、)PT- l) =ha)h()2H；TH /=ihr(l)一护_其中：=hhr(k)h伙 一1)由(222)式可得：k-P伙)=工h(,)h + h伙)h仗)=P伙-l) + h伙)h伙)(2. 2.3)令:z=z(l), z,z伙一 1)则：6（A:-l） = （H_1HA_1）-,HL1zi_1 =P（k l）fhzJ-l r于是有 P伙 一 1）6伙 一 1）=工h（7）z（，）（2. 2.4）r-l令 zk = z（i）, z,z伙）r利用（2.2.3）和（2.2.4）式，可得6伙）= P 伙）h（力r-l=P伙）P伙1）6仗 一 1） + h伙比伙）（2. 2.5）=p（r川p伙）h

9、伙）伙）6伙1）+ h伙比伙）=6伙一 1）+p伙）h伙）乙伙）-hr伙）6伙1）引进增益矩阵K伙），泄义为K伙）=P伙）h 伙）（2.2.6）则（2.2.5）式写成6伙）=6伙一） + K（k）z（k）-hr（kWk-）（2. 2. 7）进一步把（2.2.3）式写成P（k） = IP伙 -1） + h 伙）hr 伙）（2. 2.8）为了避免矩阵求逆运算，利用矩阵反演公式可将（2.2.8）式演变成P伙）=P-1 伙-1） + h伙）h伙）=I-K伙）h伙）P伙-1）（2.2.9）将（2.2.9）式代入（2.2.6）式，整理后有K伙）=P伙 一 l）h伙）If 伙）P伙 一 l）h伙）+1-1（

10、2. 2.10）综合（2.2.7）、（2.2.9）、（2. 2. 10）式便得到最小二乘参数估计递推算法。6伙）=6伙-1） + K伙）z伙）-h 伙）6伙 一 1）P伙）=p-1 伙 _ 1） + h伙）h伙）“ =I-K伙）伙）P伙 一 1）K伙）=P 伙-l）h 伙）伙）P 伙-l）h伙）+1-1利用上述公式即可求得参数的估计值63、辨识结果6 = q a2 他 S b2 如=卜 1.0005 - 0.0007 - 0.0002 0.0000 - 0.0005 - 0.0001令纟伙)=貢巧叩0(3.2.5)三.广义最小二乘法(GLS)1、数学模型设时不变SISO动态过程的数学模型为法严

11、)(3.L1)其中伙)和Z伙)分别表示过程的输入和输出，咻)是均值为零的不相关随机噪声,多 项式A(z)、B(rl)和C(z)为：/l(z_1) = l+1z1 +6/2z2 +- + ”/-* B(zl) = blZ +b2z2 + + bisz：,bC(z-,) = l + cIz-,+c2z-2+- + crtrz在本题中,在本题中,na=nb=ne=3. R卩A(z_i) = 1+6/1z1 +a2z2 +d3z (厂)=忆+山严+仅z八、!“_3(312)c(z ) = 1 + CZ +C2z +C3z 广义最小二乘发的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法

12、 对滤波后的数据进行辨识。2、辨识原理由(3.1.1)式可得 A(z)C(z-*)z(k) = Bz)C(z1 )u(k) + v(k)(3. 2.1)(3. 2.2)今 k/W = c(r1k 伙) 、L仃伙)= C(L)伙)(3. 2.3)爻 ()=勺伙一一 2),兮伙 一 3),“/ 伙 1),“ f (k 2),u f ( 3)0 = pa2,31,/?2,Z?3r(3. 2.4)可将模型(3.1.1)式化成最小二乘格式Zf(k) = h/(k)Q + v(k) 由于上式v伙)是白噪声，所以利用最小二乘法即可获得参数的无偏估计。但是，数据向 量h/伙)中的变量均需要按照(322)式计算

13、，然而噪声模型C(z)并不知道。为此需要用迭代的方法来估计C)。中血伙）= Y伙一l）r伙一2）,锹一3）r心0 I直 S（3. 2. 6）就把噪声模型3. 2.5）也化成最小二乘格式0伙）= h：伙）+卩伙）由于上式的噪声已是白噪声，所以再次利用最小二乘法可获得噪声模型参数$的无偏估讣。 但是，数据向星h伙）依然包含着不可测得噪声量e伙-1）伙-2）,0伙-3）,它可用相应的估计值代替，置h伙）=一0伙一 1）,一&伙一 2）,-&伙一 3）（3. 2. 7）其中e（k） = 0,k4 时for i=4:1:40zf (i) =z (i) +z (i-1) *oe (il, l)+z(i2)

14、*oe(i-l, 2) +z(i-3)*oe(i-l, 3); uf (i) =u(i) +u(i-1)*oe(il, l)+u(i-2)*oe(i-l, 2) +u(i-3)*oe(i-l, 3);hf(i, :) = -zf(i) -zf(i-l) -zf(i2) uf(i) uf(il) uf(i-2);kf (i, :) = pf (:, :, i-l)*hf (i, :), *inv(hf (i, :)*pf (:, :, il)*hf (i, :) +1);pf (:, :, i) = eye(6)-kf (i, :), *hf (i, :)*pf (:, :, i-1);of (i, :) = of (i-1, :), +kf (i, :), *z(i)-hf (i, :)*of (i-1, :); e(i)=z(i)-hf(i, :)*of(i,:),;he(i, :) = -e(i-l) -e(i2) -e(i-3);ke(i, :) = pe(:, :, i-l)*he(i, :), *inv(he(i