系统辨识—最小二乘法

1、最小二乘法参数辨识201403027摘 要：系统辨识在工程中的应用非常广泛,系统辨识的方法有很多种,最小二乘法是一种应用极其广泛的系统辨识方法阐述了动态系统模型的建立及其最 小二乘法在系统辨识中的应用，并通过实例分析说明了最小二乘法应用于系统辨 识中的重要意义.关键词:最小二乘法；系统辨识；动态系统Abstract : System identification in engineering is widely used, system iden tificati on methods there are many ways ， least squares method is a very w

2、ide range of applicati on of system ide ntificati on method and the least squares method elaborated establish a dyn amic system models in System Ide ntificati on applicati ons and examples an alyzed by the least squares method is applied to illustrate the importanee of system iden tificati on.Keywor

3、ds: Least Squares; system ide ntificati on; dyn amic system引言随着科学技术的不断发展，人们认识自然、利用自然的能力越来越强，对于未 知对象的探索也越来越深入.我们所研究的对象，可以依据对其了解的程度分为 三种类型：白箱、灰箱和黑箱如果我们对于研究对象的内部结构、 内部机制了解 很深入的话，这样的研究对象通常称之为“白箱”；而有的研究对象，我们对于其 内部结构、机制只了解一部分，对于其内部运行规律并不十分清楚，这样的研究对 象通常称之为“灰箱”；如果我们对于研究对象的内部结构、内部机制及运行规 律均一无所知的话，则把这样的研究对象称之为

4、“黑箱” 研究灰箱和黑箱时，将 研究的对象看作是一个系统，通过建立该系统的模型，对模型参数进行辨识来确 定该系统的运行规律对于动态系统辨识的方法有很多，但其中应用最广泛，辨识 效果良好的就是最小二乘辨识方法，研究最小二乘法在系统辨识中的应用具有现 实的、广泛的意义1. 1系统辨识简介系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。 现代控制理论中的一个分支。通过辨识建立数学模型的目的是估计表征系统行为 的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入 和输出预测系统输出的未来演变，以及设计控制器。对系统进行分析的主要问题 是根据输入时间函数和系统的特性来

5、确定输出信号。对系统进行控制的主要问题 是根据系统的特性设计控制输入，使输出满足预先规定的要求。而系统辨识所研 究的问题恰好是这些问题的逆问题。通常，预先给定一个模型类卩=M（即给定一类已知结构的模型），一类输入信号u和等价准则J=L（y, yM）（ 般情况下， J是误差函数，是过程输出y和模型输出yM的一个泛函）；然后选择使误差函数 J达到最小的模型，作为辨识所要求的结果。系统辨识包括两个方面：结构辨识 和参数估计。在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计 这两个方面并不是截然分开的，而是可以交织在一起进行的。1.2系统辨识的目的在提出和解决一个辨识问题时，明确最终使用模型

6、的目的是至关重要的。 它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。通过辨 识建立数学模型通常有四个目的。 估计具有特定物理意义的参数有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。这就需要 通过能观测到的输入输出数据，用辨识的方法去估计那些参数。 仿真仿真的核心是要建立一个能模仿真实系统行为的模型。用于系统 分析的仿真模型要求能真实反映系统的特性。用于系统设计的仿真，则强调设计 参数能正确地符合它本身的物理意义。 预测 这是辨识的一个重要应用方面，其目的是用迄今为止系统的可测 量的输入和输出去预测系统输出的未来的演变。例如最常

7、见的气象预报，洪水预 报，其他如太阳黑子预报，市场价格的预测，河流污染物含量的预测等。预测模型辨识的等价准则主要是使预测误差平方和最小。只要预测误差小就是好的预测 模型，对模型的结构及参数则很少再有其他要求。这时辨识的准则和模型应用的 目的是一致的，因此可以得到较好的预测模型。 控制 为了设计控制系统就需要知道描述系统动态特性的数学模型，建 立这些模型的目的在于设计控制器。 建立什么样的模型合适，取决于设计的方法 和准备采用的控制策略。2最小二乘方法2.1.1系统辨识最小二乘方法简介对工程实践中测得的数据进行理论分析，用恰当的函数去模拟数据原型是一 类十分重要的问题，最常用的逼近原则是让实测数

