系统辨识PPT学习教案

1、会计学12 结束第1页/共259页3第2页/共259页4一、建模的必要性一、建模的必要性课程的核心问题是建模，主要是建模。系统辨识是研究辨识建模的理论和方法。 数学模型的主要用途： 控制理论与控制工程就一直围绕着 建立模型和控制器设计这两个主题来发展， 它们相互依赖、相互渗透并相互发展。1.用来预报实际系统物理量 研究实际系统往往需要事先知道一些物理量的数值，而其中有些量可能无法直接测量或测不准，所以需要建立数学模型来预报。 第3页/共259页53.为了设计控制系统 目前,对被控系统的控制器的设计方法的选取,以及如何进行具体的控制结构和参数的设计都广泛依赖于对被控系统的理解及所建立的被控系统数

2、学模型。2.用于分析实际系统 工程上在分析一个新系统时，通常先进行数学仿真，仿真的前提必须有数学模型。 建模问题在控制器设计中起着非常重要的作用,是设计中首先需要解决的问题；是成功地进行控制器设计的关键之一。第4页/共259页系统的模型一般分物理模型与数学模型物理模型：指用物理、化学、生物等材料构成的用于 描述系统中的关系和特征的实体模型。 模型：就是把系统实体的本质信息简缩成有用的描述形式，数学模型：描述系统中一些关系和特征的数据模型。 控制领域的数学模型就是指能用来描述系统的动态或静态特性和行为的数学表达式或方程。是进行系统分析、预报、优化及控制系统设计的基础。二、模型二、模型是一种简化描

3、述。第5页/共259页71、理论建模法： 通过对系统内在机理的分析，按照已知的一些物理定律导出各物理量关系来建立数学模型。 理论建模法建立的模型称为机理模型。 一般在理论建模中,根据模型应用的目的和精度要求,仅考虑系统中起主导作用的有限的几个因素即可。缺陷：当验前信息不足时，用理论建模法会遇到很大困难。对于比较复杂的过程，必须对机理模型简化，这就使得机理建模与实际过程间有一定的误差。 三、建模方法三、建模方法第6页/共259页8 理论建模通常只能用以建立比较简单系统的模型（白箱问题）。由于许多系统的机理和所处的环境越来越复杂，因此，理论建模法的运用亦越来越困难，其局限性越来越大, 需要建立新的

4、建模方法。 在被建模的装置尚不存在（设计阶段）或虽存在但无法进行实验时，理论建模是取得模型的唯一途径，是验前问题中唯一可行的方法。 理论建模的难点在于对有关学科知识及实际经验的掌握，故不属于课程的讨论范围。 在理论建模方法难以进行或难以达到要求的情况下，系统辨识建模方法就幸运而生。第7页/共259页92、辨识建模法： 对被控系统进行测试，利用观测数据，通过辨识技术去构造系统模型的方法。系统辩识是研究怎样利用对未知系统的试验数据或在线运行数据(输入/输出数据)建立描述系统的数学模型的科学。系统辩识亦称为实验建模方法,它是“系统分析”和“控制系统设计”的逆问题。是现代控制理论的一个分支。第8页/共

5、259页1）完全辨识问题： 完全不了解系统的任何基本特性（定常时变；线性非线性；确定随机等）。 这类问题称为黑箱问题。这是一个极难解决的问题，通常需要对系统作某些主观的先验假设。2）部分辨识问题： 系统的某些基本特性假定是已知的，但不知动态模型的阶次或有关的系数。 这类问题称为灰箱问题。显然比黑箱问题容易解决。 根据对系统事先了解的程度（先验知识）可将辨识问题分成二类：完全辨识问题和部分辨识问题。第9页/共259页11 大部分工程系统及工业过程都属于灰箱问题。通常对系统的结构会有很多了解，因此可推导得系统特定的数学模型。在这种情况下只要定阶和确定模型中的一组参数。从而模型化问题简化为参数估计。

6、因此参数估计是一个最重要的问题。 有效的辨识策略：u尽可能地掌握系统的先验知识，即尽可能地使系统“白化”；有效的辨识方法：“灰箱”方法。将两种方法结合起来，互为补充。u对依然“黑”的部分，用理论建模方法不能确定的部分和参数，采用系统辨识方法。第10页/共259页12系统辨识的框图 对象 输入u(k) 测量 测量 测量噪声 测量噪声 输出y(k) 输出测量值 输入测量值 过程噪声 系统辨识 第11页/共259页13模糊数学创始人L.A. Zadeh1962年 Zadeh从数学的角度定义:辨识就是在输入输出数据的基础上，从一组给定的模型类中，确定一个与所测系统等价的模型。1978年瑞典的李龙(Lj

7、ung)提出 :系统辩识的三个要素数据、模型类和准则。系统辩识是按照一个准则，在模型类中选择一个与数据拟合得最好的模型。拟合的好坏是一个不定的概念，所以要用准则来判别。 3、系统辨识的定义 所谓辨识建模是从实验数据出发，根据辨识的目的以及对过程已有的验前知识，预先给出一个模型类（线性的、非线性的、定常的、时变的、连续的、离散的 ）进行拟合。 第12页/共259页14它是一个迭代过程。大致包括：试验设计，模型结构确定，参数估计和模型验证。四、系统辨识的内容（或步骤）四、系统辨识的内容（或步骤）满意辨识目的及先验知识确定模型结构和准则模型的参数估计模型验证最终模型不满意数据预处理输入输出数据检测试

