系统辨识—最小二乘法概要

1、最小二乘法参数辨识1引言系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。 现代控制理论中的一个分支。通过辨识建立数学模型的目的是估计表征系统行为 的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入 和输出预测系统输出的未来演变，以及设计控制器。对系统进行分析的主要问题 是根据输入时间函数和系统的特性来确定输出信号。对系统进行控制的主要问题 是根据系统的特性设计控制输入，使输出满足预先规定的要求。而系统辨识所研 究的问题恰好是这些问题的逆问题。通常，预先给定一个模型类卩=M（即给定一类已知结构的模型），一类输入信号u和等价准则J=L（y，yM）（ 般情况下，

2、 J是误差函数，是过程输出y和模型输出yM的一个泛函）；然后选择使误差函数 J达到最小的模型，作为辨识所要求的结果。系统辨识包括两个方面：结构辨识 和参数估计。在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计 这两个方面并不是截然分开的，而是可以交织在一起进行的。2系统辨识的目的在提出和解决一个辨识问题时，明确最终使用模型的目的是至关重要的。 它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。通过辨 识建立数学模型通常有四个目的。 估计具有特定物理意义的参数有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。这就 需要通过能观

3、测到的输入输出数据，用辨识的方法去估计那些参数。 仿真 仿真的核心是要建立一个能模仿真实系统行为的模型。用于系统分析的仿真模型要求能真实反映系统的特性。用于系统设计的仿真，则强调 设计参数能正确地符合它本身的物理意义。 预测 这是辨识的一个重要应用方面，其目的是用迄今为止系统的 可测量的输入和输出去预测系统输出的未来的演变。例如最常见的气象预报，洪 水预报，其他如太阳黑子预报，市场价格的预测，河流污染物含量的预测等。预 测模型辨识的等价准则主要是使预测误差平方和最小。只要预测误差小就是好的 预测模型，对模型的结构及参数则很少再有其他要求。这时辨识的准则和模型应 用的目的是一致的，因此可以得到较

4、好的预测模型。 控制 为了设计控制系统就需要知道描述系统动态特性的数学模型，建立这些模型的目的在于设计控制器。 建立什么样的模型合适，取决于设计 的方法和准备采用的控制策略。3系统辨识的方法经典方法：经典的系统辨识方法的发展已经比较成熟和完善，他包括阶跃响应法、脉冲 响应法、频率响应法、相关分析法、谱分析法、最小二乘法和极大似然法等。其 中最小二乘法（LS）是一种经典的和最基本的，也是应用最广泛的方法。但是，最 小二乘估计是非一致的，是有偏差的，所以为了克服他的缺陷，而形成了一些以 最小二乘法为基础的系统辨识方法：广义最小二乘法 （Gl S）、辅助变量法（IV）、 增广最小二乘法（El, S）

5、和广义最小二乘法（Gl S），以及将一般的最小二乘法与其 他方法相结合的方法，有最小二乘两步法（CORI S）和随机逼近算法等。经典的系统辨识方法还存在着一定的不足：（1）利用最小二乘法的系统辨识法一般要求输入信号已知，并且必须具有较丰富的变化，然而，这一点在某些动 态系统中，系统的输入常常无法保证；（2）极大似然法计算耗费大，可能得到的 是损失函数的局部极小值；（3）经典的辨识方法对于某些复杂系统在一些情况下 无能为力。现代方法：随着系统的复杂化和对模型精确度要求的提高，系统辨识方法在不断发展， 特别是非线性系统辨识方法。主要有：1、集员系统辨识法在1979年集员辨识首先出现于Fogel撰写

