建模与估计下.ppt

建模与估计讲义（二）主讲人：王欣,第三节 新息序列(Innovation sequence),定义:设(相关随机序列正交化)有二阶矩(数学期望、方差都）的随机序列，定义它的新息序列为：定义:基于k-1以前信息对y(k)的线性最小方差估计 则 规定:定理1：新息序列是正交序列（白噪声）证：射影具有无偏性 ij 不妨设ij根据射影性质,定理2:,证明:显然,注意,事实上:,重要意义:新息序列 与原序列y(k)含有相同的统计信息,新息实质上是摄影误差,定理3:递推射影公式 证明:,考虑线性离散定常随机系统 其中状态x(t),观测y(t),w(t)输入噪声.观测噪声v(t)假设1:w(t)和v(t)是零均值,方差阵各为Q和R的独立白噪声 即:,第四节 kalman滤波,假设2:x(0)与w(t)和v(t)相互独立Kalman滤波问题是基于观测(y(1)y(t)求状态 x(j)的线性最小方差估计极小化即求 滤波 平滑 预报,分析:有递推射影公式 k(t+1)称为滤波增益对取射影 因为 事实上,新息:引入定义:,即得即,即注意:,简化:即:总结:考虑随机系统问题:基于观测 求,Kalman Filter初值 协方差阵,回忆:考虑随机系统问题:基于观测 求定理:Kalman Filter,第四节,初值:协方差阵例1:已知 及 求解:注1:Kalman滤波的和预报器的不同形式封闭形式:,滤波器 预报增益预报器,注2:Riccati方程即:初值:注3:时变系统Kalman方程组同样适用,有Kalman方程组 初值:注4:Kalman滤波与RLS关系(RLS是K的特例)考虑AR(N),可写成状态空间模型应用Kalman Filter定义：具有RLS公式,注5：稳态Kalman滤波和预报器（steady-state kalman filter and pred！ctor）Kalman方程组缺点：需要在线（on-line）实时（real time）计算p问题：P（t+1|t）?（t）定常系统 即 h p为常阵定理：设系统是完全可观、完全可控或 是稳定阵即：H rank H=n rankp p*n-1p=n H*2 H*n-1则任意 p(t+1|t)=它满足Riccati方程=-H*T(H H*H+R)*(-1)H*T+PQP*T,且有关系 k（t）=k p（t|t）=p k=H H H*T+R*-1=p*T=PQP*T p=-KH 稳态Kalman滤波和预报器为 x(t+1|t+1)=x(t|t)+ky(t+1)=-KH x(t+1|t)=x(t|t-1)+y(t)=-KH=k注释：可证 为稳定阵 因而 x(t|t),x(t+|t)是渐近稳定的。即可任意选取初值 x(0|0)，或（0|0）例：考虑一维系统 x(t+1)=x(t)+bw(t)y(t)=hx(t)+v(t)Q=*2=q R=*2=r求稳态kalman滤波器和预报器 x(t+1|t),x(t+1|t+1)=x(t|t)+ky(t+1)=(1-kh)x(t+1|t)=x(t|t)-1+y(t)=-kph=kRiccati方程=-h（hh+r）*-1h+bqb=*2-h*2(h*2+r)*-1+b*2q=*2(h*2+r)-*2h*2/h*2+r+b*2q(h*2+r)=*2 r+b*2q(h*2+r)h*2*2+r=*2 r+b*2q h*2+b*2qrh*2*2+(r-*2r-b*2h*2q)-b*2qr=0作一元二次方程，取0为GT为,第五讲 ARMA时间序列预报（ARMA Time series prediction)Forcasting 1、引言 预报问题：气象、水文、经济系统、控制,最优预报，稳态线性最小方差预报稳态：=，已知无穷的观测历史线性预报：y(t+k|t)L(y(t)y(t-1)是以前历史的线性组合最小方差：minJ=E(y(t+k)-y(t+k|t)*2 y(t+k|t),Box-Jenkins递推预报器.考虑平稳可逆的ARMA过程y(t).A(q*1)y(t)=c(q*1)e(t)其中e(t)是白噪声：E e(t)=0 Ee(t)e(s)=*2A(q*1)=1+q*1+q*C(q*1)=1+q*1+q*已知（y(t)y(t-1)）求y(t+k|t)分析：.平稳性：y(t)=c(q*1)/A(q*1)e(t)=e(t-j)可逆性：e(t)=A(q*1)/C(q*1)y(t)=y(t-j)故：L(y(t)y(t-1)=L(e(t)e(t-1)y(t+k|t)=proj(y(t+k)|y(t)y(t-1)=proj(y(t+k)|e(t)e(t-1),历史：1 wiener-kdmogorov 预报方法（1940、ARMA(p、q)=MACLO）2 Box-Jenkins 递推映射新方法（1970年）射影.3 Astrom预报方法.4 kalman预报方法（1960年）.,一步预报：.y(t+1)=-y(t)-y(t-1)-y(t+1-).+e(t+1)+e(t)+e(t+1-)y(t+1|t)=-y(t)-y(t-1)-y(t+1-).+e(t)+-e(t+1-)两式相减得 y(t+1)-y(t+1|t)=e(t+1)(即对于平稳ARMA过程,e(t)为y(t)的新息),定理:Box-Jenkins逆推预报器.K级预报.y(t+k)=-y(i+k-i)+e(t+k-i)(=1)1 k y(t+k|t)=-y(t+k-i|t)+e(t+k-i)k y(t+k|t)=-y(t+k-i|t).规定:y(t+k-i|t)=y(t+k-i)(t+k-i t),例1(1-aq*1)y(t)=e(t)(AR(1)|a|1.求 y(t+k|t)=?.解：y(t)=ay(t-1)+e(t).y(t+1)=ay(t)+e(t+1).y(t|n)=ay(t).y(t+k)=ay(t+k-1)+e(t+k)(k 1).y(t+k|t)=ay(t+k-1|t)=a*k y(t).例2 y(t)=(1+(q*1)e(t)|c|1.解：y(t+1)=e(t+1)+ce(t).y(t+1|t)=ce(t)=y(t).y(t+1|t)=0(k1).y(t+1|t)+c y(t|t-1)=cy(t).初值：y(|1).,例3：ARMA（1.1)解：(1aq*-1)y(t)=(1cq-1)e(t)求y(t+k|t)=?y(t+1)=ay(t)+e(t+1)+ce(t)y(t+1|t)=ay(t)+ce(t)y(t+k|t)=ay(t+k-1|t)(k1)非递推：y(t+k|t)=a*k-1(ay(t)+ce(t)e(t)=y(t)y(t+k|t)=a*k-1ay(t)+y(t)=a*(k-1)y(t)=y(t),注6：稳态Kalman滤波器写成Wienet的滤波器 x(t+1|t+1)=-kH x(t|t)+ky(t+1)=-kHx(t+1|t+1)=x(t|t)+ky(t+1)x(t+1|t+1)=q-1 x(t+1|t+1)+ky(t+1)-q-1 x(t+1|t+1)=ky(t+1)传递少数形式：x(t|t)=-q-1 ky(t)例1:x(t+1)=0.5x(t)+w(t)(1)Y(t)=x(t)+v(t)(2)求：y(t+2|t)解1：y(t+2)=x(t+2)+v(t+2)y(t+2|t)=x(t+2|t)=2 x(t|t)2-0.25-1=0 1=1.1328(1.1328+1)-1 2=-0.8828舍 k=H(HHT+R)1=1.1328(1.1328+1)-1=0.5311 K 0.5311y(t+2|t)=0.25-y(t)=0.25-y(t)1-q-1 1-q-110.5311 0.5 0.13275=-y(t)10.2344 q-1,解2：(10.5 q*-1)y(t)=w(t1)+(10.5 q*-1)v(t)(10.5 q*-1)y(t)=w(t1)+v(t)-0.5v(t1),1+1+0.25=(1+d12)*2 d*2+4.5d+1=0 d=-0.2344-0.5=d1 2 或d=-4.2656(舍)a(a+c)0.5(0.5-0.2344)y(t+2/t)-y(t)=-y(t)1+c q*-1 10.2344 q*-1 0.1328=-y(t)10.2344 q*-1,trm预报器例1：一个启发性的例子 ARMA（1.1）（1aq-1）y(t)=(1cq-1)e(t)a 1 c 1 求y(t+2/t)=?分析：1+cq*-1 y(t+2)=-e(t+2)1-aq*-1除法综合：1+(c+a)q*-1,1-aq*-1 1+cq*-1 1-aq*-1,(c+a)q*-1(c+a)q*-1a(c+a)q*-2,a(c+a)q*-2 a(c+a),y(t+2)=(1+(c+a)q*-1-)e(t+2)1-aq*-1,a(c+a)=e(t+2)+(c+a)e(t+1)+-e(t),1-aq*-1,未知 已知,a(c+a)y(t+2|t)=-e(t)1-aq*-1 a(c+a)1-aq*-1=-.-y(t)1-aq*-1 1+cq*-1y(t+2|t)=y(t+2)-y(t+2|t)=e(t+2)+(a+c)e(t+1)Ey*2(t+2|t)=1+（a+c）*2*2,一般情况：A(q*-1)y(t)=c(q*-1)e(t)(平稳、可逆)求：y(t+k|t)=?c(q*-1)分析：y(t+k)=-e(t+k)A(q*-1)c(q*-1)G(q*-1)综合除法:-=F(q-1)+q*-k-A(q*-1).F(q*-1)=+q*-1+q*-（k-1）G(q*-1)=+q*-1+q*-ng,Diophantine方程：c(q*-1)=A(q*-1)F(q*-1)+q*-k G(q*-1)=max（na-1,nc-k）两边比较系数可求得F，G G(q*-1)y(t+k)=F(q*-1)e(t+k)+-e(t)A(q*-1)G(q*-1)trm预报器 y(t+k|t)=-y(t)C(q*-1)y(t+k|t)=y(t+k)-y(t+k|t)=F(q*-1)e(t+k)Ey*2(t+k|t)=*2*2递推trm预报器：C(q*-1)y(t+k|t)=G(q*-1)y(t),例2:ARMA(1,1)(1-aq*-1)y(t)=(1+cq*-1)e(t)|ak|ck|求 y(t+k|t)=?解:F(q*-1)=+q*-1+q*(k-1)=max(1-1,1-k)=0 G(q*-1)=(1+cq*-1)=(1+aq*-1)(+q*-1+q*-(k-1)+(a*-k)比较系数法：1=1 c=-a=c+q 0=-a=a(c+a)0=-a=a*(k-2)(a+c)0=-a=a*(k-1)(a+c)y(t+k|t)=y(t)例平稳可逆ARMA(1,2)(1-aq*-1)y(t)=(1+q*-1+q*-2)e(t)求一步trmy(t+1|t)=?及一步Box-Jenkins y(t+1|t)=?,解：k=1 F(Q*-1)=max(-k,-1)=max(2-1,1-1)=1 G(q*-1)=+q*-1Diophantine方程1+q*-1+q*-2=(1-aq*-1)+q*-1(+q*-1)1=1=-a+=+a=y(t+1|t)=y(t)Box-Jenkins方法:y(t+1|t)=ay(t)+e(t)+e(t-1)e(t)=y(t)e(t-1)=y(t-1)=y(t)y(t+1|t)=y(t),=y(t)二者等价封闭形式稳态Kalman滤波器 x(t+1|t+1)=In-KH x(t|t)+ky(t)传递函数形式 x(t|t)=In-q*-1*(-1)ky(t)封闭形式稳态Kalman预报器 x(t+1|t)=In-KH x(t|t-1)+y(t)传递函数形式 x(t+1|t)=In-q*-1 y(t)新息形式下稳态Kalman预极器 x(t+1|t)=x(t|t)x(t|t)=x(t|t-1)+k x(t+1|t)=x(t|t-1)+k=x(t|t-1)+传递函数形式 x(t+1|t)=In-q*(-1)*(-1),第十五讲 总复习,.基本概念 时间序列 时间序列预极 样本函数(实现)Box-Jenkins 时间序列分析方法:一个样本函数(一个实现)建模 估计 trm 时间序列滤波 Kalman 平稳随机过程:=0=0 r(i)=E ARMA模型(q*-1)=Q(q*-1)平稳性(x)=0的根在单位圆外(q*-1)=1-q*-1+-q*-p 可逆性(x)=0的根在单位圆外(q*-1)=1-q*-1-q*-p 相关函数=Ez(t)z(t-k)=|E=0 E=标准相关函数:方差:=,.ARMA过程展式 成形滤波比(规定:=1,=0(jq)白压滤波比(规定:=1,=0(jp)=0(jp)注意:=0 j=0 jARMA(p,q)=MA()=MA()充分大 AR()=AR()充分大也可用几何级数:(1-q*-1)z(t)=,.相关函数的计算.最小二乘法:y(t)-y(t-1)-y(t-2)-y(t-n)=y(t)=+LS结构RLS公式:=0 p(0)=I=10*5,RELS A(q*-1)y(t)=D(q*-1)A(q*-1)=1-q*-1-q*-D(q*-1)=1-q*-1-q*-定义:=y(t-1)y(t-),-LS结构:y(t)=+未知,用估值来代替=y(t-1)y(t-)-偶合,白噪声估值=y(t)-与参数估值互偶RLS-RELS=y(t)=y(t-j)=y(t-j)(q*-1)y(t)=多维多重RLS.射影理论1.射影公式 proj(x|y)=x=Ex+(y-Ey)2.新息序列:=y(k)-proj(y(k)|y(k)=y(k)-y(k|k-1),是白噪声.(ij)L(y(1)y(k)L()3.递推射影公式 proj(x|y(1)y(k)=proj(x|y(k-1)y(1)+Ex E*(-1).Kalman滤波(稳态.非稳)x(t+1)=x(t)+w(t)x(t)R*n y(t)=Hx(t)+v(t)y(t)R*m最优Kalman滤波方程Riccati方程P(t+1|t)=P(t|t-1)-P(t|t-1)H*THP(t|t-1)H*T+R*(-1)HP(t|t-1)*T+PQP*T,稳态Kalman滤波 P(t|t-1)=稳态:=-H*T(HH*T+R)*-1H*T+PQP*T K=H*T(HH*T+R)*-1 P(t+1|t+1)=In-KH.预极.最优预极 1.Box-Jenkins递推预极比 L(y(t)y(1)=L(e(t)e(1)proj(e(t+i)|e(t)e(t-1)=0 i0 e(t+i)i0 proj(y(t+i)|y(t)=y(t+i|t)i0 y(t+i)i 0 2.trm预极比 A(q*-1)y(t)=c(q*-1)e(t)平稳可逆 y(t+1|t)=y(k)解 Diophantine方程 C(q*-1)=A(q*-1)F(q*-1)+q*(-k)G(q*-1)=max(-1,-k)G(q*-1)=+q*-1+q*-,F(q*-1)=+q*-1+q*-(k-1)预极误差:(t+k|t)=y(t+k)-y(t+k|t)=F(q*-1)e(t+k)预极误差方差:E=,