8、验试验设计辨识的一般步骤第13页/共259页15 大致包括：试验设计，模型结构确定，参数估计和模型验证。1、试验设计 1）选择变量：以提取有效的信息（数据）为目的。首先根据试验对象，确定所要观测的变量。 （u是人为给定的，y是观测的，y的选取不同会改变输出矩阵C的结构和数值。） 通常为得到试验设计前的必要的知识，必须进行一些预备性试验（摸底）。四、系统辨识的内容（或步骤）第14页/共259页16预备性试验：可用一些简单方法（阶跃响应，频率响应等）获得系统的如下信息：主要时间常数（系统频宽，与试验长度有关）允许的输入信号幅度（系统的线性范围）过程的非线性与时变性（有助于模型类的选择）噪声水平（以

9、便用多大的输入，使得观测量有多大的信噪比）变量之间的延迟（滞后环节参数）2）输入信号的选择（阶跃、方波、脉冲、PRBS）。 第15页/共259页17 3）采样速度的选择（要采集数据就有采样速度选择问题）。实际上先采用较短的采样间隔，在数据分析时，可根据需要隔几个取一个数据。 4）试验长度的确定（试验时间问题）。辨识精度与试验时间的长短有关。 2、模型结构确定 根据辨识的目的及对被辨识系统的先验知识，确定系统所属的模型类 模型结构的选择主要取决于应用的目的及精度要求。通常模型精度与复杂性要折衷考虑。第16页/共259页18 常用的模型类： 参数的 或 非参数的 线性的 或 非线性的 连续的 或

10、离散的 确定的 或 随机的 I/O的 或 状态的 时变的 或 定常(时不变)的集中参数的 或 分布参数的 频率域的 或 时间域的 等等。第17页/共259页19 根据系统的空间、时间的离散化情况，模型可分为三类：1）集中参数的连续时间模型：空间变量是离散的， 时间变量连续。如常微分方程，代数方程。2）集中参数的离散时间模型：时、空变量均离散。 如差分方程，代数方程。3）分布参数模型：时、空变量均连续，如偏微分方程。 它可以在空间上离散化，简化成分块集中参数，所以对它的辨识不介绍。 a第18页/共259页203、参数估计 模型结构确定后，其中未知部分就要通过观测数据进行估计。通常未知部分是以未知

11、参数出现，故辨识工作就成了参数估计。4、模型验证 一个模型辨出来后，是否可靠必须进行多次验证。 参数估计的要求就是要辨识出来的模型与实际过程在某种意义下最“接近”。 所以必须有个准则衡量。 通常一个模型用一套数据进行辨识，然后用另一套数据来验证和修改。 第19页/共259页21第20页/共259页22模型可以有不同的形式，不同的模型适于不同的系统。古典辨识方法：采用时域法和频率法来辨识线性系统的传递函数。原则上只适用于SISO线性系统。 SISO系统通常采用传递函数。 MIMO系统通常采用状态空间表达式。 由实验来建立数学模型传递函数，可以为更复杂的系统辨识做预备性实验，它是现代系统辨识的基础

12、，属于连续系统的数学模型的辨识领域。 第21页/共259页23试验信号的选用： 对系统模型的研究方法不同，输入试验信号也相应分成非周期的和周期的两种。用时域法建模：输入信号为非周期的。 主要采用阶跃和方波（近似脉冲）函数。用频域法建模：输入信号用周期的。 主要用正弦波，二进制周期函数。它们又分为单频和多频（组合正弦波及周期方波） 第22页/共259页242 2 1 1 过渡响应法（时域法）过渡响应法（时域法） 采用非周期试验信号，通过系统的动态响应研究系统的模型。 一、非参数模型的辨识 在时域中建立线性系统非参数模型时，用很简便的方法就可得到脉冲响应曲线，阶跃响应曲线、方波响应曲线或它们的离散

