应用最小二乘一次完成法和递推最小二乘法算法的系统辨识

1最小二乘法的理论基础

1.1最小二乘法

设单输入单输出线性定长系统的差分方程表示为：

其中δ(k)为服从N(0,1)的随机噪声,现分别测出n+N 个输出输入值y(1),y(2),…,y(n+N),u(1),u(2),…,u(n+N),则可写出N 个方程，写成向量－矩阵形式

(4.1.1)

()()()()()()()()

1201121n n y k a y k a y k a y k n b u k b u k b u k n k ξ=-------+

+-+

+-+()()()()()()101122,,n n a y n n y n a n y b y n N n N b ξξθξξ??

??++????????????++?

???===??????????????++??????????

????

()()()()()()()()()

()

()()()()

()(

)()()10111121222112n n y n y n y u n u y n y n y u n u y n N y n N y N u n N u N a n a n b n N b ξξξ+--+????

????+-+-+????

=?????????+-+--+????

??

????

??+??

??????+??+???????

???+??????????

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

11y θφφξ--=-

在给定输出向量y 和测量矩阵Φ的条件下求参数θ的估计，这就是系统辨识问题。 设 表示 θ 的估计值，?表示y 的最优估计，则有 （4.1.3） 式中：

()()()10??1??2??,???n n a

y n a y n y b y n N b θ????

+????????+????==????????+??????

????

设e(k)=y(k)- ?(k), e(k)称为残差，则有e=y- ?=y-Φθ 最小二乘估计要求残差的平方和最小，即按照指数函数

（4.1.4）

求Ｊ对 的偏导数并令其等于0可得：

（4.1.5）

由式（4.1.5）可得的 θ 最小二乘估计：

（4.1.6）

J 为极小值的充分条件是：

即矩阵ΦT Φ为正定矩阵，或者说是非奇异的。 1.1.1最小二乘法估计中的输入信号

当矩阵ΦT Φ的逆阵存在是，式（1.1.6）才有解。一般地，如果u(k)是随机序列或伪随机二位式序列，则矩阵ΦT Φ是非奇异的，即（ΦT Φ）－1存在，式（1.1.6）有解。 现在从ΦT Φ必须正定出发，讨论对u(k)的要求。

y φθξ=+?θ??y

θ=Φ()()

??T

T J e e y y θ

θ==-Φ-Φ?θ()

?20?T J y θ

θ

?=-Φ-Φ=??T T y θ

ΦΦ=Φ()1

?T T y θ

-=ΦΦΦ220?T

J θ

?=ΦΦ>?1

n N yy

yu T

+-ΦΦ??

当N 足够大时有

（4.1.8）

如果矩阵ΦT Φ正定，则Ru 是是对称矩阵，并且各阶主子式的行列式为正。当N 足够大时，矩阵Ru 才是是对称的。

由此引出矩阵ΦT Φ正定的必要条件是u(k)为持续激励信号。如果序列{u(k)}的n+1阶方阵Ru 是正定的，则称{u(k)}的n+1阶持续激励信号。 下列随机信号都能满足Ru 正定的要求 1. 有色随机信号 2. 伪随机信号 3.

白噪声序列

1.1.2最小二乘估计的概率性质 最小二乘估计的概率性质 1） 无偏性

由于输出值y 是随机的，所以θ是随机的，但注意θ不是随机值。

如果{}

{}?E E θ

θθ==，则称?θ是θ无偏估计 2）一致性

估计误差θ的方差矩阵为

（4.1.9）

上式表明当N →∞时，θ以概率1趋近于θ。

当ξ(k )为不相关随机序列时，最小二乘估计具有无偏性和一致性 3)有效性

如果式（1.1.2）中的ξ的均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值 为有效值，参数估计偏差的方差达到Cramer-Rao 不等式的下界，即

（4.1.10）

式中M 为Fisher 矩阵，且

（4.1.11）

4）渐近正态性

如果式（4.1.2）中的ξ是均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值?θ

服从正态分布，即

..1

1y uy W P T yu u R R R R R N ??

ΦΦ???→=????

121?T Var E N

