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

系统辨识大作业

最小二乘法及其相关估值方法应用

学号：姓名：王家琦

基于最小二乘法的多种系统辨识方法研究

1. 最小二乘法的引出

在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。设单输入-单输出线性定长系统的差分方程为

x k +a1x k−1 +⋯+an k−n =b0u k +⋯+bnu k−n ，k=1，2，3，⋯ (5.1.1)

式中: u k 为随机干扰；x k 为理论上的输出值。x k 只有通过观测才能得到，在观测过程中往往附加有随机干扰。x k 的观测值y k 可表示为

y k =x k +n k

(5.1.2) 式中: n k 为随机干扰。由式(5.1.2)得

x k =y k −n k

(5.1.3) 将式(5.1.3)带入式(5.1.1)得

y k +a1y k−1 +⋯+any k−n

=b0u k +b1u k−1 +⋯+bnu k−n +n k + ai(k−i)

i=1

(5.1.4)

我们可能不知道n k 的统计特性，在这种情况下，往往把n k 看做均值为0的白噪声。设

ξ(k)=n(k)+ ai(k−i)

i=1

(5.1.5) 则式(5.1.4)可写成

y k =−a1y k−1 −a2y k−2 −⋯−any k−n +b0u k +b1u k−1 +⋯

+bnu k−n +ξ(k)

(5.1.6) 在观测u k 时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定ξ(k) 不仅包含了x k 的测量误差，而且包含了u k 的测量误差和系统内部噪声。假定ξ(k) 是不相关随机序列(实际上ξ(k) 是相关随机序列) 。现分别测出n+N个随机输入值y 1 ，y 2 ，⋯，y n+N ，u 1 ，u 2 ，⋯，u(n+N) ，则可写成N个方程，即

y n+1 =−a1y n −a2y n−1 −⋯−any 1 +b0u n+1 +b1u n +⋯+bnu 1

+ξ(n+1)

y n+2 =−a1y n+1 −a2y n −⋯−any 2 +b0u n+2 +b1u n+1 +⋯+bnu 2

+ξ(n+2)

⋮

y n+N =−a1y n+N−1 −a2y n+N−2 −⋯−any N +b0u n+N

+b1u n+N−1 +⋯+bnu N +ξ(n+N)

上述N个方程可写成向量-矩阵形式

⋯−y(1)u(n+1) ⋯−y(n) y(n+1)

⋯−y(2)u(n+2) ⋯−y(n+1) y(n+2)

⋮⋮⋮ ⋮ ⋮⋮y(n+N) −y(n+N−1) ⋯−y(N) u(n+N) ⋯ 设

a1 ⋮ u(1)ξ(n+1) an u(2)ξ(n+2) × b0 +

⋮

⋮ ξ(n+3) u(N)

(5.1.7)

a1 ⋮ y(n+1) ξ(n+1) an y(n+2) ξ(n+2)

y= ，θ= b0 ，ξ=

⋮⋮

⋮ y(n+N) ξ(n+N) bn ⋯−y(1)u(n+1) ⋯u(1)−y(n)

⋯−y(2)u(n+2) ⋯u(2)−y(n+1)

Φ=

⋮⋮ ⋮ ⋮⋮

−y(n+N−1) ⋯−y(N) u(n+N) ⋯u(N)

则式(5.1.7)可写为

y=Φθ+ξ

(5.1.8)

式中：y为N维输出向量；ξ为N维噪声向量；θ为(2n+1) 维参数向量；Φ为N×(2n+1) 测量矩阵。因此式(5.1.8)是一个含有(2n+1) 个未知参数，由N个方程组成的联立方程组。如果N

θ=Φ−1y

(5.1.9) 如果噪声ξ≠0，则

θ=Φ−1y−Φ−1ξ (5.1.10)

从上式可以看出噪声ξ对参数估计是有影响的，为了尽量较小噪声ξ对θ估值的影响。在给定输出向量y和测量矩阵Φ的条件下求系统参数θ的估值，这就是系统辨识问题。可用最小二乘法来求θ的估值，以下讨论最小二乘法估计。 2. 最小二乘法估计算法

表示 θ的最优估值，y 表示 y的最优估值，则有设θ

=Φθy