8、据和估计数据之间的距离平方 和最小，这即是最小二乘法。最小二乘法是一种经典的数据处理方法。在系统辨 识领域中，最小二乘法是一种得到广泛应用的估计方法，可用于动态系统，静态系统，线性系统，非线性系统。可用于离线估计，也可用于在线估计。这种 辨识方法主要用于在线辨识。在随机的环境下，利用最小二乘法时，并不要求观 测数据提供其概率统计方面的信息，而其估计结果，却有相当好的统计特性。MATLAB是一套高性能数字计算和可视化软件，它集成概念设计，算法开发，建模仿真，实时实现于一体，构成了一个使用方便、界面友好的用户环境，其强大的扩展功能为各领域的应用提供了基础。对于比较复杂的生产过程，由于过程的输入输出

9、信号一般总是可以测量 的，而且过程的动态特性必然表现在这些输入输出数据中，那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。这种建模方法就称为系统辨识。把辨识建模称作“黑箱建模”。2.1.2 最小二乘法系统辨识结构：本文把待辨识的过程看作“黑箱”。只考虑过程的输入输出特性，而不强调 过程的内部机理。图中，输入u(k)和输出z(k)是可以观测的；G是系统模型，用来描述系统的输 入输出特性；N是噪声模型，v(k)是白噪声，e(k)是有色噪声，根据表示定理： 可以表示为e(k) =N v(k)G(z)B ( z "I)A(z-)N ( z -1)D(z)C(z J)A(z'

10、;) = 1 a# azZ l| a“az 入B(z') = b才bzZ，III tnbnb；C(z-1) = 1 + ch1 + C2Z + cz-n D(z-1H dQ1 + dzz2 + dAbnb2.1.3准则函数设一个随机序列'z(k)，kE (1,2,丄)的均值是参数日的线性函数： Etz(k)'hT(k户，其中h(k)是可测的数据向量，那么利用随机序列的一个实现， 使准则函数：L2J(R 二' z(k) -hT(kp 心(式2-2)达到极小的参数估计值二称作二的最小二乘估计。最小二乘格式：z(k)二ht(k户e(k)， 日为模型参数向量，ek)为零

11、均值随机2.2广义最小二乘法2.2.1广义最小二乘数学模型A(z4)z(kB(zJ)u(kr-1v(k)C(z )式中，u(k)和z(k)表示系统的输入输出；v(k)是均值为零的不相关的随机 序列；且'a(z)=1 +aZ +a2Z，* +anaz"aBCz'Hdz'+bzz'i+bnbzH'C(z') =1 +C1Z +C2Z " +Cnczf222广义最小二乘递推算法如下攻k)=畝 _1) + K f (k)z,k) hf(k)坯k 1)K ©) = Pf(k-1)hf(k)h?(k)Pf(k-1)hf(k)+1

12、卩 f (k)=l - K f(k)h«k)Pf(k-1):? (k) "(k -1) K e(k)g(k) - he(k)玖 k-1)K e(k) = Pe(k-1) h e(k) h e(k) P°(k - 1) h e(k)1丄Pe(k)珂I -K e(k)h e(k) Pe(k_1)式中hf(k)=乙(k-1),，-乙(k rOu(k 1),，q(k n)F he(k)=Gk-1),-§(k-nc)re(k)=z(k)- h Yk)k)2.2.3广义最小二乘递推算法的计算步骤：1.给定初始条件?(0)=；(充分小的实向量)Pf(0) = a2l

13、(a为充分大的数)?e(0)= 0fe(O)=l2利用式乙(kC(zA)z(k) (k)二 C(z 丄)z(k)计算 Zf(k)和 u«k)；3利用式沪q,，比力,bjhf(k)二-乙(k-1),厂z(k rQ,Uf(k1),，u© Ob)T构造hf(k)；%k)二珂k -1) K f (k)zf(k) - hf(k)珂k - 1)4利用式 K ,k)二 Pf(k -1)hf(k)h；(k)Pf (k -1)h,k) 1戸递推计算?(k)；P f(k)-I -K f(k)h ;(k) P f (k - 1)5 利用(k)二 z(k) - h (k)讹k)和h (k)十 z(

14、k-1), ，-z(k- ra),u(k-1), ,u(k-n)计算锹);6根据 he(k) =e(k 1),，e(k - nc)广来构造 he(k)；綬(k) = (?(k -1) + K e(k)&k) - hHk)8?(k - 1)7 利用 K e(k)二 Pe(k-1)h e(k) 'h e(k) Pg(1)h g(k)什Pe(k)=l - K e(k) hj(k) R(k-1)返回第2步进行迭代计算，直至获得满意的辨识结果3工程实例3.1典型系统建模以某微循环流体系统模型的参数辨识为例我们已经得到该系统模型的差分 方程形式，取特定点的压力波作为模型的输入，以另一点的压

