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

1、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)则式(1.1.1)可写为 （4.1.2）式中：y为N维输出向量；为N为维噪声向量；为（2n+1）维参数向量；为N×（2n+1）测量矩阵。因此，式（4.1.1）是一个含有（2n+1）个未知参数，由N个方程组成的联立方程组。在给定输出向量y和测量矩阵的条件下求参数的估计，这就是系统辨识问题。设 表示 的估计值，表示y的最优估计

2、，则有 （4.1.3）式中：设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)的要求。（4.1.7）当N足够大时有（4.1.8

3、）如果矩阵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是随机的，所以是随机的，但注意不是随机值。如果，则称是无偏估计2）一致性估计误差的方差矩阵为（4.1.9）上式表明当N时，以概率1趋近于。当(k)为不相关随机序列时，最小二乘估计具有无偏性和一致性3)有

4、效性如果式（1.1.2）中的的均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值 为有效值，参数估计偏差的方差达到Cramer-Rao不等式的下界，即（4.1.10）式中M为Fisher矩阵，且（4.1.11）4）渐近正态性如果式（4.1.2）中的是均值为0且服从正态分布的白噪声向量，则最小二乘参数估计值服从正态分布，即（4.1.2）1.2递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识（4.2.1）设已获得的观测数据长度为N，则估计方差矩阵为式中于是如果再获得1组新的观测值，则又增加1个方程得新的参数估计值式中应用矩阵求逆引理，可得和的递推关系式矩阵求逆引

5、理 设A为n×n矩阵，B和C为n×m矩阵，并且A，ABCT和I+CTA-1B都是非奇异矩阵，则有矩阵恒等式（4.2.2）得到递推关系式由于是标量，因而上式可以写成（4.2.3）最后，得最小二乘法辨识公式（4.2.4）有2种给出初值的办法（1）设N0（N0>n）为N的初始值，可算出（4.2.5）（2）假定C是充分大的常数，I为(2n+1) ×(2n+1)单位矩阵，则经过若干次递推之后能得到较好的参数估计。1242两种算法的实现方案2.1最小二乘法一次完成算法实现如果把式（1.1.6）中的取好足够的输入输出数据以后一次计算出来，那么这种算法就式最小二乘法的一次完

6、成法。2.1.1最小二乘一次完成算法程序框图赋输入信号初值u定义输出观测值的长度并计算系统的输出值画出输入和输出观测值的图形给样本矩阵HL和zL赋值计算参数从中分离出并显示出被辨识参数a1, a2, b1, b2停机图3.2 最小二乘一次完成算法程序框图2.1.2一次完成法程序具体程序参见附录42.1.3一次完成算法程序运行结果ans = 1.5000 0.7000 1.0000 0.5000a1 = 1.5000a2 = 0.7000b1 =1.00002.1.4辨识数据比较1.50000.70001.0000.5000参数真值1.5000.7001.000.50误差00002.1.5程序运

7、行曲线2.2递推最小二乘法的实现2.2.1递推算法实现步骤1）输入M序列的产生程序，可以参见3.2当中M序列产生方法（具体程序见附录。）2）初始化。 一种简单的方法是选择和，其中C为尽量大的数，然后以它们为起始值进行计算(参照式2.2.3)。可以证明，当C时，按照递推公式算得的将与上面用批处理算法求得的结果相同，当N很大时，C的选择便无关紧要了。3）递推算法的停机标准：，是适当小的数。4）最小二乘估计的递推算法系统参数的步骤如下：（1）产生系统输入信号M序列，输入系统阶次，计算输入和输出数据u（i），y（i），i=1,2,n;（2）置N=n，（I单位矩阵）；（3）形成行向量（4）计算（5）输入

8、新测量数据y ( N+1 )和u（N+1）；（6）计算（7）计算（8）输出和；（9）NN+1,转（3）；2.2.2程序编制思路：（1）产生M序列，通过计算，依据算出输出值，设定递推初始值，采用4.2.4方法2设定辨识参数初值，为一个很小的量，P0为一个很大值的4×4矩阵，设定误差量，作为停机标准。（2）依据设定，计算，可以算出为一标量，依据式4.2.5算出K1,K1为4×1矩阵。（3）依据式4.2.5计算出，和，加入估计参数矩阵；（4）计算误差大小，求出相对变化率，加入误差矩阵第一列。（5）再以和为已知参数，重复（2）（3）计算出参数，并加入估计参数矩阵。（6）如此循环，计

