1 / 42

建模与估计讲义（二） 主讲人：王欣

建模与估计讲义（二） 主讲人：王欣. 第三节 新息序列 (Innovation sequence). 定义 : 设 ( 相关随机序列正交化 ) 有二阶矩 ( 数学期望、方差都　　）的随机序列，定义它的新息序列为： 定义 : 基于 k-1 以前信息对 y(k) 的线性最小方差估计 则 规定 : 定理 1 ：新息序列是正交序列（白噪声）

gaetan
Télécharger la présentation

建模与估计讲义（二） 主讲人：王欣

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 建模与估计讲义（二） 主讲人：王欣

  2. 第三节 新息序列 (Innovation sequence) 定义:设 (相关随机序列正交化) 有二阶矩(数学期望、方差都　　）的随机序列，定义它的新息序列为： 定义: 基于k-1以前信息对y(k)的线性最小方差估计 则 规定: 定理1：新息序列是正交序列（白噪声） 证： 射影具有无偏性 i≠j不妨设i〉j 根据射影性质

  3. 新息实质上是摄影误差 定理2: 证明: 显然 注意 事实上: 重要意义:新息序列 与原序列y(k)含有相同的统计信息

  4. 定理3:递推射影公式 证明:

  5. 第四节 kalman滤波 考虑线性离散定常随机系统 ① ② 其中状态x(t) ∈ ,观测y(t) ∈ ,w(t) ∈ 输入噪声.观测噪声v(t) ∈ 假设1: w(t)和v(t)是零均值,方差阵各为Q和R的独立白噪声 即:

  6. 假设2: x(0)与w(t)和v(t)相互独立 Kalman滤波问题是基于观测 (y(1)…y(t)) 求状态 x(j)的线性最小方差估计 极小化 即求 滤波 平滑 预报

  7. 分析:有递推射影公式 k(t+1)称为滤波增益 对①取射影 因为 事实上

  8. 新息: 引入定义:

  9. 即 得 即

  10. 即 注意:

  11. 简化: 即: 总结:考虑随机系统 问题:基于观测 求

  12. Kalman Filter 初值 协方差阵

  13. 第四节 回忆:考虑随机系统 问题:基于观测 求 定理:Kalman Filter

  14. 初值: 协方差阵 例1:已知 及 求 解: 注1:Kalman滤波的和预报器的不同形式 封闭形式:

  15. 滤波器 预报增益 预报器

  16. 注2: Riccati方程 即: 初值: 注3:时变系统Kalman方程组同样适用

  17. 有Kalman方程组 初值: 注4:Kalman滤波与RLS关系(RLS是K的特例) 考虑AR(N)

  18. 可写成状态空间模型 应用Kalman Filter 定义： 具有RLS公式

  19. 注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 rank[p 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

  20. 且有关系 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)

  21. x(t+1|t+1)= x(t|t)+ ky(t+1) =(1-kh) x(t+1|t) = x(t|t)-1+ y(t) =-kph = k Riccati方程 ∑= [∑-∑h（h∑h+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*2qr h*2 ∑*2+(r- *2r-b*2h*2q) ∑-b*2qr=0 作一元二次方程，取∑〉0为GT为

  22. 第五讲 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)

  23. 历史：1 wiener-kdmogorov 预报方法（1940、ARMA(p、q)=MACLO）2 Box-Jenkins 递推映射新方法（1970年）射影. 3 Astrom预报方法. 4 kalman预报方法（1960年）. Box-Jenkins递推预报器. 考虑平稳可逆的ARMA过程y(t) . A(q*－1 )y(t)=c(q*－1 )e(t) 其中e(t)是白噪声：E e(t)=0 E[e(t)e(s)]= *2 A(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) …)

  24. 一步预报：. 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)

  25. 例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 (k>1) . y(t+1|t)+c y(t|t-1)=cy(t) . 初值：y( | –1) .

  26. 例3：ARMA（1.1) 解：(1－aq*-1)y(t)=( 1－cq-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) (k>1) 非递推：y(t+k|t)=a*k-1 (ay(t)+ce(t) e(t)= y(t) y(t+k|t)= a*k-1[ay(t)+ y(t) ] = a*(k-1)[ ]y(t) = y(t)

  27. 注6：稳态Kalman滤波器写成Wienet的滤波器 x(t+1|t+1)=〔㏑-kH〕 x(t|t)+ky(t+1) =〔㏑-kH〕 ∴x(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(H∑HT+R)－1=1.1328(1.1328+1) -1 =0.5311 K 0.5311 ∴y(t+2|t)=0.25---------y(t)=0.25----------------------y(t) 1- q-1 1- q-1〔1－0.5311〕0.5 0.13275 =----------------y(t) 1－0.2344 q-1

  28. 解2：(1－0.5 q*-1)y(t)= w(t－1)+ (1－0.5 q*-1) v(t) (1－0.5 q*-1)y(t)= w(t－1)+ v(t)-0.5v(t－1) 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 1－0.2344 q*-1 0.1328 =---------------------- y(t) 1－0.2344 q*-1

  29. ÅtrÖm预报器 例1：一个启发性的例子 ARMA（1.1） （1－aq-1）y(t)=( 1－cq-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

  30. 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*-1 y(t+2|t)= y(t+2)- y(t+2|t)=e(t+2)+(a+c) e(t+1) E[y*2(t+2|t)]=[1+（a+c）*2] *2

  31. 一般情况：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

  32. 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) ∴ÅtrÖm预报器 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) E[y*2(t+k|t)]= *2∑ *2 递推ÅtrÖm预报器：C(q*-1) y(t+k|t)= G(q*-1) y(t)

  33. 例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) 求一步ÅtrÖmy(t+1|t)=? 及一步Box-Jenkinsy(t+1|t)=?

  34. 解：k=1 F(Q*-1)= =max( -k, -1)=max(2-1,1-1)=1 G(q*-1)= + q*-1 Diophantine方程 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)

  35. = 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)

  36. 第十五讲 总复习 Ⅰ.基本概念 时间序列 时间序列预极 样本函数(实现) Box-Jenkins 时间序列分析方法:一个样本函数(一个实现) 建模 估计ÅtrÖm 时间序列滤波 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 相关函数 =E[z(t)z(t-k)]= | |≤ E =0 E = 标准相关函数: 方差: =

  37. Ⅱ.ARMA过程展式 成形滤波比 (规定: =1, =0 (j<0) =0 (j>q)) 白压滤波比 (规定: =1, =0 (j>p) =0 (j>p) 注意: =0 j =0 j ARMA(p,q)=MA( )=MA( ) 充分大 AR( )=AR( ) 充分大 也可用几何级数: (1- q*-1)z(t)=

  38. Ⅲ.相关函数的计算 Ⅳ.最小二乘法: y(t)- y(t-1)- y(t-2)…- y(t-n)= y(t)= + LS结构 RLS公式: =0 p(0)=αI α=10*5

  39. 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)

  40. 是白噪声. ⊥ (i≠j) L(y(1)…y(k)) ⊥L( … ) 3.递推射影公式 proj(x|y(1)…y(k))=proj(x|y(k-1))y(1)+E[x ]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*T[HP(t|t-1)H*T+R]*(-1)HP(t|t-1)] *T+PQP*T

  41. 稳态Kalman滤波 P(t|t-1)=∑ 稳态:∑= [∑-∑H*T(H∑H*T+R)*-1H∑] *T+PQP*T K= ∑H*T(H∑H*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 i>0 e(t+i) i≤0 proj(y(t+i)|y(t)…)= y(t+i|t) i>0 y(t+i) i ≤0 2.ÅtrÖm预极比 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*-

  42. 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[ ]=

More Related