6、的文献中，1982年Fogel和Hua ng 又对其做了进一步的改进。集员辨识是假设在噪声或噪声功率未知但有界UBB（Unknown But Bounded）的情况下，利用数据提供的信息给参数或传递函数 确定一个总是包含真参数或传递函数的成员集（例如椭球体、多面体、平行六边体等）。不同的实际应用对象，集员成员集的定义也不同。集员辨识理论已广泛 应用到多传感器信息融合处理、软测量技术、通讯、信号处理、鲁棒控制及故障 检测等方面。2、多层递阶系统辨识法多层递阶方法的主要思想为：以时变参数模型的辨识方法作为基础， 在输入 输出等价的意义下，把一大类非线性模型化为多层线性模型，为非线性系统的建模给出了

7、一个十分有效的途 径。3、神经网络系统辨识法由于人工神经网络具有良好的非线性映射能力、自学习适应能力和并行信息 处理能力，为解决未知不确定非线性系统的辨识问题提供了一条新的思路。与传统的基于算法的辨识方法相比较，人工神经网络用于系统辨识具有以下 优点：（1）不要求建立实际系统的辨识格式，可以省去对系统建模这一步骤；（2）可以对本质非线性系统进行辨识；（3）辨识的收敛速度仅与神经网络的本身及所 采用的学习算法有关； 通过调节神经元之间的连接权即可使网络的输出来逼 近系统的输出；（5）神经网络也是系统的一个物理实现，可以用在在线控制。4、模糊逻辑系统辨识法模糊逻辑理论用模糊集合理论，从系统输入和输

8、出的量测值来辨识系统的模 糊模型，也是系统辨识的一个新的和有效的方法，在非线性系统辨识领域中有十分广泛的应用。模糊逻辑辨识具有独特的优越性：能够有效地辨识复杂和病态结构的系统； 能够有效地辨识具 有大时延、时变、多输入单输出的非线性复杂系统； 可以辨识性能优越的人类控 制器；可以得到被控对象的定性与定量相结合的模型。模糊逻辑建模方法的主要 内容可分为两个层次：一是模型结构的辨识，另一个是模型参数的估计。典型的 模糊结构辨识方法有：模糊网格法、自适应模糊网格法、模糊聚类法及模糊搜索 树法等。5、小波网络系统辨识法小波网络是在小波分解的基础上提出的一种前馈神经网络口，使用小波网络进行动态系统辨识，

9、成为神经网络辨识的一种新的方法。小波分析在理论上保 证了小波网络在非线性函数逼近中所具有的快速性、准确性和全局收敛性等优 点。小波理论在系统辨识中，尤其在非线性系统辨识中的应用潜力越来越大， 为 不确定的复杂的非线性系统辨识提供了一种新的有效途径，其具有良好的应用前 景。1. 设计内容设SISO系统的差分方程为：z(k) Qz(k - 1) a2z(k - 2) =bi(k - 1) b2u(k - 2) (k) 式( 1-1) 参数取真值为：丁二1.376 0.4830.57 0.42】，利用MATLAB的 M语言辨识系统中 的未知参数a1、a2、b|、b2。2. 设计过程2.1 建立系统设

