系统辨识最小二乘法

1、课程设计报告学 院:自动化学院专业名称:自动化学生姓名: 指导教师:时 间:2010年7月课程设计任务书、设计内容SISO系统的差分方程为：z(k) + a z(k -1) + a z(k - 2) = b u(k -1) + b u(k - 2) +u (k)1212参数取真值为：0 T=h.642 0.715 0.39 0.35,利用MATLAB的M语言辨识系 统中的未知参数a1、a2、b1、b2。二、主要技术要求用参数的真值及差分方程求出z(k)作为测量值，u (k)是均值为0，方差为 0.1、0.5和0.01的不相关随机序列。选取一种最小二乘算法利用MATLAB的M 语言辨识参数。三、

2、进度要求2周(6月28日-7月11日)完成设计任务，撰写设计报告3000字以上， 应包含设计过程、计算结果、图表等内容。具体进度安排：6月28日，选好题目，查阅系统辨识相关最小二乘法原理的资料。6月29日，掌握最小二乘原理，用MATLAB编程实现最小二乘一次完成 算法。6月30日，掌握以最小二乘算法为基础的广义最小二乘递推算法。7月1日，用MATLAB编程实现广义最小二乘递推算法。7月2日，针对题目要求进行参数辨识，并记录观察相关数据。7月3日-7月5日，对参数辨识结果进行分析，找出存在的问题，提出 改进方案，验证改进优化结果。7月6日-7月7日，撰写课程设计报告。7月8日，对课程设计报告进行

3、校对。7月9日，打印出报告上交。学 生 王景指导教师邢小军1.设计内容设SISO系统的差分方程为：z(k) + a z(k -1) + a z(k - 2) = bu(k -1) + b u(k - 2) +u (k)式(1-1)1212参数取真值为：T = h.642 0.715 0.39 O.35,利用MATLAB的M语言辨识系统中的未知参数a、a、b、b。12122.设计过程问题重述。设SISO系统的差分方程为：z (k) + a z(k -1) + a z (k - 2) = b u (k -1) + b u (k - 2) +u (k)式(2-1)1212参数取真值为：T = 1.6

4、42 0.715 0.39。.351利用MATLAB的M语言辨识系统中的未知参数a、a、b、b。 1212要求：用参数的真值利用差分方程求出z(k)作为测量值，U (k)是均值为0， 方差为0.1、0.5和0.01的不相关随机序列。选取一种最小二乘算法辨识。最小二乘参数辨识2.2.1、最小二乘法的概念与应用对工程实践中测得的数据进行理论分析，用恰当的函数去模拟数据原型是 一类十分重要的问题，最常用的逼近原则是让实测数据和估计数据之间的距离 平方和最小，这即是最小二乘法。最小二乘法是一种经典的数据处理方法。在 系统辨识领域中，最小二乘法是一种得到广泛应用的估计方法，可用于动态系 统，静态系统，线

5、性系统，非线性系统。可用于离线估计，也可用于在线估计。 这种辨识方法主要用于在线辨识。在随机的环境下，利用最小二乘法时，并不 要求观测数据提供其概率统计方面的信息，而其估计结果，却有相当好的统计 特性。MATLAB是一套高性能数字计算和可视化软件，它集成概念设计，算法开 发，建模仿真，实时实现于一体，构成了一个使用方便、界面友好的用户环境， 其强大的扩展功能为各领域的应用提供了基础。对于比较复杂的生产过程，由 于过程的输入输出信号一般总是可以测量的，而且过程的动态特性必然表现在 这些输入输出数据中，那么就可以利用输入输出数据所提供的信息来建立过程 的数学模型。这种建模方法就称为系统辨识。把辨识