N σθ-??????=ΦΦ??

?????

?

?2

1?lim lim 0,..1N N Var R w p N

σθ

-→∞→∞==?θ(){}

121

?T Var E M θσ--=ΦΦ=()

()

?

?ln /ln /??T

p y p y M E θθ

θθ

?????????????

?=?

??

????????

?

?

???

1

2?T -

1.2递推最小二乘法

为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识 设已获得的观测数据长度为N ，则

估计方差矩阵为 式中 于是

如果再获得1组新的观测值，则又增加1个方程

得新的参数估计值 式中

应用矩阵求逆引理，可得1N P +和N P 的递推关系式

矩阵求逆引理 设A 为n ×n 矩阵，B 和C 为n ×m 矩阵，并且A ，A ＋BCT 和I+CTA-1B 都是非奇异矩阵，则有矩阵恒等式

（4.2.2）

得到递推关系式

由于1T

N N N P ψψ+是标量，因而上式可以写成

（4.2.3）

最后，得最小二乘法辨识公式

（4.2.4）

有2种给出初值的办法

（1）设N0（N0>n ）为N 的初始值，可算出

（4.2.5）

N N N

Y θξ=Φ+()1

?N

T T N N N N

Y θ-=ΦΦΦ()1

22

T N N

N Var P θσσ-=ΦΦ=()

1

T N N N P -=ΦΦ?N

T N N N P Y

θ=Φ111

T N N N y ψθξ+++=+()1111

?T N N N N N N P Y y θψ++++=Φ+()1

1111111T

N N T

T T N N N N N N P P ψψψψ---+++++??ΦΦ?????

?==+????

???

???????

()()1

1

1111

T T T A BC A A B I C A B C A ------+=-+()1

1111T T

N N N N N N N N N

P P P I P P ψψψψ-++++=-+()1

11111T T N N N N N N N N N

P P P P P ψψψψ-++++=-+()1

00

00000

?,T

T N N N N N N N P P Y θ-=ΦΦ=Φ()

1111???T N N N N N N

K y θθψθ++++=+-()

1

1111T N N N N N N K P P ψψψ-+++=+()1

111111T T N N N N N N N N N

P P P P P ψψψψ-+++++=-

+（4.2.1）

（2）假定200

?0,,P c I θ==C 是充分大的常数，I 为(2n+1) ×(2n+1)单位矩阵，则经过若干次递推之后能得到较好的参数估计。[1][2][4]

2两种算法的实现方案

2.1最小二乘法一次完成算法实现

如果把式（1.1.6）中的?θ

取好足够的输入—输出数据以后一次计算出来，那么这种算法就式最小二乘法的一次完成法。

2.1.1最小二乘一次完成算法程序框图

2.1.2一次完成法程序 具体程序参见附录4

2.1.3一次完成算法程序运行结果 ans =

1.5000 0.7000 1.0000 0.5000 a1 = 1.5000 a2 = 0.7000 b1 =

1.0000

2.1.4辨识数据比较

2.1.5程序运行曲线

-10

1输入u (k )

0102030405060

-505

输出z (k

)

0102030405060

-5

05离散输出z (k )

2.2递推最小二乘法的实现

2.2.1递推算法实现步骤

1）输入M 序列的产生程序，可以参见3.2当中M 序列产生方法（具体程序见附录。） 2）初始化。

一种简单的方法是选择?0n

θ=和2n P C I =，其中C 为尽量大的数，然后以它们为起始值进行计算(参照式2.2.3)。可以证明，当C →∞时，按照递推公式算得的将与上面用批处理算法求得的结果相同，当N 很大时，C 的选择便无关紧要了。 3）递推算法的停机标准：

()()

()

1??1max

?i i i

i

k k k θθεθ-?--<，ε是适当小的数。

4）最小二乘估计的递推算法系统参数的步骤如下：

（1）产生系统输入信号M 序列，输入系统阶次，计算输入和输出数据u （i ），y （i ），i=1,2,…,n;

（2）置N=n ，10?0,10N N

P I θ==（I …单位矩阵）； （3）形成行向量[]1()(1)()

