系统辨识—最小二乘法汇总

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

2、s a very wide range of application of system identification method and the least squares method elaborated establish a dynamic system models in System Identification applications and examples analyzed by the leastsquares method is applied to illustrate the importance of system identification.Keyword

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

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

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

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

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)1、1G ( z 1) B(z J N

10、 ( z 1) D(z JA(z 1)C(z 1)A(z 1) 1 4z1 a2Z2 LanaZ naB(z1) Dz1 b2z2 L 吐112C(z ) 1 cz c?zCnazD(z 1) d逐1 d2z 2dnbzn2.1.3准则函数设一个随机序列z(k),k (1,2, ，L)的均值是参数 的线性函数: Ez(k) hT(k),其中h(k)是可测的数据向量，那么利用随机序列的一个实现， 使准则函数：L2J( )z(k) hT(k)k 1(式 2-2)达到极小的参数估计值称作的最小二乘估计。最小二乘格式：z(k) ht(k) e(k),为模型参数向量，ek为零均值随机2.2广义最小二乘法

11、2.2.1 广义最小二乘数学模型A(z1)z(k)11B(z1)u(k)-v(k)C(z1)式中,u(k)和z(k)表示系统的输入输出；v(k)是均值为零的不相关的随机 序列；且A(z 1)1a1z 1a2z2a%z naB(z 1)z1 b2z 2bnbz nbC(z 1)1Gz 1C2z2Cncz nc2.2.2 广义最小二乘递推算法如下?(k)?(k 1) Kf(k)zf(k) hf(k)?(k 1)1Kf(k) Pf(k 1)hf(k)hf(k)Pf(k 1)hf(k) 1 1Pf (k)I Kf(k)hf(k)Pf(k 1)?e(k)?e(k 1) Ke(k)e?(k) he(k)

12、?e(k 1)1Ke(k)Pe(k 1)he (k) he(k)Pe(k 1)he(k) 1 1Pe(k) IKe(k)he(k)Pe(k 1)式中hf (k)zf(k 1), ,zf(kna),uf (k1),uf (knb)he(k) e?(k 1), , e?(k nc)e?(k) z(k) h (k) ?(k)2.2.3 广义最小二乘递推算法的计算步骤：?(0)(充分小的实向量)1. 给定初始条件2 利用式zf (k)uf(k)P(0) a2I (a为充分大的数)?e(0) 0Pe(0) IC(z 1)z(k)i，计算 Zf(k)和 Uf(k);C(z1)z(k) 0 a ,a,b,

13、,b3利用式1 na 1nb，构造 hf(k) ；hf (k) Zf(k 1), , Zf (k na),uf(k 1), ,uf (k nb)?(k)?(k1) Kf(k)Zf(k)hf(k)?(k1)4利用式 Kf(k) Pf(k 1)hf(k)hf(k)Pf(k 1)h,(k) 11 递推计算 &k)；Pf (k)IKf(k)hf(k)Pf(k1)5利用 e(k) z(k) h(k)4k)和h(k) Z(k 1), , Z(k na),u(k 1), ,u(k nb) 计算 e?(k) ；6根据 he(k) e?(k 1), e?(k nc) 来构造 he(k) ；?(k)7利用Ke(k

14、)Pe(k)?(k 1) Ke(k)砥k) he(k)(k 1)_1Pe(k 1)he(k) he(k)Pe(k 1)h,(k) 1IKe(k)he(k)Pe(k 1)返回第2步进行迭代计算，直至获得满意的辨识结果3工程实例3.1 典型系统建模以某微循环流体系统模型的参数辨识为例.我们已经得到该系统模型的差分 方程形式，取特定点的压力波作为模型的输入，以另一点的压力波作为模型的输 出.由于我们采集的数据是实时的，因此采用在线辨识方法.由于建立的微循环流 体系统模型是一个单输入、单输出的模型，为使参数估计的结果很好地跟踪参数 真值的变化，我们采用渐消记忆的最小二乘法对系统模型参数进行辨识，即强调

15、新数据的作用，贬低老数据的作用.抽象出的SISO系统的差分方程为：z(k) az(k 1) a?z(k 2) hu(k 1) b?u(k 2)(k) 式(1-1)参数取真值为：T 1.376 0.483 0.57 0.42 ,利用MATLA由勺M语言辨识 系统中的未知参数a1、a2、b1、b2。要求：用参数的真值利用差分方程求出z(k) 作为测量值，(k)是均值为0,方差为0.1、0.5和0.01的不相关随机序列。 使用最小二乘算法辨识。3.2 广义最小二乘递推算法的 MATLA的真(程序源代码见附录)考虑仿真对象z(k)= -1.376z(k-1)-0.483z(k-2)+0.57u(k-1

16、)+0.42u(k-2)+v(k)式中，v(k)是均值为0,方差为0.01、0.1和0.5的不相关随机序列。输入信 号采用4阶M序歹1，幅度为1。v(k)u(k)0.57z 1 0.42z 21 1.376 z 1 0.483z 2y(k)z(k)选择如下形式的辨识模型图3.1广义最小二乘法辨识实例结构图其中取 c1=0,c2=0.4结果分析及算法优化由于辨识算法中输入或噪声信号为不相关随机序列，所以每次辨识结果都不 完全相同。但是，在相同输入、相同的噪声、相同的步长条件下，精度大体相同。算法优化方案：(1)使用M序列(具有近似白噪声的性质)为输入信号；(2)增加数据长度去L;(3)减小噪声信