9、算出参数估计。（7）判断误差变化率满足要求否，如果满足，则退出循坏，不再进行计算。（8）分离估计参数，绘出图形。2.2.3递推最小二乘法程序框图如下所示：具体程序参见附录Y工作间清零产生输出采样信号给被辨识参数和P赋初值按照式(4.2.4)的第三式计算P(k)计算被辨识参数的相对变化量参数收敛满足要求？停机按照式(4.2.4)的第二式计算K(k)按照式(4.2.4)的第一式计算(k)第四个移位寄存器的输出取反，并将幅值变为0.03得到辨识系统的输入信号样本值给M序列的长度L和移位寄存器的输入赋初始值画出被辨识参数的各次递推估计值图形分离参数画出被辨识参数的相对误差的图形画出辨识的输入信号径线图

10、形图5 最小二乘递推算法辨识的Malab程序流程图2.2.4程序运行曲线图6M序列信号的波形2.2.5测试结果见附录62.2.6地退数据表递推次数误差误差误差误差11.500.0011.4990.70.0010.6991.00.0010.9990.50.0010.499921.501.50.700.71.001.00.500.531.50.0011.4990.70.0010.6991.00.750.250.50.750.2541.5-1.500.70.0010.6991.00.750.250.50.750.2551.5-1.500.70.701.00.750.250.50.750.2561.5

11、-1.500.70.701.01.000.50.5071.51.500.70.701.01.000.50.50从以上数据可以看出，随着递推次数的增加，参数估计的误差逐渐减少，最终趋于稳定。下面是辨识误差的分布趋势。（误差=）56结论最小二乘法是应用广泛的辨识方法之一。可用于动态系统，也可用于静态系统；可用于线性系统，也可用于非线性系统。通过本课程设计，了解和掌握了基于最小二乘法原理的一次性完成算法和递推算法的实现方法，完成了Matlab下的仿真计算，能够精确的估计辨识参数。最小二乘一次性完成算法是离线算法，需要采集大量数据，一次性完成计算，因此，数据计算量大，当数据量很大时，数据输入不方便，但

12、在本课程设计过程当中，考虑到了此问题，运用相应的办法，解决了矩阵输入的问题。递推算法适合于在线算法，利用原有参数估计进行下一步估计，可以做到运算量小，实时进行估计，根据仿真结果图示，可以看出计算一定量的数据后，参数估计趋于稳定，只要满足误差要求，不必计算完所有数据。两种算法不足之处，在考虑噪声干扰时，看成了白噪声，具有不相关性。如果噪声引起的方程误差是相关的，而不是白噪声，可以通过下面两种方法进行估计先估计未受噪音污染的系统输出，然后再估计模型的参数，就是辅助变量法；使用一种滤波器，将相关的方程的误差转换为白噪声再进行估计，就是增广矩阵法。7参考文献1李言俊 张科,系统辨识理论及应用,国防工业

13、出版社,2006.7 .2刘宏才,系统辨识与参数估计,北京,冶金工业出版社,1999.3黄道平 ,MATLAB与控制系统的数字仿真及CAD化学工业出版社, 2004.104侯媛彬 汪梅 王立奇 系统辨识及其MATLAB仿真科学出版社 2004.25蔡季冰，系统辨识，北京，北京理工大学出版社，1991.3附录附录1 M序列产生程序function Out=Mserise(length,plength,a)%plength=4;length=15; %M序列周期=2plength-1 ,%置M序列总长度和级数 ，a是幅度 for n=1:plength %移位寄存器输入X(i)初T态（1111）

14、X(n)=1;endfor 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=ik=1:1:length;plot(k,Out,k,Out,'rx')xlabel('k')ylabel('M序列')title('移位

15、寄存器产生的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*16dd=zeros(1,4);gg=zeros(14,1);for k=n+1:L+1 z(k)=-1.5

16、*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值endsubplot(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; %显示输入信号和

17、输出观测信号%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)*2nend%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(1

18、1) -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赋值%c

19、alculating 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 PARAMETERSa1=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;%四个移位积存器的输出初

20、始值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;%当

21、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.00

22、1'%直接给出被辨识参数的初始值,即一个充分小的实向量p0=106*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=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;%求参数当前值与上一次的值的差值