15、力波作为模型的输 出由于我们采集的数据是实时的，因此采用在线辨识方法.由于建立的微循环流 体系统模型是一个单输入、单输出的模型，为使参数估计的结果很好地跟踪参数 真值的变化，我们采用渐消记忆的最小二乘法对系统模型参数进行辨识，即强调新数据的作用，贬低老数据的作用抽象出的SISO系统的差分方程为：z(k) Qz(k - 1) a2z(k - 2) =bi(k - 1) b2u(k - 2)(k) 式( 1-1)参数取真值为：J二376 0.483 0.57 0.42】，利用MATLAB的M语言辨识 系统中的未知参数3、a2、“、b2。要求：用参数的真值利用差分方程求出z(k) 作为测量值，：(k

16、)是均值为0，方差为0.1、0.5和0.01的不相关随机序列。 使用最小二乘算法辨识。3.2广义最小二乘递推算法的 MATLA仿真(程序源代码见附录)考虑仿真对象z(k)= -1.376z(k-1)-0.483z(k-2)+0.57u(k-1)+0.42u(k-2)+v(k)式中，v(k)是均值为0,方差为0.01、0.1和0.5的不相关随机序列。输入信 号采用4阶M序列，幅度为1。选择如下形式的辨识模型图3.1广义最小二乘法辨识实例结构图其中取 c1=0,c2=0.4结果分析及算法优化由于辨识算法中输入或噪声信号为不相关随机序列，所以每次辨识结果都不 完全相同。但是，在相同输入、相同的噪声、

17、相同的步长条件下，精度大体相同。算法优化方案：（1）使用M序列（具有近似白噪声的性质）为输入信号；（2）增加数据长度去L；（3）减小噪声信号v（k）的方差。4.1广义最小二乘递推算法的的 MATLA仿真结果及分析（1）、输入选用题目给出的30个随机数，即数据长度去L=30，噪声选用均值 零，方差分别为0.5、0.1和0.01的随机序列，辨识结果如表表 4-1。表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差表4-1真值噪声方差为0.5噪声方差为0.1噪声方差为0.01估计值相对误差估计值相对误差估计值相对误差a11.3761.5534-2.01051.47300.1231

18、1.37230.0063a20.4830.6884-2.00540.58920.10610.48150.0231b10.570.7916-2.02160.67900.18210.57990.0510b20.420.6213-2.00130.53610.11610.42180.0234（2）、输入均采用M序列，噪声选择均值为零，方差为0.5、0.1和0.01 的随机序列，辨识步长均为300步，辨识结果如表4-2。表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差.表4-2真值噪声方差为0.5噪声方差为0.1噪声方差为0.01估计值相对误差估计值相对误差估计值相对误差a11.3

19、761.3660-0.02801.37500.00791.37290.0005a20.4830.4749-0.07010.48590.02920.4710-0.0068b10.570.6613-0.12490.58270.05820.57200.0051b20.420.4112-0.08230.43440.09830.4183-0.0049(3) 数据结果分析：输入采用 M序列比采用随机序列得到的辨识效果更 好。噪声均值相等时，方差越大，辨识效果越差，反之，方差越小辨识效果越 好。可以通过增加步长的方法提高辨识精度。下面给出以M序列作为输入，噪声均值为零，方差为 0.01的随机序列，数 据长度

20、取L=30,得到的变化曲线图：下面给出以M序列作为输入，噪声均值为零，方差为 0.01的不相关随机序 列，数据长度去L=300,得到的变化曲线图：File Edi t View Insert Toob Desktop Window HelpD Q S k貝Q的瓷盟h 螯皱变化曲线1. i户L-L11 1 1 05Q100150200290300参考文献1 李言俊，张科，系统辨识理论及应用，国防工业出版社，2006年2 方崇智，萧德云，过程辨识，清华大学出版社，2002年3 贾秋玲，袁冬莉，栾云凤，基于 MATLAB7.x/Simulink/Stateflow 系统仿真、分析及设计，西北工业大学