10、SISO系统的差分方程为：z(k) a1z(k -1) a2z(2Hb1u(k -1) b2u(k -2) (k)式(2-1)参数取真值为：”二1.376 0.483 0.57 0.421利用MATLAB勺M语言辨识系统中 的未知参数a1、a2、b1、b2。要求：用参数的真值利用差分方程求出z(k)作为测量值，：(k)是均值为0, 方差为0.1、0.5和0.01的不相关随机序列。选取一种最小二乘算法辨识。2.2 最小二乘简介2.2.1 最小二乘法的概念与应用对工程实践中测得的数据进行理论分析，用恰当的函数去模拟数据原型是一 类十分重要的问题，最常用的逼近原则是让实测数据和估计数据之间的距离平方

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

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

13、)C(z )A(z°) = 1 + 卒一1 +| + an< aB(zj = biz” b2z-2gz叫:C(zT) = 1 + epT + C2Z + CnazfD(zJ 二 dh1 d2z2dnbzf2.2.3准则函数设一个随机序列Nk),k(1,2，丄)'的均值是参数丁的线性函数： ENk)'hT(k戸，其中h(k)是可测的数据向量，那么利用随机序列的一个实现， 使准则函数：L2J(R 二' z(k) hT(k)d7(式 2-2)达到极小的参数估计值二称作二的最小二乘估计。最小二乘格式：z(k)二ht (k：h e(k)，押为模型参数向量，e k为

14、零均值 随机噪声。2.3 广义最小二乘法2.3.1 广义最小二乘数学模型A(z')z(k) =B(z')u(k)v(k)式中，u(k)和z(k)表示系统的输入输出；v(k)是均值为零的不相关的随机序列；且'a(z)=1 +aiz，+a2Z，* +anaza2.3.2 广义最小二乘递推算法如下玫k)=珂k -1) K f (k)Zf (k) - hf (k)珂k - 1)K f(k) = Pf(k-1)hf(k)h«k)Pf(k-1)hf(k) +1$ f(k)=l - K f(k)h«k)Pf(k-1)£(k)二纟(k -1) K e(k)

15、gk) - he(k户?(k-1) K e(k)二 Pe(k -1)he(k) hJ(k)P°(k -1)h°(k) + 1 Pe(k)珥I K e(k)h e(k) PeT)式中h(k)二-乙(k-1),厂乙(k - nju (k -1),，q(k - a) *he(k)=£k-1),-e(k-nc)r貝k)=z(k)- h Yk)致k)233广义最小二乘递推算法的计算步骤:1.给定初始条件邀0) = g(充分小的实向量) Pf(0) = a2l (a为充分大的数)?e(0)= 0Pe(O)"2利用式Zf(k)二 C(z 丄)z(k)Uf(k)二 C(

16、z)z(k)计算 Zf(k)和 u«k)；3利用式0羽，构造h f(k)；n(k)=-乙(k-1),，-z(k-rQ,Uf(k-1),u(k-nJF珂k)卫k-1) K f(k) zjk)- hf(k)?(k-1)4利用式 K f(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 利用 e(k)二 z(k) -h (k)坯k)和h(k)十z(k-1),，-z(k- na),u(k-1), ,u(k- nb)计算e(k);6根据 he(k)二-e(k - 1),-e(k

17、 - nJ来构造 h°(k)；玖k)?(k-1) Ke(k)gk)- he(k)玖k-1)7 利用 K e(k)= Pe(k-1)h e(k) °(k) Pe(1)h °(k) 1 Pe(k)=l - K e(k) h J(k) Pe(kT)返回第2步进行迭代计算，直至获得满意的辨识结果。2.3.4广义最小二乘递推算法的 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的不相关随机序列。输入信 号采

