1、（2003 年教案） 辨识与自适应 第三章 1第三章 参数估计的改进算法31 最小二乘方法的改进 辅助变量法（I V法 Instrument Variable ）辅助变量法是在模型误差为相关噪声的情况下，通过引进辅助变量矩阵，对线性最小二乘估计的一种改进。一、参数的辅助变量估计考虑模型 式（3 - 1- 1）其中：当进行了 k = 1-n, 2-n, , 0, 1, 2, , N 共计（ N+n）次采样，得到N 个方程：用矩阵表示成y N = N + e N 式（3 -1- 2）其中 e N = e(1),.,e(N)T最小二乘估计为： LS =（ N T N） - 1 N T y N)(),
2、.1(),(),.1( ,.,.22 nkuknykybaTknk )()(kekyTk（2003 年教案） 辨识与自适应 第三章 2设（ N T N）满秩，将过程模型式（2-3-2）带入上式，得出： LS =（ N T N） - 1 N T（ N + e N）式（3-1- 3）有：和令： = E k kT2n2n 式（3 -1- 4）和 e= E k e(k) 2n1 式（3 -1- 5）可以证明：当N 时式（3-1- 6）式（3-1- 7）依据 Frechek 定理： 若随机序列 xK ， eNkpKTke1)(NTNTLS e11)(11keNekTnkTkTNTN 121,2, （20
3、03 年教案） 辨识与自适应 第三章 3则 f(xK) f( ) 故有： LS + - 1 e 式（3 -1- 8）如果： 满秩，且 e = 0 ，则 LS 为一致性无偏估计，如 e 0 ，则 LS + 为有偏估计。分析 e= E k e(k) k = -y(k-1),-y(k-n),u(k-1),u(k-n) T, 即 e不仅依赖于e(k)，也依赖于e(k-1)、e(k-2).，如果 e(k) 是白噪声，e(k) 与e(k-1)、e(k-2).不相关，则 e = 0 ， 如果 e(k) 是有色噪声, 则 e 0， + 问题：如何在有色噪声干扰条件下得到无偏估计？通过构造一个与 N相似，但是能
4、满足以下关系的矩阵Z N2n , 使得当N 时能满足满秩 式（3-1- 9）式（3 -1- 10）则： IV =（Z N T N） - 1 Z N T y N (无偏估计)其中： IVI V 估计；Z NI V 矩阵。01zeTN1 TNzZ1（2003 年教案） 辨识与自适应 第三章 4对应的关键在于如何具体构造出辅助变量矩阵Z N 。人们已经得出几种不同的方案，要严格证明它们满足上述条件是很难的，但是在实际应用中已经取得很好的结果。方案之一是用 LS 得出的模型的计算输出量代替中的实测输出组成 z 和Z ,即：z kT= -y(k-1),-y(k-n),u(k-1),u(k-n) 式（3-
5、1-11 ）其中：y(k)= z kT LS 式（3-1-12）对比： kT= -y(k-1),-y(k-n),u(k-1),u(k-n)二、I V 法的计算步骤a) 先用 LS法估计出 LSb) 由 LS 用式（2-3-12 ）计算出 y(k) k=1,2,N 组成矩阵 Z Nc) 计算 IV =（Z N T N） - 1 Z N T y N 式（3-1-12 ）三、递推辅助变量算法（R I V 法）R I V 法是在 I V 法思路启发下得出的，它并不与I V法等价。R I V 法由以下几个递推公式组成： N+1 = N+ K N+1(y(N+1) N+1T N) 式（3-1-13）（20
6、03 年教案） 辨识与自适应 第三章 5式（3-1-14）式（3-1-15）N+1T = - y(N)， - y(N- n+1) , u(N) ， , u(N- n+1) zN+1T = - y(N)， - y(N- n+1) , u(N) ， , u(N- n+1) 式（3-1-16）y(N) = zNT N 式（3-1-17）递推初值： 0 = 0 ；P 0= 106 I ;y(0) = y(-1) = = y(1-n) = 0四、仿真结果 I V 法与LS 法的比较参数 a1 a2 b1 b1真值 1.5 - 0.7 1.0 0.5LS 法 1.320 - 0.531 1.021 0.6
7、90I V 法 1.496 -0.697 1.017 0.504R I V法与RLS法的比较111 NTNNzPzP111TNNzK（2003 年教案） 辨识与自适应 第三章 6参数真值 a1=0.6 A2= - 0.7 b1=1.0 B2=0.5n=500 (0.6682)0.6598 (- 0.7144)- 0.7379 (1.0387)1.0298 (0.3679)0.4423n=1000 (0.6652)0.6237 (- 0.7214)- 0.7247 (1.0594)0.9855 (0.3691)0.5049n=1500 (0.6700)0.6156 (- 0.7291)- 0.7
8、129 (1.0556)1.0176 (0.3706)0.4953n=2000 (0.6840)0.6013 (- 0.7297)- 0.7111 (1.0568)1.0067 (0.2762)0.5122注：( ) 中为RLS法的估计结果五、IV 法评价a) I V 法计算量比线性最小二乘估计增加不多，而在相关噪声情况下估计精度有明显改善。b) I V 法与R I V 法思路相似，但不等价。c) R I V 法计算量与 RLS 法接近，在相关噪声情况下估计精度优于RLS 法，为了提高 R I V 法的可靠性，在递推的前几十步最好用RLS 法过渡。32 广义最小二乘估计方法（ GLS 法）Ge
9、neralized Least Squares Esimate一、算法思路考虑模型与最小二乘方法不同，由以下两个方程联立组成：A (z-1) y (k) = B (z-1) u (k) + (k) 式（3-2-1 ）和 C (z-1) (k) = (k) 式（3-2-2） (k) 为有色噪声， (k) 为白噪声。其中：A (z-1) = 1 + a1 z -1 + a n z -n（2003 年教案） 辨识与自适应 第三章 7B (z-1) = b1 z -1 + b n z -nC (z-1) = 1 + c1 z -1 + c r z -r由式（2-4-1）和式（2-4-2）有：A (z
10、-1) C (z -1) y (k) = B (z -1) C (z -1) u (k) + (k) 式（3-2-3）对 y (k) 和 u (k) 做如下变换：y (k) = C (z -1) y (k) 式（3-2-4）u (k) = C (z -1) u (k) 式（3-2-5）则模型成为：A (z -1) y (k) = B (z -1) u (k) + (k) 式（3-2-6 ）讨论：a) 如果C (z-1) 的参数c = c1 ,., c r T已知，可由 y (k) 和 u (k) 用式（3-2-4）和式（3-2-5）计算出 y (k) 和 u (k) ，由式（2-4-6）用 L
11、S 法求出 = a1,a n,b1,bn T ,可是 c 未知；b) 如果 已知，可由式（3-2-1）计算出 (k) ，式（3-2-2）为u (k) 0 的自回归 AR 模型，也可利用LS 法估计出 c , 但是 未知。GLS 法的思路是用迭代法逐次逼近地估计出 和 c（松弛算法） 。二、GLS 算法（1）先假设 C 0 (z-1) = 1 ，并置初值 0 = 0 ；（2）y (k) = C (z-1) y (k) 和 u (k) = C (z-1) u (k) 计算出 y (k) 和 u (k) ；（3）由 y (k) 和 u (k) 用 LS 法估计出 1 ；（2003 年教案） 辨识与自
12、适应 第三章 8 1= T - 1 T y 式（3-2-7 ）其中：（4）计算出 e (k) e (k) = y (k) - kT + 1 ; k=1,N 式（3-2-8）k = -y(k-1)-y(k-n),u(k-1)u(k-n) T(5) 对 C (z-1) e (k) = (k) 用LS 法估计：c1 = c1 . cr Tc1 = E T E -1 E T e 式（3-2-9）其中：（6）用下式校验是否继续迭代，若需要继续则转至步骤（2 ）再次估计。 足够小的正数 式（3-2-10）TNnNuunyy)(,.1 )(.)1()(. 0)0( TrNeN)(,.1).(.)0()(ma
13、x10ii（2003 年教案） 辨识与自适应 第三章 9三、GLS 计算框图开始置初值C(z -1)=1和 0 = 0 0 N计算 y (k) 和 u (k) 式（3-2-7）估计 N+1 N+1N由式(3-2-8)计算 e (k) 由式（3-2-9）估计c N+1否终点判断是输出模型（2003 年教案） 辨识与自适应 第三章 10四、仿真结果二阶过程A (z-1) y (k) = B (z-1) u (k) + (k) C (z-1) (k) = (k) 其中：A (z-1) = 1 + 1.5 z- 1 0.7 z- 2B (z-1) = 1.0 z- 1 + 0.5 z- 2C (z-1
14、) = 1 + 1.5 z- 1 0.8 z- r真值 1.5 -0.7 1.0 0.5 1.5 -0.8LS法 1.733 -0.963 1.329 0.137 GLS法六次迭代1.502 -0.730 1.038 0.471 1.488 -0.767五、GLS 算法的评价a) 有色干扰下估计精度较高b) 迭代收敛较快，但是收敛性未得到证明c) 能同时得到过程参数和噪声参数的估计d) 计算量大，费机时33 递推增广最小二乘法（ RELS 亦称增广矩阵法）一、算法思路考虑 CARMA 模型A (z-1) y (k) = B (z-1) u (k) + C (z-1) (k) 式（3- 3-1）
15、（2003 年教案） 辨识与自适应 第三章 11其中：A (z-1) = 1 + a1 z- 1 + a n z- nB (z-1) = b1 z- 1 + b n z- nC (z-1) = 1 + c1 z- 1 + c r z- r还可表示为 ：式（3- 3 -2）其中：以上两向量均由 2n 维扩展为 (2n+r) 维。在式（3- 3 -2）中若 ( k-1)、 、 、 ( k-r) 是已知量，则可直接用 LS 法估计出 ，但是 ( k) 是不可测的未知量。一个简单可行的方法是用计算的 ( k-1) 代替 ( k-1) ,. ( k-r) , ( k) 由下式递推得出： (k) = y(
16、k) - k k-1 (计算残差) 式（3- 3 3）其中：k = -y(k-1),-y(k-n) , u(k-1),u(k-n) , (k-1), (k-r) T)()(kkyTk 1)2(11,.,.,. )(),.1(,.,),.( rnTn TKcba rkkuky （2003 年教案） 辨识与自适应 第三章 12二、递推算式 N+1 = N+ K N+1(y(N+1) N+1T N) 式（3-3-4 ）式（3-3-5）式（3-3-6） (N+1) = y(N) - N+1 N 式（3 37）N = -y(N-1),-y(N-n),u(N-1),u(N-n),(N-1), (N-r)
17、T式（3-3-8 ）递推初值： 0= 0 ; P 0= 10 6 I; (0) = (-1) = = (1-r) = 0三、仿真结果2 阶过程， ( k) 为零均值 1方差的伪白噪声序列， u ( k) 为周期15、幅度5的方波，遗忘因子为 0.99 。真 值 a1=0.5 a2=- 0.7 b1=1.0 b2=0.5 c1=0.9 c2=0.7N a1 估计 a2 估计 b1 估计 b2 估计 c1 估计 c2 估计200 0.450 -0.658 1.044 0.518 0.560 0.165400 0.479 -0.710 1.005 0.559 0.748 0.421600 0.503
18、 -0.704 0.991 0.535 0.924 0.697800 0.502 -0.703 0.968 0.535 0.896 0.635111 NTNNPP111TNN（2003 年教案） 辨识与自适应 第三章 131000 0.514 -0.695 0.987 0.473 0.846 0.63034 相关最小二乘两步法 （Cor - LS 法）一、思路和算法离散随机序列u(k) 和 y(k) (平稳遍历) 有：式（3-4-1）式（3-4-2）考虑过程：y(k)+a1y(k-1) +any(k-n) = b1u(k-1)+bnu(k-n)+ (k)式（3-4-3）设 u(k) 与 (k)
19、 不相关，即 Ru 0 , 将式（3-4-3）左右两边同乘以 u(k-) u(k-)y(k)+a1u(k-)y(k-1) +an u(k-)y(k-n) = b1 u(k-)u(k-1)+ bn u(k-)u(k-n)+ u(k-) (k)令k=0 , N-1，共得N个等式，将N 个等式相加并除以 N , 得出 ：Ruy()+a1 Ruy(-1) +an Ruy(-n) = b1 Ruu(-1)+ + bn Ruu(-n)+ h()式（3-4-4）注：R u () 0 ，h() 项考虑由于利用有限 N 值计算R uy 和R uu 的误差。1010)()( )()(Nknuyku kyuRLim
20、i（2003 年教案） 辨识与自适应 第三章 14如果从样本数据计算出 = - (n-1),0,1,L,共L组相关函数，由式（3-4-4）得出 L个由R uy（*）和R uu（*）组成的方程组：L 1 L 2n 2n 1 L 1式（3-4-5 ）上式可表达成如下矩阵形式：式（3-4-6） 上式与最小二乘的 yN = N N + N 相似，故用 LS法可得出参数估计 为： = (RLT RL) - 1 RLT+gL 式（3-4-7）为保证 (RLT RL) 满秩，要求 L 2n做为特例：如 u(k) 为白噪声序列，R uy() = g()，上式可简化为 )()1(*)()1()()( )2()(
21、2()1( 100)()2( 1hhbanRnRnnRuuuyuy uuuyuyuyuy lll hRg )()2(1*)(.)1(.0)(2.1.)()( 1)( LhhbanLgLgOIg nnnLn（2003 年教案） 辨识与自适应 第三章 15式（3-4-8）二、计算步骤1 由u(k) 和 y(k)用式（3-4-1）和式（3-4-2）计算出 Ruu() 和Ruy() , = - (n-1),0,1,L 。可由维纳何甫方程可得出g()2由式（3-4-7 ）估计出 三、两步法的特点1、只要求 (k) 是与 u(k) 不相关的零均值平稳噪声，则 (k) 不影响辨识结果。并不要求 (k) 为白
22、噪声。2、该方法能同时获得非参数模型 g(t) 和参数模型 3、计算量不太大，乐于为工程界采用。四、应用举例例1热交换器的在线辨识 ：效果比RLS 法好例2电厂锅炉动态特性在线辨识（南京工学院）附：最小二乘及其改进算法的比较1） LS、IV 、GLS、ELS 和 Cor-LS 五种算法同属LS法体系。2） 五种算法的比较算法 模型获得一致性无偏估计对噪声的要求计算量 特点luuy Rg1)()(（2003 年教案） 辨识与自适应 第三章 16LS Ay(k)=Bu(k)+(k) 白噪声 最小算法可靠应用广泛IV Ay(k)=Bu(k)+e(k) 有色噪声中等RIV小没有辩识噪声模型GLSAy(
23、k)=Bu(k)+(k)R(k)= (k)有色噪声模型为 AR大 辩识了噪声模型ELS Ay(k)=Bu(k)+C(k)有色模型为 CARMA较小辩识噪声模型,只有递推算法CorLSAy(k)=Bu(k)+e(k) 有色噪声 中等 同时得到两类模型1 实际应用并不要求无偏估计，在线辩识中以 RLS 法应用最多，离线辩识往往用多种方法，通过比较择优。2 LS 法的局限性3 其它算法：随机逼近法（STA）,极大似然法（LM）35 应用举例例1. 弹射加速度作用下人体动态响应的模型辨识和几种方法的比较（自动化学报）（2003 年教案） 辨识与自适应 第三章 17例2. 中医脉象辨识 “信息与控制”8
24、2.5思路：用系统辨识方法解决中医脉象的图象识别问题y(k)脉象数据； u(k)虚拟的输入函数（2003 年教案） 辨识与自适应 第三章 18将脉象视为在虚拟输入函数作用下的二阶系统的响应典型缓脉 滑脉 弦脉 芤脉（2003 年教案） 辨识与自适应 第三章 19（2003 年教案） 辨识与自适应 第三章 21（2003 年教案） 辨识与自适应 第三章 22例 3、基于系统辨识方法的电动车动力蓄电池放电性能研究（ 2001，EVS-18, Berlin）电动车动力蓄电池充放电过程的工作机理复杂，工作性能受到诸多因素的影响，它们之间的定量关系至今还不很清楚，如若使用不当，特别是过度充电和过度放电，
25、将严重影响蓄电池的使用寿命。国内外不少学者，从不同技术路线研究蓄电池的放电过程的数学模型的研究。目前能够被电动车采用的动力蓄电池主要有铅酸蓄电池、胶体铅酸蓄电池、镍氢蓄电池和锂离子蓄电池四种，它们的放电性能有明显区别。可以通过实验数据，利用系统辨识方法得出四种蓄电池充放电过程电流与端电压之间的动态响应的数学模型，进行比较研究，为蓄电池管理系统的设计提供依据。蓄电池放电状态的数学模型等值电路为（2003 年教案） 辨识与自适应 第三章 23蓄电池放电状态的等值电路蓄电池空载端电压为 U0 ，当蓄电池对负载 Ra 放电时，由于电池内阻抗Z 0 的压降，端电压由U 0下降至 Un ，两者之差为UU
26、= U0 Ua 蓄电池的内阻抗 Z 0 并不是纯电阻性质的，而是有一定容性的，而且与电池工作状态，特别是与剩余电量有密切关系，这正是我们要通过大量实验，进行研究的内容。经过实验比较，我们选择采样间隔 T0 = 1/4 秒，对 U n 和 I a 采样，通过对蓄电池在不同放电条件下的实验采样，获得了一批实验数据样本。利用系统辨识的时域建模方法，借助于线性最小二乘估计获得模型参数的估计。蓄电池放电状态的数学模型结构图为蓄电池放电状态内阻抗模型的结构图（2003 年教案） 辨识与自适应 第三章 24采用以下双输入单输出的自回归（CAR）线性的、时间离散差分方程模型y(k) = a y(k-1) +
27、b0 u(k-1)+b1 u(k-2) + c (k-k0) +(k) 其中：采样序列 u(k) 为可测的输入，当蓄电池负载回路开关 K 接通时 u(k) = Uo，而开关 K 断开时 u(k) = 0； y(k) 为可测的输出量，即蓄电池内阻抗压降的采样值；而 (k-k0) 可被视为另一个可测输入，参数 c 用于辨识模型中的衰减系数 ，k 0 为功率开关K 接通时刻的采样序号， (k) 为不可测的随机干扰量（测量干扰和模型误差） 。下图为用 RLS 方法的参数估计的结果曲线之一例。下图为镍氢蓄电池在不同放电安时后的动态压降的比较，其中红线为模型输出，黑线为实测输出，可见两者拟合得很好，表明用
28、系统辨识和参数估计的方法得到的模型是可信的。（2003 年教案） 辨识与自适应 第三章 25蓄电池放电模型输出与实测值的比较和收放电安时的影响可见，经过不同放电安时后的动态压降波形有非常明显的区别，也就是模型参数的大小受到放电安时数的强大影响。以下几图给出几种蓄电池放电模型时间常数 T 与放电 Ah 的关系的比较，可见区别之明显。镍氢蓄电池时间常数 T 与放电 Ah 的关系00.511.522.533.5秒 时间常数 T2 6 9.25 放 电 Ah0.5C放 电 （2003 年教案） 辨识与自适应 第三章 26锂离子蓄电池时间常数 T 与放电 Ah 的关系铅酸蓄电池时间常数 T 与放电 Ah
29、 的关系36 极大似然法（ ML法）一、极大似然估计与预报误差估计极大似然估计是与最小二乘估计相平行的另一种参数估计方法，由 Fisher (1912 年) 提出。该方法估计性质较好，既使噪声较强时也能得到较好的估计结果，但是计算复杂，而且要求已知输出量的条件概率密度为先验知识。预报误差方法是极大似然法的一种推广，某种意义下两者是等价的，而且通过预报误差方法可以看出 ML 法00.511.522.533.5秒 时间常数T1 3 5 7 8 放 电 Ah0.5C放 电 00.511.522.5秒 时间常数T0 2 4 6 8 放 电 Ah0.5C放 电 （2003 年教案） 辨识与自适应 第三章
30、 27与 LS 法的联系。二、极大似然估计构造一个似然函数 L(y,), 它是观测数据与未知参数的函数。Y 是随机变量的观测值， 为未知参数。似然函数是总体 Y 的条件概率密度L(y,)= p(y) 式（3-6-1 ） 。使得 L(y,) 达到最大的参数估计 LM 称为极大似然估计，即：p ( y LM )= max p( y ) 式（3-6-2 ） LM 意味着使模型输出的概率分布，最大可能地接近实测输出的概率分布，即是使出现观测值 y 出现的可能性最大的估计。由于 ln( )是单调增函数，还可引出对数似然函数 ln L( y, )。极大似然估计 LM 必须满足以下似然方程：式（3-6-3）
31、式（3-6-4）通过求解上述非线性代数方程，得出 LM 是很难的。以下仅讨论正态分布的情况。1 独立观测（静态）观测值 y(1),y(2)y(N);它们的联合概率密度可以简单地表示为每个0),(ln),(MLLyL（2003 年教案） 辨识与自适应 第三章 28观测的概率密度的乘积（Bayes 公式） ，即似然函数Ly(1),y(N) ; =p(y(1); )p(y(1); )p(y(1); )式（3-6-5）正态分布：则其中 (k) = y(k)- 将上式两边同取对数并反号，得出：式（3-6-6）L 的极大化即为上式（3-6-6）的极小化，若 为已知，即相当于使下式的极小化：即等价于最小二乘
32、估计。Niiy1);(NkNNyL12)(exp)2();(,.1NkN122lnl)(lnNkJ12)(2)(1);(yp（2003 年教案） 辨识与自适应 第三章 292 序贯观测（动态）这时的似然函数（联合概率密度函数）不能如同式（3-5-5 ）那样简单。可借助于一种更一般化的方法“预报误差估计”来解。三、预报误差模型和参数的预报误差估计引入一个更为一般的预报方程，描述受随机干扰的系统： y(k)= fYk-1 ,Uk ,k; + (k) 式（3-6-7）其中：Yk-1= y(k-1),y(k-2), 输出历史数据的集合Uk-1= u(k-1),u(k-2), 输入历史数据的集合F 是确
33、定的，y(k)和 (k)为随机量。若 E (k) Yk-1 ,Uk = 0 k称 (k) 为新息序列。对于 = 时的预报误差为 w (k, )w (k, ) = y(k)- fYk-1 ,Uk ,k; = y(k)-y(kk-1) 式（3-6-8）对于多维过程 y(k) Rm , w 为 m 维向量。N 次观测的预报误差样本协方差阵为 ：式（3-6-9）D() 为 mm 方阵。选择与 D() 有关的正值的标量函数做估计准则，以下为其中常用的两种：),(),(1)(kNDT（2003 年教案） 辨识与自适应 第三章 30J1() = tr R D() 式（3-6-10 ）J2() = log d
34、et D() 式（6-6-11 ）上式中的 R 为正定加权矩阵使 J1 或 J2 达到极小的 称为预报误差估计。可以证明：若 R = I , 则预报误差估计等同于 LS 估计；若 w 为正态分布，预报误差估计即 ML 估计。所以，LS 和 ML 估计在一定条件下都可视为预报误差估计的特例。用最优化算法，通过 J1() 或者 J2()的极小化求解 。四、极大似然估计的松弛算法介绍一种估计有色噪声干扰下多变量模型参数的算法，可视为GLS 法在多变量系统的推广。考虑以下多输入多输出线性系统A(z-1) y(k) = B(z-1) u(k) + (k) 式（3-6-12 ）其中：y Rm Rm u R
35、pA(z-1)= I + a1 z-1 Im + an z-n Im 式（3-6-13）B(z-1)= Im+ B1 z-1Im + Bn z-nIm 式（3-6-14）ai 、 bi 为标量, I 为 mm 方阵,B I 为 mp 方阵。(k) 可表示为白噪声 (k) 的自回归模型：F (z-1) (k) = (k) 式（3-6-15 ）其中 (k) 为 m 维正态分布白噪声序列，具有零均值、同分布，协方差阵为 mm 。（2003 年教案） 辨识与自适应 第三章 31F(z-1)= Im + F1 z-1 Im + Fr z-r Im 式（3-6-16）其中：F i ( i = 1,r ) 为 mm 常参数阵。例：双输入双输出的 2 阶系统和采用以下模型拟合式（3-6-12）过程的观测数据A(z-1)y(k) - B(z-1)u(k) = (k) 式（3-6-17）预报误差w (k, ) = F(z-1) (k) (m1) 式（3-6-18）对数似然函数 L() = log P YNUN , 的极大化等价于以下目标的极小化：L() = - log P YNUN , =式(3-6-19)其中 是 (k) 的协方差阵 )(2(.)1(.)2()1() 1122121 kubkubkyaky )().( 21222 kkFkNkTkwmN11),(),(2detlog2l