(1N y n y N n u N u N n ψ+=--+-+-

（4）计算()()()1

11111;T

T

N N N K P N N P N ψ

ψψ-+??=++++??

（5）输入新测量数据y ( N+1 )和u （N+1）；

（6）计算()11???(1)1;N N N N

K y N N θθψθ++??=++-+?

?

（7）计算[]11

1;N N N N P I K P ψ+++=- （8）输出1

?N θ+和1N P +； （9）N ←N+1,转（3）； 2.2.2程序编制思路：

（1）产生M 序列，通过计算，依据算出输出值，设定递推初始值，采用4.2.4方法2设定辨识参数初值，0

?θ为一个很小的量，P0为一个很大值的4×4矩阵，设定误差量，作为停机标准。 （2）依据设定，计算[]

1

(2)(1)(1)(2)T

y y u u ψ=--，1011T P ψ

ψ??+可以算出为一标量，依

据式4.2.5算出K1,K1为4×1矩阵。

（3）依据式4.2.5计算出1?θ，和1

P ，加入估计参数矩阵； （4）计算误差大小，求出相对变化率，加入误差矩阵第一列。

（5）再以1?θ和1

P 为已知参数，重复（2）－（3）计算出参数2?θ，并加入估计参数矩阵。 （6）如此循环，计算出参数估计?n

θ。 （7）判断误差变化率满足要求否，如果满足，则退出循坏，不再进行计算。 （8）分离估计参数，绘出图形。 2.2.3递推最小二乘法程序框图

如下所示：

具体程序参见附录

2.2.4程序运行曲线

-0.5

-0.4-0.3-0.2-0.100.10.20.30.4

0.5

图6M 序列信号的波形

10

20

30

40

50

60

70

80

90

100

-0.4

-0.200.20.40.60.811.21.41.6Parameter Identification with Recursive Least Squares Method

0102030405060708090100

-400

-200020040060080010001200

1400Identification Precision

2.2.5测试结果 见附录6

2.2.6地退数据表

差的分布趋势。（误差=

-新估计值旧估计值旧估计值

）[5]

6结论

最小二乘法是应用广泛的辨识方法之一。可用于动态系统，也可用于静态系统；可用于线性系统，也可用于非线性系统。通过本课程设计，了解和掌握了基于最小二乘法原理的一次性完成算法和递推算法的实现方

法，完成了Matlab下的仿真计算，能够精确的估计辨识参数。最小二乘一次性完成算法是离线算法，需要采集大量数据，一次性完成计算，因此，数据计算量大，当数据量很大时，数据输入不方便，但在本课程设计过程当中，考虑到了此问题，运用相应的办法，解决了矩阵输入的问题。递推算法适合于在线算法，利用原有参数估计进行下一步估计，可以做到运算量小，实时进行估计，根据仿真结果图示，可以看出计算一定量的数据后，参数估计趋于稳定，只要满足误差要求，不必计算完所有数据。两种算法不足之处，在考虑噪声干扰时，看成了白噪声，具有不相关性。如果噪声引起的方程误差是相关的，而不是白噪声，可以通过下面两种方法进行估计①先估计未受噪音污染的系统输出，然后再估计模型的参数，就是辅助变量法；②使用一种滤波器，将相关的方程的误差转换为白噪声再进行估计，就是增广矩阵法。

7参考文献

[1]李言俊张科,系统辨识理论及应用,国防工业出版社,2006.7 .

[2]刘宏才,系统辨识与参数估计,北京,冶金工业出版社,1999.

[3]黄道平 ,MATLAB与控制系统的数字仿真及CAD化学工业出版社, 2004.10

[4]侯媛彬汪梅王立奇系统辨识及其MATLAB仿真科学出版社 2004.2

[5]蔡季冰，系统辨识，北京，北京理工大学出版社，1991.

3附录

附录1 M序列产生程序

function [Out]=Mserise(length,plength,a)

%plength=4;length=15; %M序列周期=2^plength-1 ,%置M序列总长度和级数，a是幅度

for n=1:plength %移位寄存器输入X(i)初T态（1111）

X(n)=1;

end

for i=1:length

for j=plength:-1:2

X(j)=X(j-1);

end

X(1)=xor(X(plength-1),X(plength)); %异或运算

if X(plength)==0

Out(i)=-1;

else

Out(i)=X(plength);

end

end

%save ('mdata','Out');

%绘图

i1=i

k=1:1:length;

plot(k,Out,k,Out,'rx')

xlabel('k')

ylabel('M序列')

title('移位寄存器产生的M序列')

end

附录4程序运行结果如下图所示

图7移位寄存器产生的M序列

附录2 最小二乘一次完成法

%Oncefinish

%model is y(k)+1.5y(k-1)+0.7(k-2)=u(k-1)+0.5(k-2)+v(k);

L=50;A=1;

n=2; %output number;

u=mseries(A,L);

%u=Mserise(L,M,A); %系统辨识的输入信号为一个周期的M序列,信号长度L；幅度A；z=zeros(1,16); %定义输出观测值的长度1*16

dd=zeros(1,4);

gg=zeros(14,1);

for k=n+1:L+1

z(k)=-1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值

end

subplot(3,1,1) %画三行一列图形窗口中的第一个图形

stem(u) %画出输入信号u的经线图形

ylabel('输入u(k)')

subplot(3,1,2)

i=1:1:L+1;

plot(i,z)

ylabel('输出z(k)')

subplot(3,1,3) %画三行一列图形窗口中的第三个图形

stem(z),grid on %画出输出观测值z的经线图形，并显示坐标网格

ylabel('离散输出z(k)')

%u,z; %显示输入信号和输出观测信号

%L=L-1%数据长度

for i=n+1:L+1

h=[-z(i-1) -z(i-2) u(i-1) u(i-2)];

dd(i-2,:)=h; %dd=(L-1)*2n

end

%HL=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14)] %给样本矩阵HL赋值