17、号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.12311.37230.0063a20.4830.6884-2.00540.58920.10610.48150.0231b10.570.7916-2

18、.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.3761.3660-0.02801.37500.00791.37290.0005a20.4830.4749-0.07010.48590.029

19、20.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的随机序列，数 据长度取L=30,得到的变化曲线图：下面给出以M序列作为输入，噪声均值为零，方差为 0.01的不相关随机序 列，数据长度去L=300,得到的变化曲线

20、图：J J参考文献1 李言俊，张科，系统辨识理论及应用，国防工业出版社， 2006年2 方崇智，萧德云，过程辨识，清华大学出版社， 2002年3 贾秋玲，袁冬莉，栾云凤，基于 MATLAB7.x/Simulink/Stateflow 系统仿真、分析及设计，西北工业大学出版社， 2006年4侯媛彬，汪梅，王立琦，系统辨识及其MATLAB&真，科学出版社，2004年附录广义最小二乘递推算法的MATLA防真程序源代码：clear % 清理工作间变量L=300; % M 序列的周期y1=1;y2=1;y3=1;y4=0; % 四个移位寄存器的输出初始值for i=1:L;% 开始循环，长度为 Lx1=

21、xor(y3,y4); % 第一个移位寄存器的输入是第三个与第四个移位寄存器的输出 的“或”x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出y(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; % 为下一次的

22、输入信号做准备end %大循环结束，产生输入信号ufigure(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*ze(k-2)+v(k-1);%C(z1)=1, 即取 c1=0,c2=0 end z(2)=0;z(1)=0; % 设 z 的前两个初始值为零 for k=3:301; % 循环变量从3 到 301z(

23、k)=-1.376*z(k-1)-0.483*z(k-2)+57*u(k-1)+0.42*u(k-2)+ze(k-1); %输出采样信号(测量值)end%RGLST义最小二乘辨识c0=0.0001 0.0001 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); % 相对误差

24、的初始值及大小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(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 的值d1=zf(k)-hf1*c0; c1=c0+k1*d1; %求被辨识参数ce1=c1-c0; % 求参数当前值与

25、上一次的值的差值e2=e1./c0; %求参数的相对变化e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列c0=c1; % 新获得的参数作为下一次递推的旧参数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)人2;%求准则函数ee(k)=z(k)-h1*c1;he1=-ee(k-1),-ee(k-2);x=he1*pe0*he1+1

26、; x1=inv(x);k1=pe0*he1*x1;d1=ee(k)-he1*ce0;ce1=ce0+k1*d1;pe1=pe0-k1*he1*pe0;ce0=ce1;ce(:,k)=ce1;pe0=pe1;end %大循环结束c%W识参数变化矩阵%显示被辨识参数及其误差( 收敛 ) 情况%分离参数a1=c(1,1:300); a2=c(2,1:300);b1=c(3,1:300);b2=c(4,1:300);c1=ce(1,1:300);c2=ce(2,1:300);ea1=e(1,1:300); ea2=e(2,1:300);eb1=e(3,1:300); eb2=e(4,1:300);f

27、igure(2); % 第二个图形i=1:300; % 横坐标从 1 到 300plot(i,a1,r,i,a2,k,i,b1,b,i,b2,c,i,c1,b,i,c2,r) %画出 a1 ，a2,b1 ， b2,c1,c2 的各次辨识结果title( 参数变化曲线) % 图形标题figure(3); % 第三个图形i=1:300; % 横坐标从 1 到 300plot(i,ea1,r,i,ea2,k:,i,eb1,b,i,eb2,k:)%画出a1， a2， b1， b2, 的各次辨识结果的收敛情况 title( 误差曲线 ) % 图形标题考虑仿真对象 z(k)+1.5z(k-1)+0.7z(

28、k-2)=u(k-1)+0.5u(k-2)+v(k),其中v(k)是服从正态分布的白噪声N(0,1)。输入信号采用4阶M序歹U。选择如下形式的辨识模型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;m length = r*Np;a=1;for i=1:1:m_lengthY4=X4;Y3=X3;Y2=X2;Y1=X1;X4=Y3;X3=Y2;X2=Y1;X1=xor(Y3,Y4);if Y4=0M(i)=-a;elseM(i)=a;end

29、endfigure;i=1:1:m_length;plot(i,M);%白噪声noise = zeros(1,m length);for i=1:1:m lengthtemp = noise + 0.5*rands(1,m length);noise = temp;endnoise = noise/12;%noise = temp;% parameter of systemn=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5;S U0=0.2;S Y0=0.2;% generate u,y u0=ones(1,m length)*S U0;U = M + u0 + noise;figure;i=1:1:m length;plot(i,U);%formulationy(1)=0;y(2)=0;y(3)=0;Y(1)=S Y