(5.1.11)

a 1 ⋮ y (n+1)

a ny (n+2) = y ，θ=

⋮ b0 y (n+N) ⋮

写出式(5.1.11)的某一行，则有

y k =−a y y y1 k−1 −a2 k−2 −⋯−an k−n +b0u k +b1u k−1 +⋯

iu(k−i) ，k+b iy(k−i) + bnu k−n +ξ k =− a

i=1

i=0

=n+1，n+2，⋯，n+N

设e(k) 表示y(k) 与y (k) 之差，即 e k = y k -y k =

(5.1.12)

iu k−i = y k −y k =y k — a iy k−i + b

i=1

i=0

−1−n −n −1+⋯+b 1+a z+⋯+a zyk− b u k 1n0+b1znz

z−1 u k ，k=n+1，n+2，⋯，n+N =a z−1 y k −b

(5.1.13)

式中

−1−n

a z−1 =1+a z+⋯+a z 1n

−n z−1 =b −1+⋯+b b 0+b1znz

e k 成为残差。把k=n+1，n+2，⋯，n+N分别代入式(5.1.13)可得残差e n+

1 ，e n+2 ，⋯，e n+N 。设

e=[e n+1 e n+2 ⋯e n+N ]T 则有

=y−Φθe k = y−y

最小二乘估计要求残差的平方和为最小，即按照指数函数

) T(y−Φθ ) J=eTe=(y−Φθ

。求J对θ 的偏导数并令其等于0可得为最小来确定估值θ

∂J

=0 =−2ΦT y−Φθ

∂θ

(5.1.16) (5.1.17)

=ΦTy ΦTΦθ =(ΦTΦ) −1ΦT y θ

(5.1.14)

(5.1.15)

由式(5.1.17)可得θ的最小二乘估计

(5.1.18)

3. 递推最小二乘法

为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。设已获得的观测数据长度为N，将式(5.1.8)中的y和ξ分别用Y N ，ΦN， ξN来代替，即

N Y N =ΦNθ+ξ

(5.3.1)

N表示 θ的最小二乘估计，则用θ

N=(ΦNTΦN) −1ΦNTYN θ

(5.3.2)

设

PN=（ΦNΦN）

(5.3.5)

于是

N=PNΦNTYN θ

(5.3.6) 如果再获得1组新的观测值u n+N+1 和y(n+N+1) ，则又增加1个方程

yN+1=ψTN+1θ+ξN+1

(5.3.7)

−1

式中

yN+1=y n+N+1 ，ξN+1=ξ(n+N+1)

ψTN+1=[−y n+N ⋯ −y N+1 u n+N+1 ⋯ u(N+1)]

将式(5.3.1)和式(5.3.7)合并，并写成分块矩阵形式，可得

ΦNYξ

N = T θ+ N yN+1ψN+1ξN+1

−1

(5.3.8)

根据上式可得到新的参数估值

ΦNTΦN

N+1= T T θ

ψN+1ψN+1

式中

ΦNTYNΦNTYN T =PN+1 T ψN+1yN+1ψN+1yN+1

(5.3.9)

=PN+1(ΦNTYN+ψN+1yN+1)

−1

ΦNTΦN

PN+1= T T

ψN+1ψN+1

=（ΦNΦN+ψN+1ψN+1）

−1

根据矩阵求逆引理可以求得递推最小二乘法辨识公式

N+1=θ N+KN+1(yN+1−ψN+1Tθ N) θ

(5.3.19) (5.3.20)

KN+1=PN+1ψN+1(1+ψN+1TPNψN+1) −1

PN+1=PN−PNψN+1(1+ψN+1TPNψN+1) −1ψN+1TPN

(5.3.21)

N和PN初值P0和θ 0，通过计算证明，可以取初值：θ 0=0，由于进行递推计算需要给出θ

P0=c2I，c 是充分大的常数，I为 2n+1 (2n+1) 单位矩阵，则经过若干次递推之后能够得到较好的参数估计。 3. 辅助变量法

辅助变量法是一种可克服最小二乘有偏估计的一种方法，对于原辨识方程

y=Φθ+ξ

(5.4.1) 当ξ(k) 是不相关随机序列时，最小二乘法可以得到参数向量θ的一致无偏估计。但是，在实际应用中ξ(k) 往往是相关随机序列。

假定存在着一个(2n+1) ×N的矩阵Z满足约束条件

Zξ=0N→∞

limZTΦ=QN→∞lim

(5.4.2)

式中Q是非奇异的。用ZT乘以式(5.4.1)等号两边得

ZTy=ZTΦθ+ZT ξ

由上式得如果取

θ=(ZTΦ) −1ZTy−(ZTΦ) −1ZTξ

IV =(ZTΦ) −1ZTy θ

(5.4.4)

(5.4.3)

(5.4.5)

IV 为辅助变量估值，作为θ估值，则称θ矩阵Z成为辅助变量矩阵，Z中的元素称为辅助变量。

常用的辅助变量法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。 4. 广义最小二乘法

广义最小二乘法是能克服最小二乘法有偏估计的另一种方法，这种方法计算比较复杂但效果比较好。下面直接介绍广义最小二乘法的计算步骤： (1)应用得到的输入和输出数据u(k) 和y k (k=1，2，3，⋯，n+N) ，按模型

a z−1 y k =b z−1 u k +ξ k

求出θ 的最小二乘估计

(1)a 1 ⋮ (1) a n

= (1) b 0 ⋮(1) b n

(1) θ

(2)计算残差e 1 (k)

(1) z−1 u k e 1 k =a (1) z−1 y k −b

(1) =[(Ω(1) ) TΩ(1) ]−1(Ω(1) ) Te 1 (3)用残差e 1 k 代替 ξ k ，计算f

(4)计算y 1 (k) 和u 1 (k)

1 1 y k−1 +⋯+f y 1 k =y k +fmy(k−m) 1

1 1 u k−1 +⋯+f u 1 k =u k +fmu(k−m) (5)应用得到的y

(k) 和u (k) 按模型

a z−1 y 1 k =b z−1 u 1 k +ε(k)

(2) 。然后按步骤(2)计算残差e 2 (k) ，按步骤用最小二乘法重新估计θ，得到θ的第2次估值θ

2 。再按照步骤(4)计算y(3)重新估计f，得到估值f 2 (k) 和u 2 k ，按照步骤(5)求θ的第3 (3) 。重复上述循环，之道θ的估值θ (i) 收敛为止。次估值θ

5. 一种交替的广义最小二乘法求解技术(夏式法)

这种方法是夏天长提出来的，又称夏式法。以上讨论过的广义最小二乘法的特点在于系统的输入和输出信号反复过滤。一下介绍的夏式法是一种交替的广义最小二乘法求解技术，

它不需要数据反复过滤，因而计算效率较高。这种方法可消去最小二乘估计中的偏差，而且由这种方法导出的计算方法也比较简单。

基于以上的几种方法，有

y=Φθ+ξ ξ=Ωf+ε

因而有

(5.7.1)

y =Φθ+Ωf+ε= ΦΩ +ε

应用最小二乘法可得到参数估值

−1TT θΦ = T ΦΩ ΦT y

fΩΩ

可以推出

(5.7.2)

(5.7.3)

=(ΦTΦ) −1ΦT y−(ΦTΦ) −1ΦΩfθ

(5.7.11)

LS，第2项是偏差项θ B，所以必须准确计算θ B。上式中的第1项是θ最小二乘估计θ

B，可采用迭代的方法。为了准确计算θ6. 专题解答

设但输入-单输出系统的差分方程为

y k =−a1y k−1 −a2y k−2 +b1u k−1 +b2u k−2 +ξ(k)

ξ k =ε k +a1ε k−1 +a2ε k−2

取真实值θT

用θ的真实值利用查分方程求出yk作为测量值，εk为均值为0，方差为0.1，0.5的不相关随机序列。

(1) 用最小二乘法估计参数θT= a1a2b1b2 。 (2) 用递推最小二乘法估计θT= a1a2b1b2 。 (3) 用辅助变量法估计参数θT= a1a2b1b2 。

(4) 设ξ k +2ξ k =ε k ，用广义最小二乘法估计参数θT= a1a2b1b2 。 (5) 用夏式法估计参数θT= a1a2b1b2 。

(6) 详细分析和比较所获得的参数辨识结果，并说明上述参数便是方法的优缺点。根据题目要求的解法，利用Matlab 编程实现系统辨识的估值

利用最小二乘法估计的结果如下：

部分程序运行结果

部分程序运行结果：

结论：通过编程计算，获得在噪声方差比较小的情况下，各种方法所获得的估值比较理想，但随着噪声方差的增大，估值的偏差随之增大，横向比较看来夏式法与广义最小二乘法能够

更好地还原参数值，当观测值足够多时，各种方法都能很好地反映参数真实值。 Matlab 源程序：

%最小二乘估计% 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.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 k=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)]; for k=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