21、出版社，2006年侯媛彬,汪梅,王立琦,系统辨识及其MATLAB真,科学出版社,2004年附录广义最小二乘递推算法的MATLA仿真程序源代码：clear %清理工作间变量L=300; % M序列的周期 y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值for i=1:L;% 开始循环，长度为 Lx1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出y(

22、i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即 M序列if y(i)>0.5,u(i)=-1; %如果M序列的值为"1",辨识的输入信号取“ -1”else u(i)=1; % 如果M序列的值为"0",辨识的输入信号取“ 1”end %小循环结束y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备end %大循环结束，产生输入信号ufigure(1); % 第一个图形stem(u),grid on %显示出输入信号 M序列径线图并给图形加上网格v=normrn

23、d(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*ze(k-2)+v(k-1);%C(z1)=1,即取 c仁0,c2=0endz(2)=0;z(1)=0; %设z的前两个初始值为零for k=3:301; %循环变量从 3到301z(k)=-1.376*z(k-1)-0.483*z(k-2)+57*u(k-1)+0.42*u(k-2)+ze(k-1); %输出采样信号(测量值)end%RGL旷义最小二乘辨识c0=0.0001 0.00

24、01 0.0001 0.0001' %直接给出被辨识参数的初始值,即一个充分小的实向量pf0=10A6*eye(4,4); %直接给出初始状态P0,即一个充分大的实数单位矩阵ce0=0.001 0.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);

25、uf(k)=u(k)+ce(1,k-2)*u(k-1)+ce(2,k-2)*u(k-2); hf1=-zf(k-1),-zf(k-2),uf(k-1),uf(k-2)'x=hf1'*pf0*hf1+1; x1=inv(x); %开始求 K(k)k1=pf0*hf1*x1;% 求出 K 的值d仁 zf(k)-hf1'*c0; c仁 c0+k1*d1; %求被辨识参数 ce仁c1-c0; %求参数当前值与上一次的值的差值e2=e1./c0; %求参数的相对变化e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列c0=c1; %新获得的参数作为下一次递推的旧参

26、数c(:,k)=c1; %把辨识参数c列向量加入辨识参数矩阵的最后一列pf1=pf0-k1*hf1'*pf0; %求出 p(k)的值pf0=pf1; %给下次用h1=-z(k-1),-z(k-2),u(k-1),u(k-2)'s=s+(z(k)-h1'*1.642 0.715 0.39 0.35')A2;%求准则函数ee(k)=z(k)-h1'*c1;he1=-ee(k-1),-ee(k-2)'x=he1'*pe0*he1+1; x1=i nv(x);k仁 pe0*he1*x1;d仁 ee(k)-he1'*ce0;ce仁 ce0+

27、k1*d1;pe1=pe0-k1*he1'*pe0;ce0=ce1;ce(:,k)=ce1;pe0=pe1;end %大循环结束c%辨识参数变化矩阵%显示被辨识参数及其误差(收敛)情况%分离参数a1=c(1,1:300); a2=c(2,1:300);b 仁c(3,1:300);b2=c(4,1:300);c仁 ce(1,1:300);c2=ce(2,1:300);ea仁 e(1,1:300); ea2=e(2,1:300);eb 仁e(3,1:300); eb2=e(4,1:300); figure(2); %第二个图形画出al，b1，b2,的各次i=1:300; % 横坐标从1到3

28、00Plot(i,a1,'r',i,a2,'k',i,b1,'b',i,b2,'c',i,c1,'b',i,c2,'r') %a2,b1，b2,c1,c2 的各次辨识结果title(' 参数变化曲线')%图形标题figure(3); %第三个图形i=1:300; % 横坐标从1到300 plot(i,ea1,'r',i,ea2,'k:',i,eb1,'b',i,eb2,'k:')%画出 al，a2,辨识结果的收敛情况ti

29、tle('误差曲线')%图形标题考虑仿真对象 z(k)+1.5z(k-1)+0.7z(k-2)=u(k-1)+0.5u(k-2)+v(k),其中v(k)是服从正态分布的白噪声N(0,1)。输入信号采用4阶M序列。选择如下形式的辨识模型z(k)+a1z(k-1)+a2z(k-2)=b1u(k-1)+b2u(k-2)+v(k)，试利用递推最小二乘法辨识参数al、a2、bl、b2。Np=15;r=4;X1=1;X2=1;X3=1;X4=1;men gth = r*Np;a=1;for i=1:1: m_len gthY4=X4;Y3=X3;Y2=X2;Y 仁X1;X4=Y3;X3=Y

30、2;X2=Y1;X仁 xor(Y3,Y4);if Y4=0M(i)=-a;elseM(i)=a;endendfigure;i=1:1:m_le ngth;plot(i,M);%白噪声no ise = zeros(1,m_le ngth);for i=1:1: m_len gthtemp = no ise + 0.5*ra nds(1, m_len gth);no ise = temp;endno ise = no ise/12;%no ise = temp;% parameter of systemn=2;d=1;a 1= -1;a2=0.5;b1=1;b2=0.5;S_U0=0.2;S_Y0=0.2;% gen erate u,yuO=o nes(1,me ngth)*S_UO;U = M + u0 + noise;figure;i=1:1:m_le ngth;plot(i,U);%formulatio ny(1)=0;y(2)=0;y(3)=0;Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise (2);Y(3)=S_Y0