18、用4阶M序列，幅度为1。选择如下形式的辨识模型v(k) 11 +c1z_L+c2 z-2 e(k)+0.57z-40.42z u(k)1 41.376Z二*.483z 丄y(kz(k)图2广义最小二乘法辨识实例结构图其中取 c1=0,c2=0.3. 结果分析及算法优化由于辨识算法中输入或噪声信号为不相关随机序列，所以每次辨识结果都不 完全相同。但是，在相同输入、相同的噪声、相同的步长条件下，精度大体相同。算法优化方案：（1）使用M序列（具有近似白噪声的性质）为输入信号；（2）增加数据长度去L；（3）减小噪声信号v（k）的方差。3.1广义最小二乘递推算法的的MATLA仿真结果及分析（1）、输入选

19、用题目给出的30个随机数，即数据长度去L=30,噪声选用 均值为零，方差分别为0.5、0.1和0.01的随机序列，辨识结果如表 3-2-1表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差表 3-2-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.02160.67900.18210.57990.0510b20.420.62

20、13-2.00130.53610.11610.42180.0234输入M序列，30步，噪声方差0.5时: 错误！未找到引用源。=12.556;输入M序列，30步，噪声方差0.1时:错误！未找到引用源。=2.5822;输入M序列，30步，噪声方差0.01时: 错误！未找到引用源。=0.1706;(2) 、输入均采用M序列，噪声选择均值为零，方差为 0.5、0.1和0.01 的随机序列，辨识步长均为300步，辨识结果如表3-2-2。表中给出了三种情况下辨识参数结果即表中的估计值，估计值与真值的相对误差.表 3-2-2真值噪声方差为0.5噪声方差为0.1噪声方差为0.01估计值相对误差估计值相对误差

21、估计值相对误差a11.3761.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输入采用M序列，噪声方差0.5时：错误！未找到引用源。=135.288 输入采用M序列，噪声方差0.1时：错误！未找到引用源。=28.917 输入采用M序列，噪声方差0.01时：错误！未找到引用源。=2.7374(3) 数据结果分析：

22、输入采用 M序列比采用随机序列得到的辨识效果更 好。噪声均值相等时，方差越大，辨识效果越差，反之，方差越小辨识效果越 好。可以通过增加步长的方法提高辨识精度。下面给出以M序列作为输入，噪声均值为零，方差为 0.01的随机序列，数 据长度去L=30,得到的变化曲线图：回凶下面给出以M序列作为输入，噪声均值为零，方差为 0.01的不相关随机序 列，数据长度去L=300,得到的变化曲线图：Figure 2 XPile EditView Insert Tools Desktop Window Help%D 3 R总溟因匚口4. 参考文献1）侯媛彬,汪梅,王立琦,系统辨识及其MATLA仿真,科学出版社,

23、2004年2）方崇智，萧德云，过程辨识，清华大学出版社，1988年3）贾秋玲，袁冬莉，栾云凤，基于MATLAB7.x/Simulink/Stateflow系统仿真、分析及设计，西北工业大学出版社，2006年4）李言俊，张科，系统辨识理论及应用，国防工业出版社，2006年课程心得上完这门课的第一次课，使我懂得了：辨识就是在输入和输出数据的基础上， 从一组给定的模型中，确定一个与所测系统等价的模型。了解了辨识的一些基本 概念。知道了辨识的内容和步骤。我们要判断一个模型是否和原有系统等价或者 近似，就好像我们要认识一个人一样，必须通过多种方法，方式去测试他。这就 需要我们给这模型充分的最优激励。这个

24、形象的比喻让我记忆深刻。接着我们学习了随即过程的概念及相关知识， 谱密度函数等内容，这些为学 习白噪声，M序列等典型激励提供了基础。有了激励U(k)必然有相应的输出Z(k) 与之对应，通过一批一批的数据 U(k)、Z(k)，通过用Matlab编程，能够算出 模型的参数。最后学习的是模型的辨识方法，辨识方法有很多，但根据模型形式可以分成 两类。一类是非参数模型辨识方法，另一种是参数模型辨识方法。第一种方法又 称经典辨识方法，不必事先确定模型具体结构，可适用于复杂的线性过程。第二 种方法又称现代辨识方法，必须假定一种结构模型，通过极小化模型与过程之间 的误差准则函数来确定模型的参数。这种方法中我们

25、主要学得是最小二乘类方 法。这次作业选择了一个较简单的仿真实例，通过这次作业知道了一些函数和指令 的使用方法。以后我想我会努力用一些时间来继续学习这两方面的知识内容。最后感谢张老师的辛勤教育与谆谆教诲!附录广义最小二乘递推算法的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; %

26、第三个移位寄存器的输入是第二个移位寄存器的输出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; %为下一次的输入信号做准备end %大循环结束，产生输入信号ufigure(1

27、); % 第一个图形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,即取 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*

28、u(k-2)+ze(k-1); %输出采样信号(测量值)end%RGL旷义最小二乘辨识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); %相对误差的初始值及大小ee=zeros(2,300);s=0;%广义最小二乘递推算法的计算步骤f

29、or 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 的值d仁 zf(k)-hf1'*c0; c仁 c0+k1*d1; %求被辨识参数 ce仁c1-c0; %求参数当前值与上一次的值的差值e2=e1./c0; %求参数的相对变化e(:,k)=e2; %把当