%递推最小二乘法估计

n=normrnd(0, sqrt(0.1), 1, 31);

for k=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.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值, 即一个充分小的实向量

p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵 E=0.000000005; %取相对误差E=0.000000005

c=[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; %求被辨识参数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; %给下次用 if e2

break ; %如果参数收敛情况满足要求，终止计算 end end

%display('方差为0.0001递推最小二乘法辨识后的结果是：'); a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);

%display('a1,a2,b1,b2经过递推最小二乘法辨识的结果是：'); for i=3:31; if (c(1,i)==0)

q1=c(1,i-1); break ; end end

for i=3:31; if (c(2,i)==0)

q2=c(2,i-1); break ; end

end

for i=3:31; if (c(3,i)==0)

q3=c(3,i-1); break ; end end

for i=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.642 0.715 0.39 0.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;

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-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.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.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)]; 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); 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); end

for i=1:29;

fai(i,:)=[-e(i+1) -e(i)]; end

f=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); end

yk(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); end

uk(2)=uk(2)+f(1)*uk(1); for j=1:30

for i=1:29;

A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)]; end

siit=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); end

for i=1:29;

fai(i,:)=[-e(i+1) -e(i)]; end

f=inv(fai'*fai)*fai'*e(3:31)'; k1(j)=f(1); k2(j)=f(2); for i=3:31;

yk(i)=yk(i)+f(1)*(yk(i-1))+f(2)*(yk(i-2)); end

yk(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); end

uk(2)=uk(2)+f(1)*uk(1); end

siit';

%用夏氏偏差修正法估计参数 clear;

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)]; for k=3:30