for j=n+1:L+1

zz=[z(j)];

gg(j-2,:)=zz;

end

%c=inv(dd'*dd)*dd'*gg;

%ZL=[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)]% 给样本矩阵zL赋值

%calculating parameters%计算参数

c1=dd'*dd; c2=inv(c1); c3=dd'*gg; c=c2*c3;c'

%c1=dd'*dd; c2=inv(c1); c3=dd'*ZL; c=c2*c3 %计算并显示

%DISPLAY PARAMETERS

a1=c(1), a2=c(2), b1=c(3), b2=c(4); %从中分离出并显示a1 、a2、b1、b2

%End

附录3递推最小二乘算法程序

clear;clc;%清理工作间变量

L=100;% M序列的周期

y1=1;y2=1;y3=1;y4=0;%四个移位积存器的输出初始值

for i=1:L;%开始循环，长度为L

x1=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或”

x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出

x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出

x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出

y(i)=y4;%取出第四个移位积存器幅值为"0"和"1"的输出信号，

if y(i)>0.5,u(i)=-0.5;%如果M序列的值为"1"时,辨识的输入信号取“-0.03”

else u(i)=0.5;%当M序列的值为"0"时,辨识的输入信号取“0.03”

end%小循环结束

y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备

end%大循环结束，产生输入信号u

figure(1);%第1个图形

stem(u),grid on%以径的形式显示出输入信号并给图形加上网格

z(2)=0;z(1)=0;%取z的前两个初始值为零

for k=3:L;%循环变量从3到25

z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%给出理想的辨识输出采样信号

end

%RLS递推最小二乘辨识

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,99)];%被辨识参数矩阵的初始值及大小

e=zeros(4,100);%相对误差的初始值及大小

for k=3:L; %开始求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<=E break;%若参数收敛满足要求，终止计算

end%小循环结束

end%大循环结束

c%显示被辨识参数

e%显示辨识结果的收敛情况

%分离参数

a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); figure(2);%第2个图形

i=1:L;%横坐标从1到15

plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出a1，a2，b1，b2的各次辨识结果

title('Parameter Identification with Recursive Least Squares Method')%图形标题

figure(3); %第3个图形

i=1:L; %横坐标从1到15

plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('Identification Precision') %图形标题