13、采样数据表。 对于线性系统，脉冲响应，阶跃响应和方波响应之间是可以相互转换的。 脉冲响应：可以采用幅值相当大，宽度很窄的方波来近似函数 。 第23页/共259页25二、由阶跃响应曲线辨识传函1、试探法 工业中常用的模型类：（即便是高阶系统也用低阶模型去逼近） 由非参数模型转变成参数模型，包括确定传函的结构及参数。 1TsK) s (Gse1TsK) s (G1Ts2sTK) s (G22s22e1Ts2sTK) s (G 先观察试验所得响应曲线的形状特征，据此判断，从模型类中确定一种结构。然后进行参数估计，最后验证数据拟合程度，反复多次，直至误差e(t)最小（验证数据拟合可只取若干点）。 第2

14、4页/共259页261）若阶跃响应曲线特征为：0)t (ymax)0(yconst)(y曲线逐渐上升到稳态值：1TsK) s (G可采用结构： 待估参数为：K，T稳态增益： 0U)(yK)(y) t (y) t (y1)(y将试验曲线标么化，即 , 第25页/共259页27要确定 T ，只要一对观测数据：y*(t1)，t1 则标么化后响应：Tte1) t (y可得： )t (y1 ntT11Tte1) t (y由若取 y*(t1) = 0.63 ，则 T = t1验证数据拟合如何， 可在 t=T/2 和 t=2T 二点进行：39. 0)2T(y87. 0)T2(y若拟合不好，则应另选模型结构类

15、。1第26页/共259页28待估参数为：K，T， 稳态增益： 0U)(yK将试验曲线标么化，即)(y) t (y) t (y1)(y2）实验曲线是一条S形非周期曲线 ）可选用模型类： se1TsK) s (G te1t 0) t (yTt则为了确定 T 和，必须将两个坐标值（观测值）代入，122Tt21Tt1tt y e1)t (yy e1)t (y21则第27页/共259页29两边同取对数得：)t (y1 nTt)t (y1 nTt2211)t (y1n)t (y1n)t (y1nt)t (y1nt)t (y1n)t (y1nttT2121122112根据两对观测值 y*(t1) 和 y*(

16、t2) ，可求出 T 和 。 122Tt21Tt1tt y e1)t (yy e1)t (y21第28页/共259页30若选y*(t1) = 0.39，y*(t2) = 0.63，则Tt2Tt21模型验证：0.87)(ty T2t0.55)(ty T8 . 0t0)(ty t554433tT/20.8TT2Ty*(t)0.390.550.630.87Tte1) t (y由则2112tt2)tt (2T122Tt21Tt1tt y e1)t (yy e1)t (y21第29页/共259页31待估参数为：K，T， 究竟选一阶惯性带延时的模型结构，还是选二阶模型，事先无法确定，完全看两种模型与试验曲

17、线拟合程度，哪个精度高，选哪个。 由于大多数工业过程的试验曲线是过阻尼的，即1， 只讨论此种情况，而1的传函辨识比较麻烦。 ）也可选用模型类： 1Ts2sTK) s (G22 S形曲线本身就说明是过阻尼（ 1）。 若 0 且为实数。 21212121T第31页/共259页33)s)(s()s(G2121代入可得： t121t12221ee1) t (y它的单位阶跃响应为： 改写为：t121t12221ee) t (y1令 2 = 1 ( 1 ) ，代入上式得：e11 e1e11e1) t (y1t)1(ttt1111两边同取对数得：e11 nt1n)t (y1 nt)1(11211T第32页/

18、共259页bktt1n)t (y1 n 1可见，当t 时，是一条直线。斜率：k = 1 ，截距： 1nb 在坐标纸上，根据数据y*(t)，画出t较大时)t (y1 n的图形，作其渐进线，即可得斜率 k 和截距 b 。bbb1e111eek则可得： 2 T， 当t 时,0 e11 nt)1(1第33页/共259页35若用常用对数，则当t 时：t4343.01gtge1g)t (y1 g11则：1gb4343. 0k1bbb1101111010k303. 2缺点：计算G(s)时采用的点都是 t 较大时的点，而当 t 较大时，往往 1y*(t) 的值较小，这就会产生较大的误差。bbtk第34页/共2

19、59页362、Laplace 变换的极限定理法（终值定理法） 利用 Laplace 变换的极限定理，由非参数模型的单位阶跃响应，求参数模型传递函数。它克服了试探法需选择模型类的不足，但它仅适用于下述一种模型类。 设线性SISO定常系统的传函结构为： 1sasasaK) s (G11n1nnn特点：系统只有极点、无零点。 第35页/共259页37当输入u(t) = 1(t) 时，输出 y(t)为：终值定理为：) s (Ysim) t (yim0sts1)s (G)s (U)s (G)s (Y对于阶跃响应: 代入上式得： )(yK)s (Gim) t (yimK0st0 K0 第36页/共259页

20、38 在已存在的系统G(s)的基础上构造一个新系统G1(s)，当输入 u(t) = 1(t) 时，其单位阶跃响应为： t001d)(yK)t(y（y1(t)与 y(t)的关系） 求G1(s)的稳态增益K1 :0011dt)t (yK)(yK K1 当输入u(t) = 1(t) 时，输出 y(t)为：第37页/共259页39 G1(s) 求拉氏变换：s1) s (G) s (Ys) s (GsKs1 )(yKLs1d)(yKL)t (yL1100t001)s(GKs1)t (1L)t (yL)s(G011求得 G1(s) 与 G(s) 的关系: 当输入u(t) = 1(t) 时，输出 y(t)为

21、：t001d)(yK)t(y第38页/共259页40再利用终值定理可求得G1(s)的稳态增益K1：) s (Ysim) t (yimK10s1t1 a1 )s (GK(s1im) s (Gim00s10s1sasasaKKs1im11n1nnn00s1sasasaKKsaKsaKs1im11n1nnn010nn00s10aK当输入u(t) = 1(t) 时，输出 y(t)为： 把y1(t)定义成与y(t)有联系，当然G1(s)也与G(s)有联系，而输入均为u(t)=1(t)，再利用终值定理求G1(s)中的参数，从而也就求出G(s)中的参数。 求得G1(s)的稳态增益K1与G(s)中的参数关系。

22、第39页/共259页41 K2 同理，在系统G1(s)基础上构造一个新系统G2(s) 。t0112d)(yK) t (yG2(s)的单位阶跃响应为： 01122dt)t (yK)(yK求拉氏变换：)(yKLs1d)(yKL)t (yL11t0112s) s (GsKs111s1) s (G) s (Y22当输入u(t) = 1(t) 时，输出 y(t)为：第40页/共259页42再利用终值定理可得：20112aKaKK a2 ) s (Ysim) t (yimK20s2t2)s (GK(s1im) s (Gim110s20s G2(s) )s(GKs1)s(G112求得 G2(s) 与 G1(

23、s) 的关系: 当输入u(t) = 1(t) 时，输出 y(t)为：第41页/共259页43当输入u(t) = 1(t) 时，输出 y(t)为： Kr 同理，在系统Gr-1(s)基础上构造一个新系统Gr(s) 。t01r1rrd)(yK) t (yGr(s)的单位阶跃响应为： 01r1rrrdt)t (yK)(yK再用终值定理，由数学归纳法可得： KKn),1,2,(r aK) 1(aKaKK0r01r22r11rr ar 第42页/共259页44KKn),1,2,(r aK) 1(aKaKK0r01r22r11rr特点：）每求一次Ki，要计算一次面积，所以计算量大，而且误差随着积分次数增大而

24、增大。 故仅适用于低阶模型的辨识。 ）使用过程受到一定的限制，仅适用于特定的模型结构（即传函G(s)只有极点，而没有零点的情况）。 由上述（n+1）个方程可求出（n+1）个待估参数： K，a1， ，an第43页/共259页45K1的物理意义：由 可知，001dt)t (yKKK1为阴影部分的面积。 （几何意义）0001dt)t (y) t (y dt)t (yKK将 G(s) 改写成极点形式：n1ii)1sT(K)s (G显然： n1ii1Tan1ii101TKaKK所有时间常数之和系统增益1K （物理意义） 1第44页/共259页46三、由脉冲响应曲线辨识传函 1、矩法 脉冲响应g(t)可由

25、单位阶跃响应微分后求得，也可用窄方波响应来近似。（方波宽度 T 时， g() 0 02T)g(T)g( 对于第一个周期的激励而言， （时间 在 0T 区间）有： )(RK1)(g uy第86页/共259页88 采用具有上述二性质的输入信号后，即可保持采用白噪声信号所具有的优越性，又可以解决 Ruy() 的积分时间太长的问题，理论上只要在一个周期 T 内积分就可以了。 )(RK1)(g uy 思路： 维纳霍甫方程 解决了抗干扰问题 引出了解积分方程难的问题 采用白噪声作为u(t)解决之，使得Ruy()与 成比例 尚存在求Ruy()积分时间长问题 采用周期性，近似白噪声伪随机信号解决之，并仍保持R

26、uy()与 成比例)(g )(g 第87页/共259页89一、M序列产生的方法及性质：随机地掷一枚硬币的随机试验，结果：正面：+1 ；反面：1 反复试验 得到以+1，1两元素组成的随机序列u(k)。 当实验次数N相当大时，该序列u(k)具有以下两性质： 序列中+1与1出现的次数几乎相等；（ Eu=0 ） 随机序列的自相关函数 Ruu(0) = max， 离开原点时， Ruu() = 0 。（Ruu ） 显然该序列接近于白噪声，最好它还应该是一个周期序列，在一个周期内具有上述白噪声性质。它在一个周期内观测时是一个随机信号；若观测时间很长时是一个周期信号。 由于序列只有+1，1两元素，称为伪随机二

27、位式序列 PRBS（Pseudo Random Binary Signal）序列 (它有规律性，故称伪随机，且可以人为产生和复制。)第88页/共259页90n = 4 ，k = 2 （初态为：1111）时， 码数 NP = 6n = 4 ，k = 3 （初态为：1111）时，M序列 最大长度的伪随机二位式序列。 由 n 个双稳态触发器顺序组成 n 级移位寄存器，将第 k 级与第 n 级状态“异或”后，反馈到第一级输入端。 究竟k选哪一级呢？它将影响输出的性质。若k选择合适，将得到一个M序列。（初态不能为全零，否则输出总是零） 码数 NP = 151、M序列的产生第89页/共259页91一个 n

28、 级移位寄存器的输出序列的最大长度 = ？ 除各级全0的状态外，共有（2n 1）种不同的组合状态。 若 NP = 2n 1，则该序列为最大长度序列或M序列。 当 n 12 时，大约半数序列要用二级反馈产生，其他的则要用 4 级反馈来产生M序列。2、M序列的性质1）是一个确定的周期性序列，它的周期长度 NP = 2n 1 2）一个周期内。“0”状态比“1”状态少1个。 （避免出现全“0”状态） “1”状态： 21N222p1nn“0”状态： 21N12122p1nn第90页/共259页923）若将序列中相邻状态不变的那一部分长度称“游程”（或“段”），则在一个周期内的游程总数为 m 。 不允许

29、n 个全零状态， 只有一个 n 个码为全“1”1)2m(m1n1ii 游程总数 21N2mp1n1211)211(21m1n12mm1n 如：n = 4 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0m = 8根据概率论可知，若游程总数为 m ，则： 个个 1 ) 1n(i1 2m i游程长度为 i 的为： 游程长度为 n 的为： 第91页/共259页934）移位相加性： 若将一个M序列与将其延迟了r个码以后的序列，按模2加法原则相加，所得的新序列还是M序列，不过延迟了q个码，r、q均为整数，且1 r，q Np 1 。， 1 r，q Np 1 )qk(u) rk(u)k(u例：5）

30、M序列具有近似离散的白噪声性质。下面将详细讨论M序列的自相关函数和功率密度谱。 b第92页/共259页943、M序列的相关函数和功率密度谱：1）相关函数： 定义:“1”状态的逻辑电平为“a”，“0”状态的逻辑 电平为“+a”（为负逻辑关系，反之结论一样）。 由Ruu()定义可知：T0Tuudt)t (u) t (uT1im)(RpT0pndt)t (u) t (unTnimTp=Np ppnT0pnTdt)t (u) t (unT1impT0pdt)t (u) t (uT1第93页/共259页95(1) 离散情况( = / ): 其中： = / , 设 是 的整数倍，此时 Ruu() 取值于

31、= 0，1，2，Np1Ruu()写成离散形式为： 1N0kpuup)k(u)k(uN1)(R 2a)k(u)k(u 显然模2乘法的结果与模2加法的结果在逻辑上是完全一样的，都为异或关系。 即： )k(u)k(u同号码 +a2（0） u(k)与u(k + )码的电平符号相同。异号码 a2（1） u(k)与u(k + )码的电平符号相异。 （同号码个数）（异号码个数）p2uuNa)(R1N0kpp)k(u)k(uN1第94页/共259页96 Na21N21NNa)(Rp2ppp2uu =1，2，Np1 当 =1，2，Np1 时： u(k)u(k + )在逻辑状态上相当于原序列u(k)与另一延迟序列

32、u(k+)按摸2加法原则相加。根据M序列的移位相加性质可知，所得的结果在逻辑状态上仍是一个M序列。 同号码个数 = 新序列“0”状态个数 =异号码个数 = 新序列“1”状态个数 =21Np21Np当新的序列为“0”状态时，说明u(k)与u(k + )是同号当新的序列为“1”状态时，说明u(k)与u(k + )是异号第95页/共259页97)Na(N1)k(u)k(uN1)0(Rp2p1N0kpuup 当 = 0 时： ( 同一个M序列自乘 ) p2uuNa)(R =1，2，Np1 2uua)0(R第96页/共259页98（2）连续情况（ 不是 的整数倍）：) 1N,1, 0 ( ) 1( pp

33、T0puudt)t (u) t (uN1)(R（为平均面积值） 为一个周期内曲线 u(t)u(t+) 所围成的面积。pT0dt)t (u) t (u即：一个周期内曲线 u(t)u(t+)所围成的： 正面积 负面积第97页/共259页99 当 = 0 时： ( 0 )M序列每出现一次状态转换，积分将出现一个负面积( a2)2a 由M序列的性质3可知，在一个周期内：状态转换的次数 M序列游程总数 21Np负面积)a(21N2p正面积)a (21NNa2pp2第98页/共259页100负面积 )a(21N2p正面积 )a(21NNa2pp2pT0puudt)t (u) t (uN1)(R )a(21

34、N)a(21NNa N12p2pp2p N1N1 a)(Rpp2uu 负面积正面 积 N1p第99页/共259页101 N1N1 a)(Rpp2uu它是 的线性函数，因此可确定其两点：当= 0 时，Ruu()= a2 ; 当 = 时，Ruu()= p2Na Ruu() = Ruu(-) 为偶函数， 在 内，Ruu() 为一个 波 。 第100页/共259页102 当 Ts（系统的过渡过程时间），即： 当 t (Np - 1) 时， g(t) 0 维纳霍甫方程的离散形式为： 1N0kuuuyp)k(R)k(g)(R1N0puup)(u)k(uN1)k(R1N0kpuyp)k(u)k(yN1)(R

35、其中：3 3 用用M序列序列第122页/共259页124在离散情况下M序列的Ruu()为： 1N0k2uu1N0k1uu1N0kuuuyppp)k(R)k(g)k(R)k(g )k(R)k(g)(R) 1N( ,1,2,k Na 0k a )k(RPp22uu1N0kp2pp2p)k(gNa)(gN)1N(a 它是由Ruu2()产生的 令1N0kp2p)k(gNaC第123页/共259页125 脉冲响应函数 g(k) 有界， C 为有界常数，且 C 0 1、作图法C)(R 0)(g uy则时，当)(gN) 1N(app2)(gN) 1N(app2将 Ruy() 上移C 就可以得到： - C 一

36、般可以通过对Ruy() 的稳态值的目测得到。)(g 从而得到： )k(gNa)(R)1N(aN 1N0kp2uyp2pp )C)(R()1N(aN)(guyp2p第124页/共259页1262、解析法)(g 通过精确计算公式得到：C)(gN)1N(a)(Rpp2uy1N01N0pp21N0uypppC)(gN)1N(a)(RC)k(gNa1N0kp2p )k(R)(R)1N(aN)(g 1N0kuyuyp2pp方程两边同求和，可得：1N0kp2p)k(gNaC1N0kp2p1N0pp2pp)k(gNaN)(gN)1N(a 第125页/共259页127二、估计量 的统计特征)(g )(R)k(g

37、Na)(gN)1N(auy1N0kp2pp2p1、 是无偏的，即)(g )(g)(g E方程两边同取均值，可得：)(RE)k(gENa)(gEN)1N(auy1N0kp2pp2p)k(u)k(yN1E1N0kpp)k(u)k(n)k(zN1E1N0kpp)(R)(REunuz第126页/共259页128)(R)(RE)(REunuzuy )k(u )k(u)(gN1E1N0k1N0ppp 此时，又回到理论的维纳霍甫方程，所以： )k(gNa)(gN)1N(aE)(RE1N0kp2pp2uyp=01N0kp2pp2p)k(gNa)(gN)1N(a此时g()已经是一个确定的量了！)k(u )k(u

38、N1)(gE1N0kp1N0pp)(R)(gE1N0uup )k(u)k(zN1E1N0kpp第127页/共259页1291N0kp2pp2uy1N0kp2pp2pp)k(gNa)(gN) 1N(a )(RE )k(gENa)(gEN) 1N(a )(g)(gE 2、 是一致估计量，即)(g 1)(g E)(g P)(R)k(gNa)(gN)1N(auy1N0kp2pp2p方程两边同取二阶原点矩，可得：)(RE )k(gNa)(gN) 1N(aE2uy21N0kp2pp2p0)(g D第128页/共259页130)(g)(gE 22同证明无偏性一样，可得：)(RE )k(gNa)(gN)1N(

39、aE2uy21N0kp2pp2p)(g)(gE 无偏性 ：)(g)(g)(g 2)(g E)(g E)(g E)(g D222可得 的方差为:)(g 0)(g)(g2)(g222自然说明了估计量 的有效性。0)(g D)(g )k(gNa)(gN)1N(aE21N0kp2pp2p21N0kp2pp2 )k(gNa)(gN)1N(ap是一致估计量。)(g 第129页/共259页131三、提高估计精度的方法的估计精度取决于Ruy()的精度。)(g 1、提高采样速率提高Ruy()的精度：取采样周期To= / (=14),用更多的y(t)数据计算Ruy()。2、采用多个周期的M序列输入r+1个周期的M

40、序列，测得r个周期的y(t)计算Ruy()。1N,0,1,2, )k(u)k(yNr1)(Rp1Nr0kpuyp通常取 r = 144第130页/共259页132四、计算 的方法（采用多个周期）)(g 1、一次完成法（离线计算法）定义：) 1N(g ) 1 (g ) 0 (g gp) 1N(R ) 1 (R ) 0 (R Rpuyuyuyuy1N,0,1,2, )k(R)(R)1N(aN)(g p1N0kuyuyp2pp1N,0,1,2, )k(u)k(yNr1)(Rp1Nr0kpuyp )k(R)(R2)k(R)1N(aN1N1kuyuy10kuyp2ppuyp2pR 211121112)

41、1N(aNg上式写出向量矩阵形式：第131页/共259页133) 1rN(y ) 1 (y ) 0(y )NrN(u) 2N(u) 1N(u) 2rN(u) 0(u) 1(u) 1rN(u) 1 (u) 0(uNr1R ppppppppuyUNp x r NpYr Np x 1YUNr1p R 211121112) 1N(aNguyp2pYU 211121112) 1N(ar1p2 1N,0,1, )k(u)k(yNr1)(R ) 1N(R ) 1 (R )0(R R p1Nr0kpuypuyuyuyuyp第132页/共259页1342）需要输入r+1个周期的u(k)：u(-Np+1) u(r

42、 Np-1) 。特点： 1）一次离线求出 （=0,1，Np-1)。)(g 3）精度要求较高时，Ruy()的计算精度要高，r的数目要大， 所以数据存储量大。4）不是递推公式，无法在线辨识。 YU 211121112) 1N(ar1g p2第133页/共259页1352、递推算法 设已获得M对I/O数据，且M （ = 0,1，Np-1 ， M Np-1），即必须先观测至少一个周期。 1N,0,1,2, )k(u)k(y1M1)M,(RpM0kuyu(My(M) )k(u)k(yM1M1M11M0k）uyuypM1 R( ,M)R( ,M1) y(M) u(M M1M1 (0,1,2,N1) ） R

43、uy()的递推公式，即全部的M对I/O数据的Ruy(,M) 可以用过去的(M-1)对I/O数据算得的Ruy(,M-1)和第M次观测的最新数据y(M)和u(M)递推地计算出。 第134页/共259页136可得向量矩阵形式：令)M, 1N(R )M, 1 (R)M, 0(R RpuyuyuyM,uy) 1NM(u ) 1M(u)M(u UpMM,uyp2pMR 211121112)1N(aNg Uy(M)1M1R1MMR M1Muy,M,uyuyuypM1 R( ,M)R( ,M1) y(M) u(M M1M1 (0,1,2,N1) ）第135页/共259页137M,uyp2pMR 2111211

44、12) 1N(aNgUy(M)1M1 R1MM 211121112) 1N(aNM1Muy,p2p 可以通过过去的 和最新的数据y(M)和u(M)在线求得，随着I/O数据的增加， 的精度不断地提高。Mg 1Mg Mg Uy(M)211121112) 1N(aN1M1 g1MM g Mp2p1MM递推公式第136页/共259页138 用逆重复M序列辨识线性系统的脉冲响应函数，比采用M序列所得的算式更简单。 pT0l ly ldt) t(R) t (g)(RpT02 l l1 l ldt)t(R) t(R)t (gpT02 l ldt) t(R) t (gC 对于周期为2Tp的逆重复M序列，仍要求

45、Tp满足 (Tp- ) Ts 。可得维纳霍甫方程为： Cdt) t(R) t (gpT01 l l第137页/共259页139 Ts（系统的过渡过程时间）。工程上一般取： Np = (1.21.5) Ts / 在理论推导时，假设：当 t (Np-1) 时，g(t) 0 一般取：既保证系统的线性，又不超出设备允许公差的最大幅值 a 。第144页/共259页146五、计算Ruy() 在生产现场做试验，一般是在系统的正常工作状态 uo上再附加一个PRBS 输入 u(t) 。 当 t Ts 后，系统的非零初始条件将消失， y(t) 中的零输入响应消失了，只剩下由 u(t) 激起的强制响应，此时 y(t

46、) 已经是一个平稳随机过程了。 所以，必须对系统先加一个周期 u(t) 的预激励，从第二周期开始再量测 I/O 数据，用以计算Ruy() 。b系统的实际输入： u*(t) = uo+ u(t) 系统的实际输出： y*(t) = yo+ y(t) 所以，计算 Ruy() 时必须从实际输出 y*(t) 中将稳态数据 yo 除去。第145页/共259页147 设MIMO系统:输入:u1(t),u2(t),uJ(t) 输出:y1(t),y2(t),yI(t)JIIJ1IJ221J111)(g)(g)(g)(g)(g)(g)(G) t (u) t (u) t (u) t (UJ21G() t (y) t

47、 (y) t (y) t (YI21为了辨识多变量系统 G() ，需要辨识IJ个子系统： gij() (i=1,2,I ; j=1,2,J)第146页/共259页148式中T1 是所有子系统 gij()中最大的系统调节时间，即： 0)(g T T maxTijssj , i1时，J1jT0jiji1d)t (u)(g) t (y第i个输出：显然，yi(t) 与 J 个子系统的脉冲响应有关。 设uj(t) (j=1,2,J) 均为周期性试验信号，且周期T T1 ,与SISO系统相似，可得MIMO系统的维纳霍甫方程：多变量系统卷积公式为向量矩阵形式:1T00d)t (U)(Gd)t (U)(G)

48、t (Y1T0uuuydt) t(R) t (G)(RIJ IJ JJ第147页/共259页149JIyuyuyuyuyuyuT0J1I1)(R)(R)(R)(R)(R)(Rdt)t (u)t (u) t (y) t (yT1IJI12J211J11T0Tuudt)t (U) t (UT1)(R式中:T0Tuydt)t (U) t (YT1)(RJJuuuuuuuuuuuuT0J1J1)(R)(R)(R)(R)(R)(Rdt)t (u)t (u) t (u) t (uT1JJJ12J211J11第148页/共259页150第j个输入、第i个输出的维纳霍甫方程：1Jj1jijT0uuuuiJ1

49、iyudt) t(R) t(R) t (g) t (g)(Rdt) t(R) t (gdt) t(R) t (g11mjjjT0Jjm1mT0uuimuuij 第 j 列、第 i 行自相关函数若uj(t)（j=1,2,J）为互不相关，即：)J, 2 , 1m , jm( 0)(Rmjuu 且 uj(t) 又最好都是 M 序列或逆重复 M 序列 ,则可以用类似前面 SISO 系统的相关法求得：)()(Rjjuu) t (gj i由于存在第二项互相关函数，使得求 产生困难。) t (gj i=0互相关函数第149页/共259页151 取 uj(t) = M 序列， um(t) = Um = Con

50、st ( m j )1、逐个试验法) J , 2 , 1m , jm( C)a(TUdt)t (uTU)(RmpmT0jpmuupmjI),1,2,(i C dt) t(R) t (g)(R 1jjijT0uuijyu缺点： 1）若输入多，J 较大时，则试验时间较长。 2）各个输入之间的交叉作用不能充分反映出来，误差较大。根据 uj(t) 和 y1(t), y2(t),yI(t) 可求得 G() 的第 j 列：)(g, )(g, )(gIjijj1 继续改变 uj(t) , 让 j = 1, 2,J ，可求得 G() 的全部 IJ 个元素：J),2, 1j I;,2, 1i ( )(gij“1

51、”比“0”的状态多1个第150页/共259页1522、联合试验法 为了消除式中的第二项互相关函数，最好的办法是使得: uj(t) (j=1,2,J) 为互不相关，即在一个T周期内：) jm( 0)(R)(K)(Rmjjjuuuu2）当 J 2 时： 利用 M 序列的移位相加性质，可以从一个 M 序列得到J 个 M 序列，它们的相位依次后移了Np J 个码。1）当 J = 2 时： 在 2Tp 内：M序列与逆重复M序列是不相关的。 所以，两个输入可以分别采用M序列和逆重复M序列。同时可求得 G() 的全部 IJ 个元素：J),2, 1j I;,2, 1i ( )(gij第151页/共259页15

52、3)T 1J t (u)JN 1J t (u) t (u JT JN T其中： )T2t (u)JN2t (u) t (u )Tt (u)JNt (u) t (u ) t (u) t (u JpJppJJp3Jp21pjT01jjp1juudt)t (u) t (uT1)(RpT0Jjjpdt)Tt (u) t (uT1对于 M 序列：) 1N( ,2 0Na 0 a )(RPp22uu第152页/共259页154ppjT0JjjpT01jjp1juudt)Tt (u) t (uT1dt)t (u) t (uT1)(R 取 0)(g T T maxTTijssj , i1J时，在小于一个 TJ

53、 内：uj(t) 与 uj+1(t) 以及其它的 um(t) 都不相关。第153页/共259页155 上述两种方法仅适合 J 4 时的系统，否则试验 时间会太长，可能破坏 y(t) 为平稳过程的假设。 只需要一个Tp 就可求得 G() 的全部 IJ 个元素：J),2, 1j I;,2, 1i ( )(gij J 试验时间增长 。 逐个试验法：每个 uj(t) 都要试验两个 Tp（预激励 + 试验） 0)(g T T maxJTJTijssj , i1p时，联合试验法：求 Ruu 和 Ruy 时要试验两个 Tp ，而且：5第154页/共259页156第155页/共259页)nk( C) 1k(

54、C)k( )nk(uB) 1k(uB )k(uB )nk(y a) 1k(y a)k(y n1n10n1)k( e)nk(uB)k(uB)nk(y a) 1k(y a)k(y n0n1（过程噪声） R)k( e )nk( C) 1k( C)k()k( e mn1令：mmin1iiim1rmjn0jjj111in1iii1RC zCI)z(CRB zB)z(BRa za1)z(a )k()C(ze(k) 1)k()z(C)k(u)z(B)k(e)k(u)z(B)k(y)z(a 1111引入延迟算子 z-1 ： ) 1k(x z)k(x 1)白噪声 (R)k( ,R)k(y ,R)k(u其中：mm

55、r157第156页/共259页158CARMA模型（可控自回归滑动平均模型） 1、CARMA： （Controlled Auto-Regression Moving Average）对于SISO系统：)k()z(C)k(u)z(b)k(y)z(a 1112、CAR（可控自回归）： )n,1,2,i ( 0Ci)k()k(u)z(b)k(y)z(a 113、ARMA（自回归滑动平均）： )n,1,2, 0j ( 0bj)k()z(C)k(y)z(a 11第157页/共259页1596、静态自回归：n = 0 （回归分析）5、MA（滑动平均）： )n,1,2, 0j ( 0bj)n,1,2,i (

56、 0ai)k()z(C)k(y1)k()k(ub)k(ub )k(y r14、AR（自回归）： )n,1,2,0j ( 0bj)n,1,2,i ( 0Ci)k()k(y)z(a 1第158页/共259页160一、最小二乘法（ L S ） SISO的随机差分方程描述（CAR模型可控自回归 ）： )n,1,2,i ( 0c , 0bi0（白噪声）)k()k(e )k( e)nk(ub) 1k(ub)nk( y a) 1k( y a)k( y n1n1其中 u(k) 和 y(k) 为实际测量的 I/O 序列， e(k) 为过程噪声（观测噪声）。 假设 e(k) 为独立同分布的随机变量序列，具有零均值

57、和方差2 。 上述方程可写成：)k( e)nk(ub) 1k(ub)nk( y a) 1k( y a)k( y n1n1待估未知参数 2n 个： b,b,a,a n1n1T第159页/共259页161 为了估计 2n 个参数，必须观测 N+n 次，N 2n，从而得到 N 个方程组。定义：1Nx)Nn(y)2n(y) 1n(yY1Nx)Nn( e)2n( e) 1n( een2NxTTT)N(u) 1Nn(u)N(y) 1Nn(y) 2(u) 1n(u) 2(y) 1n(y) 1 (u)n(u) 1 (y)n(y)Nn() 2n() 1n()nk(u ) 1k(u )nk(y ) 1k(y )k

58、(T每一个观测方程：Nn,1,nk )k(e)k()k(yT)k( e)nk(ub) 1k(ub)nk( y a) 1k( y a)k( y n1n1第160页/共259页162每一个观测方程：Nn,1,nk )k(e)k()k(yT可得观测方程组的向量矩阵方程：eY最小二乘的标准格式 特点：输出关于参数是线性的。 最小二乘估计准则： 在最小二乘模型类（CAR模型）中，找出这样一个模型，在这个模型中，系统参数向量 的估计量 ，使得性能指标函数（标量函数）： )Y()Y( ee ) in(eJTTN1i2因为 J 是 的二次函数，所以 J 存在极值。 min第161页/共259页163 极小化估

59、计准则 J 的必要条件：正规方程式 T eeJT)Y()Y(0)Y(2T从而得：YTTYT1T）（解的表达式 T)Y(J)Y(T)Y(2)(0Ja估计量 的解： 第162页/共259页164 极小化估计准则 J 的充分条件： TJ为正定的。 T TT )Y(2(J为实对称阵为实矩阵 T 若T是正则的(非奇异)， 则0)(Ti 0(A)中至少有一个(A)为奇异阵A iiT是准正定的0)(Ti0)(TiT是正定的满足 J 极小的充分条件。T2)Y(2)(JT第163页/共259页165 因为 T 与 无关（为常矩阵），所以 J 只有一个局部极小值存在，当然也是全局极小值。 故最小二乘估计量 是唯一

60、的。 Y 残差： LS 的估计准则实际上是: 使得残差平方和为观测误差平方和的极小值。)Y()Y()(JTT)J(min eemin)Y()Ymin(TT 由于在 中， 是在取好足够数据后一次计算出来的，所以称之为一次完成估计式。Y)(T1T第164页/共259页166二、递推最小二乘法 递推算法可以减少内存存贮量和计算量，同时还可以实现在线辨识。在一次完成法中： I/O数据越多 结果精度越高 存贮量增大，计算量增大 设观测 N+n 对 I/O 数据后，获得参数最小二乘估计为)N()N(Y)N()N()N()N(T1T其中：n2NTTT)Nn()2n() 1n()N(1N)Nn(y)2n(y)