h1=[-z(k) -z(k-1) u(k) u(k-1)]';

HLT(:,k-1)=h1; end

HL=HLT';

fai1_xiashi=HL; siit1_ls_guzhi=c;

ka_ma=inv(fai1_xiashi'*fai1_xiashi)*fai1_xiashi'; M_xiashi=eye(29)-fai1_xiashi*ka_ma; siit_N0_xshi=siit1_ls_guzhi; siit_bxs=[10,10,10,10]'; for n=1:30

e=zeros(29,1);

e(3)=y(5)+siit_N0_xshi(1,1)*y(4)+siit_N0_xshi(2,1)*y(3)-siit_N0_xshi(3,1)*u(4)-siit_N0_xshi(4,1)*u(3);

e(4)=y(6)+siit_N0_xshi(1,1)*y(5)+siit_N0_xshi(2,1)*y(4)-siit_N0_xshi(3,1)*u(5)-siit_N0_xshi(4,1)*u(4);

e(5:29) = y(5:29) - fai1_xiashi(5:29,:) * siit_N0_xshi; ou_mei_ka=zeros(29,1); for i=1:29

ou_mei_ka(i,:)=-e(i); end

D_xiashi=ou_mei_ka'*M_xiashi*ou_mei_ka;

f_xiashi=inv(D_xiashi)*ou_mei_ka'*M_xiashi*y(1:29); siit_b=ka_ma*ou_mei_ka*f_xiashi; siit_xiashi=siit1_ls_guzhi-siit_b;

if norm(siit_b-siit_bxs)/norm(siit_bxs)

siit_N0_xshi=siit_xiashi; siit_bxs=siit_b; end

siit1_xiashi=siit_xiashi;

%display('利用夏氏修正法得出的辨识结果：'); a1=siit1_xiashi(1);

a2=siit1_xiashi(2); b1=siit1_xiashi(3); b2=siit1_xiashi(4);

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

相关内容

热门内容