6、建模称作“黑箱建模”。2.2.2、最小二乘法系统辨识结构：本文把待辨识的过程看作“黑箱”。只考虑过程的输入输出特性，而不强调 过程的内部机理。e(k)图1 SISO系统辨识“黑箱”结构图v(k)e(k)图1 SISO系统辨识“黑箱”结构图v(k)u(k)+y(k)z(k)G( z-1)N (z-1)图中，输入u(k)和输出z(k)是可以观测的；G (二-)是系统模型，用来描述系 统的输入输出特性；N (二-)是噪声模型，v(k)是白噪声，e(k)是有色噪声，根据表示定理：可以表示为e(k) =N (厂-)v(k)G (z -1) = B(z-Q N (z -1) = D (z-1)+ a z

7、- nana+ b z - nb. . n+ a z - nana+ b z - nb. . nbA( z-1) = 1 + a z-1 + a z -2 + 12B( z-1) = bz-1 + b z -2 +12C ( z -1) = 1 + c z -1 + c z -2 + c z-气12气D(z-1) = d z-1 + d z-2 + d z-nbI12n,b2.2.3、准则函数设一个随机序列(k)，k e (1，2，.,异的均值是参数0的线性函数：(kR=阮(k)0，其中h(k)是可测的数据向量，那么利用随机序列的一个实现， 使准则函数：J (0) = * z(k) - hT

8、(k)0 2k=1(式 2-2)达到极小的参数估计值#称作0的最小二乘估计。最小二乘格式：z(k)= ht(k)0+ e(k)，0为模型参数向量，上)为零均值 随机噪声。最小二乘一次完成算法2.3.1、最小二乘问题的解考虑系统模型：（式 3-1）（式 3-2）（式 3-3）z (k) = ht (k )0 （式 3-1）（式 3-2）（式 3-3）准则函数可写成：J (0)=(z - HQ a l(zl - H?)极小化准则函数得到：0“=(Htah )】Hta zWLS L L L L L通过极小化式3-2计算0 WLS的方法称作加权最小二乘法，0 WLS为加权最小 二乘估计，若取A L =

9、 1，则退化为一般最小二乘估计值，对应方法叫最小二乘 法：L；(HlH )-iHlz(式 3-4)当获得一批数据后，利用式3-3或3-4可一次求得相应参数估计值，对应的0称为最小二乘估计值，这样处理问题的方法就称作一次完成算法。LS输入信号：随机序列(如白噪声)；伪随机序列(如M序列或逆M序列)；离散序列2.3.2、最小二乘一次完成算法的MATLAB仿真(程序源代码见附录)考虑仿真对象z(k)= -1.642z(k-1)-0.715z(k-2)+0.39u(k-1)+0.35u(k-2)+v(k) 式中，v(k)是均值为0，方差为0.01. 0.1和0.5的不相关随机序列。输入信 号采用4阶M

10、序列，幅度为1。选择如下形式的辨识模型z(k) + a z(k -1) + a z(k - 2) = b u(k -1) + b u(k - 2) +u (k)，1212构造ZL和Hl，数据长度去L=30；加权阵取I，利用式(式3-4)计算参数估计值 人0LS赋输入信号初值u赋输入信号初值u0 =aa 2， z =-z(3) z，H =-z(2) -z(1) u(2) u(1)- -z(3) - z(2) u(3) u(2)LS* b2L一 z (30).L _-z(29) -z(28) u(29) u(28)_程序框图如图3.2所示。Matlab6.0仿真程序如下:定义输出观测值的长度并计算

11、系统的输出值显示输入和输出观测值给样本矩阵hl和zL赋值根据式(3-4 )计算参数0。I LS从0lS中分离出并显示出被辨识参数气a b1, b停机广义最小二乘法2.4.1、广义最小二乘数学模型A( z -i) z(k) = B( z - i)u(k) +1v(k)C( z -1)式中，u(k)和z(k)表示系统的输入输出；v(k)是均值为零的不相关的随机 序列；且 ， 一、 - . . 一.A(z-i) = 1 + %z-i + a2z-2 + a z-na B(z-i) = b z-i + b z-2 + b z-nb12nC ( z -1) = 1 + c z -1 + c z - 2

12、+ c z - ncI12nc2.4.2、广义最小二乘递推算法如下6(k) =0(k-1) + K (k)z (k)-h(k)6(k -1)K (k) = P (k - 1)h(k )k(k )P (k - 1)h (k) +11P (k) = I - K (k)K(k)f (k -1) f TOC o 1-5 h z ff6 (k) =6 (k -1) + K (k)e(k) -hT (k)6 (k -1) K (k) = P (k - 1)h (k)k (k)P (k - 1)h(k) +11P (k) = I - K (k)丘(k)P (k -1) eV eeee式中h (k) = -z

13、 (k-i),(k -n ),u (k-1),u (k -n)hf(k) = -e(k-1),-,-e(kf -n )af f be八c(k) = z(k) - h (k)6 (k)2.4.3、广义最小二乘递推算法的计算步骤:0(0) f (充分小的实向量)P(0) = a2I (a为充分大的数)1.给定初始条件Lf0 (0) = 0eP(0) = Iez (k) = C(z_i)z(k)2利用式/，计算z (k)和u (k)；u (k) = C(z_i)z(k)f ffo = a，,a , b，,b 1构造h (k)；f TOC o 1-5 h z 3 利用式 ,11 构造h (k)；fh

14、(k) = -z (k1),z (k-n ),u (k1),u (k-n)fffa ffb京446(k) =6(k -1) + K (k)z (k) -h(k)6(k -1)4 利用式Jk (k) = P (k-1)h(k)lh(k)P (k - 1)h (k) +11 递推计算0(k)；ffffffP (k) = I - K (k)h(k)P (k -1)L ff ff5 利用 e(k) = z(k) - hT(k)6 (k)和h(k) = z(k -1),，z(k n ), u(k -1),u(k n )t 计算 2(k)； ab根据h (k) = 一e(k - D，一e(k 一 n )T

15、来构造h (k)；6 (k) =6(k -1) + k (k)e(k) -h: (k)6 (k -1)利用K (k) = P (k - 1)h (k)lh(k)P (k - 1)h (k) +1-11 eee_ee P (k) = I - K (k)h (k)P (k -1)eeee返回第2步进行迭代计算，直至获得满意的辨识结果。2.4.4、广义最小二乘递推算法的MATLAB仿真(程序源代码见附录)考虑仿真对象z(k)= -1.642z(k-1)-0.715z(k-2)+0.39u(k-1)+0.35u(k-2)+v(k)式中，v(k)是均值为0，方差为0.01、0.1和0.5的不相关随机序列

16、。输入信 号采用4阶M序列，幅度为1。选择如下形式的辨识模型其中取 c1=0,c2=0.3.结果分析及算法优化由于辨识算法中输入或噪声信号为不相关随机序列，所以每次辨识结果都不 完全相同。但是，在相同输入、相同的噪声、相同的步长条件下，精度大体相同。算法优化方案：（1）使用M序列（具有近似白噪声的性质）为输入信号；（2）增加数据长度去L；（3）减小噪声信号v（k）的方差。3.1、最小二乘一次完成算法仿真结果及分析：采用给定的30组随机数据作为输入，即数据长度去L=30，噪声v（k）选用 均值为零，方差分别为0.5、0.1和0.01辨识结果如表3-1-1：表 3-1-1真值噪声方差为0.5噪声方

17、差为0.1噪声方差为0.01估计值相对误差估计值相对误差估计值相对误差a11.6421.6184-0.01441.64280.00051.64130.0004a20.7150.6851-0.04180.72540.01450.73960.0344b10.390.45350.16280.41490.06380.39410.0105b20.350.37270.06490.37780.07940.37250.0643由辨识结果可知，v(k)的方差越小，辨识效果越好。但由于v(k)噪声是均 值为零，方差分别为0.5、0.1和0.01的不相关随机序列，故一次完成算法每次 辨识结果都不一样，得到的参数与真

18、值间的误差较大，不太理想。如V(k)选用均值为零，方差为0.1的不相关随机序列。增加数据长度L=300, 效果如表3-1-2：表 3-1-2真值噪声方差为0.1，且输A为M序列， 采集数据30个噪声方差为0.1，且输A为M序列，米 集数据300个估计值相对误差估计值相对误差a11.6421.64280.00051.64390.0011a20.7150.72540.01450.7125-0.0034b10.390.41490.06380.39630.0162b20.350.37780.07940.35480.0137通过数据比较(见上表3-1-1和3-1-2)，可知增加数据长度L效果明显， 但一

19、次完成算法计算数据量很大；需要计算矩阵的逆，当数据很多时，会出还 现“数据饱和”现象。3.2、广义最小二乘递推算法的的MATLAB仿真结果及分析、输入选用题目给出的30个随机数，即数据长度去L=30，噪声选用 均值为零，方差分别为0.5、0.1和0.01的随机序列，辨识结果如表3-2-1表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相 对误差表 3-2-1真值噪声方差为0.5噪声方差为0.1噪声方差为0.01估计值相对误差估计值相对误差估计值相对误差a11.6421.4555-0.11351.56300.12311.65230.0063a20.7150.5884-0.17710

20、.63920.10610.73150.0231b10.390.49160.26050.31900.18210.40990.0510b20.350.42130.20140.34610.01140.34180.0234输入M序列，30步，噪声方差0.5时:8 =二：=二k 一 二、一 日一 =12.556;输入M序列，30步，噪声方差0.1时:,9 二二二-广=2.5822;输入M序列，30步，噪声方差0.01时：=&_ 二三一1乙二0.1706;、输入均采用M序列，噪声选择均值为零，方差为0.5、0.1和0.01 的随机序列，辨识步长均为300步，辨识结果如表3-2-2。表中给出了三种情况下辨识

21、参数结果即表中的估计值，估计值与真值的相 对误差.表 3-2-2真值噪声方差为0.5噪声方差为0.1噪声方差为0.01估计值相对误差估计值相对误差估计值相对误差a11.6421.5960-0.02801.65500.00791.64290.0005a20.7150.6649-0.07010.73590.02920.7101-0.0068b10.390.3413-0.12490.41270.05820.39200.0051b20.350.3212-0.08230.38440.09830.3483-0.0049输入采用M序列，噪声方差0.5时：9 =!：=.=-0.5,u(i)=-1; %如果M序

22、列的值为1,辨识的输入信号取“-1” else u(i)=1; %如果M序列的值为0,辨识的输入信号取“1” end %小循环结束 y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备 end %大循环结束，产生输入信号u figure(1); %第一个图形 stem(u),grid on %显示出输入信号M序列径线图并给图形加上网格 v=normrnd(0, sqrt(0.01), 1, 300);%均值为零的，方差为 0.01 或 0.5 或 0.1 不相 关的随机噪声 ze(2)=0;ze(1)=0; for k=3:301;ze(k)=0*ze(k-1)+0*z

23、e(k-2)+v(k-1);%C(z1)=1,即取 c1=0,c2=0 end z(2)=0;z(1)=0; %设2的前两个初始值为零 for k=3:301; %循环变量从3到301 z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.351*u(k-2)+ze(k-1);嘛 出采样信号(测量值) end%RGLS广义最小二乘辨识c0=0.0001 0.0001 0.0001 0.0001; %直接给出被辨识参数的初始值，即一个充分 小的实向量pf0=106*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵 ce0=0.001 0

24、.001;pe0=eye(2,2);c=c0,zeros(4,299); %被辨识参数矩阵的初始值及大小 ce=ce0,zeros(2,299);e=zeros(4,300); %相对误差的初始值及大小 ee=zeros(2,300); s=0;%广义最小二乘递推算法的计算步骤 for k=3:300;zf(k)=z(k)+ce(1,k-2)*z(k-1)+ce(2,k-2)*z(k-2); uf(k)=u(k)+ce(1,k-2)*u(k-1)+ce(2,k-2)*u(k-2); hf1=-zf(kT),-zf(k-2),uf(kT),uf(k-2);x=hf1*pf0*hf1+1; x1=inv(x); %开始求 K(k)k1=pf0*hf1*x1;%求 出 K 的值 d1=zf(k)-hf1*c0; c1=c0+k1*d1; %求被辨识参数 c e1=c1-c0; %求参数当前值与上一次的值的差值 e2=e1./c0; %求参数的相对变化e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列c0=c1; %新获得的参数作为下一次递推的旧参数c(:,k)=c1; %把辨识参数