参数估计的改进算法.doc

（2003年教案） 辨识与自适应 第三章 32第三章 参数估计的改进算法31最小二乘方法的改进辅助变量法 （I V法 Instrument Variable ） 辅助变量法是在模型误差为相关噪声的情况下，通过引进辅助变量矩阵，对线性最小二乘估计的一种改进。一、参数的辅助变量估计考虑模型 式（3- 1- 1）其中：当进行了 k = 1-n,2-n,.,0,1,2,N 共计（N+n）次采样，得到N个方程：用矩阵表示成 y N = FN q + e N 式（3 -1- 2）其中 e N = e(1),.,e(N)T最小二乘估计为：q LS =（F N T F N）-1 F N T y N设（F N T F N）满秩，将过程模型式（2-3-2）带入上式，得出：q LS =（F N T F N）-1 F N T（FN q + e N） 式（3-1- 3）有：和令： Yjj= E jk jkT2n2n 式（3 -1- 4） 和 Yje= E jk e(k) 2n1 式（3 -1- 5）可以证明：当N 时 式（3-1- 6）式（3-1- 7）依据 Frechek 定理： 若随机序列 xK x ， 则 f(xK) f( x ) 故有：q LS q + Yjj -1 Yje 式（3 -1- 8）如果：Yjj 满秩，且Yje = 0 ，则 q LS q 为一致性无偏估计，如Yje 0 ，则 q LS q +D q为有偏估计。分析Yje= E jk e(k) jk = -y(k-1),-y(k-n),u(k-1),u(k-n) T, 即 Yje不仅依赖于e(k)，也依赖于e(k-1)、e(k-2).，如果 e(k) 是白噪声，e(k) 与e(k-1)、e(k-2).不相关，则 Yje = 0 ， q q如果 e(k) 是有色噪声, 则 Yje 0， q q+D q问题：如何在有色噪声干扰条件下得到无偏估计？通过构造一个与F N相似，但是能满足以下关系的矩阵ZN2n , 使得当N 时能满足 满秩 式（3-1- 9） 式（3 -1- 10）则：q IV =（Z N T F N）-1 Z N T y N q (无偏估计)其中：q IVI V 估计；Z NI V 矩阵。对应的关键在于如何具体构造出辅助变量矩阵ZN 。人们已经得出几种不同的方案，要严格证明它们满足上述条件是很难的，但是在实际应用中已经取得很好的结果。方案之一是用q LS 得出的模型的计算输出量代替F中的实测输出组成 z 和Z ,即：z kT= -y(k-1),.,-y(k-n),u(k-1),u(k-n) 式（3-1-11） 其中：y(k)= z kT q LS 式（3-1-12）对比：j kT= -y(k-1),.,-y(k-n),u(k-1),u(k-n)二、I V 法的计算步骤a) 先用LS法估计出 q LSb) 由 q LS 用式（2-3-12）计算出 y(k) k=1,2,N 组成矩阵 Z Nc) 计算 q IV =（Z N T F N）-1 Z N T y N 式（3-1-12）三、递推辅助变量算法（R I V 法）R I V 法是在 I V 法思路启发下得出的，它并不与I V法等价。R I V 法由以下几个递推公式组成： q N+1 = q N+ K N+1(y(N+1) jN+1T q N) 式（3-1-13） 式（3-1-14） 式（3-1-15） jN+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 qN 式（3-1-17）递推初值：q 0 = 0 ；P 0= 106 I ; y(0) = y(-1) = = y(1-n) = 0四、仿真结果 I V 法与LS法的比较参数a1a2b1b1真值1.5- 0.71.00.5LS 法1.320- 0.5311.0210.690I V 法1.496-0.6971.0170.504 R I V法与RLS法的比较参数真值a1=0.6A2= - 0.7b1=1.0B2=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.7129(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 法） Generalized Least Squares Esimate一、算法思路考虑模型与最小二乘方法不同，由以下两个方程联立组成： A (z-1) y (k) = B (z-1) u (k) + x (k) 式（3-2-1）和 C (z-1) x (k) = e (k) 式（3-2-2） x (k) 为有色噪声， e (k) 为白噪声。其中：A (z-1) = 1 + a1 z -1 + a n z -n B (z-1) = b1 z -1 + b n z -n C (z-1) = 1 + c1 z -1 + c r z -r由式（2-4-1）和式（2-4-2）有：A (z -1) C (z -1) y (k) = B (z -1) C (z -1) u (k) + e (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) + e (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）用 LS 法求出 q = a1,a n,b1,bn T ,可是 c 未知； b) 如果 q 已知，可由式（3-2-1）计算出 x (k) ，式（3-2-2）为u (k) 0 的自回归 AR 模型，也可利用LS 法估计出 c , 但是 q 未知。GLS 法的思路是用迭代法逐次逼近地估计出 q 和 c（松弛算法）。二、GLS 算法（1）先假设C 0 (z-1) = 1 ，并置初值 q 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 法估计出 q 1 ； q 1= f*T f* - 1 f*T y * 式（3-2-7）其中：（4）计算出 e (k) e (k) = y (k) - jkT + q1 ; k=1,N 式（3-2-8） jk = -y(k-1)-y(k-n),u(k-1)u(k-n) T(5) 对 C (z-1) e (k) = e (k) 用LS法估计：c1 = c1 . cr T c1 = E T E -1 E T e 式（3-2-9）其中：（6）用下式校验是否继续迭代，若需要继续则转至步骤（2）再次估计。 e 足够小的正数 式（3-2-10）三、GLS 计算框图 开始 置初值C(z-1)=1和 q 0 = 0 0 N 计算 y * (k) 和 u * (k) 式（3-2-7）估计q N+1 N+1N 由式(3-2-8)计算 e (k) 由式（3-2-9）估计c N+1 否 终点判断 是 输出模型四、仿真结果二阶过程A (z-1) y (k) = B (z-1) u (k) + x (k) C (z-1) x (k) = e (k) 其中：A (z-1) = 1 + 1.5 z- 1 0.7 z- 2 B (z-1) = 1.0 z- 1 + 0.5 z- 2 C (z-1) = 1 + 1.5 z- 1 0.8 z- r真值1.5-0.71.00.51.5-0.8LS法1.733-0.9631.3290.137GLS法六次迭代1.502-0.7301.0380.4711.488-0.767五、GLS 算法的评价a) 有色干扰下估计精度较高b) 迭代收敛较快，但是收敛性未得到证明c) 能同时得到过程参数和噪声参数的估计d) 计算量大，费机时33递推增广最小二乘法（RELS 亦称增广矩阵法）一、算法思路考虑 CARMA 模型A (z-1) y (k) = B (z-1) u (k) + C (z-1) e (k) 式（3- 3-1）其中：A (z-1) = 1 + a1 z- 1 + a n z- n B (z-1) = b1 z- 1 + b n z- n C (z-1) = 1 + c1 z- 1 + c r z- r还可表示为 ：式（3- 3 -2）其中： 以上两向量均由 2n 维扩展为 (2n+r) 维。在式（3- 3 -2）中若e ( k-1)、e ( k-r) 是已知量，则可直接用 LS 法估计出 q ，但是 e ( k) 是不可测的未知量。一个简单可行的方法是用计算的e ( k-1) 代替e ( k-1) ,. e ( k-r) , e ( k) 由下式递推得出： e (k) = y(k) - jk q k-1 (计算残差) 式（3- 3 3）其中： jk = -y(k-1),.,-y(k-n) , u(k-1),.,u(k-n) , e (k-1),., e (k-r) T二、递推算式 q N+1 = q N+ K N+1(y(N+1) jN+1T q N) 式（3-3-4） 式（3-3-5） 式（3-3-6） e (N+1) = y(N) - jN+1 q N 式（337） jN = -y(N-1),.,-y(N-n),u(N-1),u(N-n),e(N-1),.,e (N-r) T 式（3-3-8）递推初值： q 0= 0 ; P 0= 10 6 I; e (0) = e (-1) = = e (1-r) = 0三、仿真结果 2 阶过程， e ( k) 为零均值1方差的伪白噪声序列， u ( k) 为周期15、幅度5的方波，遗忘因子为0.99 。真 值a1=0.5a2=- 0.7b1=1.0b2=0.5c1=0.9c2=0.7Na1 估计a2 估计b1 估计b2 估计c1 估计c2 估计2000.450-0.6581.0440.5180.5600.1654000.479-0.7101.0050.5590.7480.4216000.503-0.7040.9910.5350.9240.6978000.502-0.7030.9680.5350.8960.63510000.514-0.6950.9870.4730.8460.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)+ e(k) 式（3-4-3） 设 u(k) 与 e(k) 不相关，即 Ru e 0 , 将式（3-4-3）左右两边同乘以 u(k-t) u(k-t)y(k)+a1u(k-t)y(k-1) +an u(k-t)y(k-n) = b1 u(k-t)u(k-1)+ bn u(k-t)u(k-n)+ u(k-t) e(k)令k=0 , N-1，共得N个等式，将N 个等式相加并除以 N , 得出 ：Ruy(t)+a1 Ruy(t-1) +an Ruy(t-n) = b1 Ruu(t-1)+ + bn Ruu(t-n)+ h(t) 式（3-4-4） 注：Ru e (t) 0 ，h(t) 项考虑由于利用有限N 值计算Ruy 和Ruu 的误差。如果从样本数据计算出t = - (n-1),.,0,1,L,共L组相关函数，由式（3-4-4）得出 L个由Ruy（*）和Ruu（*）组成的方程组： L 1 L 2n 2n 1 L 1 式（3-4-5）上式可表达成如下矩阵形式： 式（3-4-6） 上式与最小二乘的 yN = FN qN + eN 相似，故用LS法可得出参数估计 q 为：q = (RLT RL) -1 RLT+gL 式（3-4-7）为保证 (RLT RL) 满秩，要求 L 2n做为特例：如 u(k) 为白噪声序列，Ruy(t) = g(t)，上式可简化为 式（3-4-8）二、计算步骤1 由u(k) 和 y(k)用式（3-4-1）和式（3-4-2）计算出 Ruu(t) 和Ruy(t) , t = - (n-1),.,0,1,L 。可由维纳何甫方程可得出g(t)2由式（3-4-7）估计出 q三、两步法的特点1、只要求 e(k) 是与 u(k) 不相关的零均值平稳噪声，则 e(k) 不影响辨识结果。并不要求 e(k) 为白噪声。2、该方法能同时获得非参数模型 g(t) 和参数模型 q3、计算量不太大，乐于为工程界采用。四、应用举例例1热交换器的在线辨识 ：效果比RLS法好例2电厂锅炉动态特性在线辨识（南京工学院）附：最小二乘及其改进算法的比较1） LS、IV、GLS、ELS 和 Cor-LS 五种算法同属LS法体系。2） 五种算法的比较算法模型获得一致性无偏估计对噪声的要求计算量特点LSAy(k)=Bu(k)+e(k)白噪声最小算法可靠应用广泛IVAy(k)=Bu(k)+e(k)有色噪声中等RIV小没有辩识噪声模型GLSAy(k)=Bu(k)+x(k)Rx(k)= e(k) 有色噪声模型为 AR大辩识了噪声模型ELSAy(k)=Bu(k)+Ce(k) 有色模型为 CARMA较小辩识噪声模型,只有递推算法Cor|LSAy(k)=Bu(k)+e(k)有色噪声中等同时得到两类模型1 实际应用并不要求无偏估计，在线辩识中以 RLS 法应用最多，离线辩识往往用多种方法，通过比较择优。2 LS 法的局限性3 其它算法：随机逼近法（STA）,极大似然法（LM）35 应用举例例1. 弹射加速度作用下人体动态响应的模型辨识和几种方法的比较（自动化学报）例2. 中医脉象辨识 “信息与控制”82.5思路：用系统辨识方法解决中医脉象的图象识别问题y(k)脉象数据；u(k)虚拟的输入函数将脉象视为在虚拟输入函数作用下的二阶系统的响应 典型缓脉 滑脉 弦脉 芤脉例3、基于系统辨识方法的电动车动力蓄电池放电性能研究（ 2001，EVS-18, Berlin）电动车动力蓄电池充放电过程的工作机理复杂，工作性能受到诸多因素的影响，它们之间的定量关系至今还不很清楚，如若使用不当，特别是过度充电和过度放电，将严重影响蓄电池的使用寿命。国内外不少学者，从不同技术路线研究蓄电池的放电过程的数学模型的研究。目前能够被电动车采用的动力蓄电池主要有铅酸蓄电池、胶体铅酸蓄电池、镍氢蓄电池和锂离子蓄电池四种，它们的放电性能有明显区别。可以通过实验数据，利用系统辨识方法得出四种蓄电池充放电过程电流与端电压之间的动态响应的数学模型，进行比较研究，为蓄电池管理系统的设计提供依据。蓄电池放电状态的数学模型等值电路为蓄电池放电状态的等值电路蓄电池空载端电压为 U0 ，当蓄电池对负载 Ra 放电时，由于电池内阻抗Z0 的压降，端电压由U0下降至 Un ，两者之差为UU = U0 Ua 蓄电池的内阻抗 Z0 并不是纯电阻性质的，而是有一定容性的，而且与电池工作状态，特别是与剩余电量有密切关系，这正是我们要通过大量实验，进行研究的内容。经过实验比较，我们选择采样间隔 T0 = 1/4 秒，对U n 和 I a 采样，通过对蓄电池在不同放电条件下的实验采样，获得了一批实验数据样本。利用系统辨识的时域建模方法，借助于线性最小二乘估计获得模型参数的估计。蓄电池放电状态的数学模型结构图为蓄电池放电状态内阻抗模型的结构图采用以下双输入单输出的自回归（CAR）线性的、时间离散差分方程模型 y(k) = a y(k-1) + 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用于辨识模型中的衰减系数，k0 为功率开关K 接通时刻的采样序号， e(k) 为不可测的随机干扰量（测量干扰和模型误差）。下图为用RLS方法的参数估计的结果曲线之一例。下图为镍氢蓄电池在不同放电安时后的动态压降的比较，其中红线为模型输出，黑线为实测输出，可见两者拟合得很好，表明用系统辨识和参数估计的方法得到的模型是可信的。 蓄电池放电模型输出与实测值的比较和收放电安时的影响可见，经过不同放电安时后的动态压降波形有非常明显的区别，也就是模型参数的大小受到放电安时数的强大影响。以下几图给出几种蓄电池放电模型时间常数T 与放电Ah 的关系的比较，可见区别之明显。 镍氢蓄电池时间常数T 与放电Ah 的关系锂离子蓄电池时间常数T 与放电Ah 的关系 铅酸蓄电池时间常数T 与放电Ah 的关系 36 极大似然法（ML法）一、极大似然估计与预报误差估计 极大似然估计是与最小二乘估计相平行的另一种参数估计方法，由Fisher (1912年) 提出。该方法估计性质较好，既使噪声较强时也能得到较好的估计结果，但是计算复杂，而且要求已知输出量的条件概率密度为先验知识。预报误差方法是极大似然法的一种推广，某种意义下两者是等价的，而且通过预报误差方法可以看出ML法与LS法的联系。二、极大似然估计构造一个似然函数L(y,q), 它是观测数据与未知参数的函数。Y是随机变量的观测值，q 为未知参数。似然函数是总体Y的条件概率密度L(y,q)= p(yq) 式（3-6-1）。使得L(y,q) 达到最大的参数估计qLM 称为极大似然估计，即： p ( yq LM )= maxq p( yq ) 式（3-6-2）q LM 意味着使模型输出的概率分布，最大可能地接近实测输出的概率分布，即是使出现观测值y出现的可能性最大的估计。由于ln( )是单调增函数，还可引出对数似然函数 ln L( y,q )。极大似然估计 q LM 必须满足以下似然方程： 式（3-6-3）式（3-6-4） 通过求解上述非线性代数方程，得出 q LM是很难的。以下仅讨论正态分布的情况。1 独立观测（静态）观测值y(1),y(2)y(N);它们的联合概率密度可以简单地表示为每个观测的概率密度的乘积（Bayes 公式），即似然函数Ly(1),.y(N) ; q=p(y(1); q)p(y(1); q)p(y(1); q) 式（3-6-5）正态分布：则 其中 e(k) = y(k)- m将上式两边同取对数并反号，得出：式（3-6-6） L的极大化即为上式（3-6-6）的极小化，若s为已知，即相当于使下式的极小化：即等价于最小二乘估计。2 序贯观测（动态）这时的似然函数（联合概率密度函数）不能如同式（3-5-5）那样简单。可借助于一种更一般化的方法“预报误差估计”来解。三、预报误差模型和参数的预报误差估计引入一个更为一般的预报方程，描述受随机干扰的系统： y(k)= fYk-1 ,Uk ,k; q+ e(k) 式（3-6-7）其中： Yk-1= y(k-1),y(k-2), 输出历史数据的集合 Uk-1= u(k-1),u(k-2), 输入历史数据的集合 F 是确定的，y(k)和 e(k)为随机量。 若 E e(k) Yk-1 ,Uk = 0 k称 e(k) 为新息序列。对于 q = q 时的预报误差为w (k, q)w (k, q) = y(k)- fYk-1 ,Uk ,k; q= y(k)-y(kk-1) 式（3-6-8） 对于多维过程y(k) Rm , w 为m维向量。N次观测的预报误差样本协方差阵为 ： 式（3-6-9）D(q) 为mm 方阵。 选择与D(q) 有关的正值的标量函数做估计准则，以下为其中常用的两种： J1(q) = tr R D(q) 式（3-6-10） J2(q) = log det D(q) 式（6-6-11）上式中的R 为正定加权矩阵使J1 或J2达到极小的q 称为预报误差估计。可以证明：若R = I , 则预报误差估计等同于LS估计；若w为正态分布，预报误差估计即ML估计。所以，LS和ML估计在一定条件下都可视为预报误差估计的特例。用最优化算法，通过J1(q) 或者J2(q)的极小化求解 q 。四、极大似然估计的松弛算法 介绍一种估计有色噪声干扰下多变量模型参数的算法，可视为GLS法在多变量系统的推广。考虑以下多输入多输出线性系统A(z-1) y(k) = B(z-1) u(k) + x(k) 式（3-6-12）其中：y Rm x Rm u RpA(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 方阵,BI为mp 方阵。x(k) 可表示为白噪声e (k) 的自回归模型： F (z-1) x(k) = e (k) 式（3-6-15）其中 e (k) 为m维正态分布白噪声序列，具有零均值、同分布，协方差阵为Smm 。F(z-1)= Im + F1 z-1 Im + Fr z-r Im 式（3-6-16）其中：Fi ( i =