数学建模社区-数学中国

标题: 卡尔曼算法的matlab程序 [打印本页]

作者: 工科男    时间: 2011-11-26 16:16
标题: 卡尔曼算法的matlab程序
求卡尔曼算法的matlab预测程序
作者: 厚积薄发    时间: 2011-11-28 10:48
matlab下面的kalman滤波程序
  s: L1 s9 N9 u8 `; q! k& lclear  N=200; w(1)=0;   4 c5 S8 G! S, \4 K" l
w=randn(1,N) / S3 y; F: {! A& j# C) Q7 V
x(1)=0;   $ ~0 _# i* m" ]
a=1;   
  x, T0 U6 w, Cfor k=2:N;   2 F; ~9 Y# k5 w) o
x(k)=a*x(k-1)+w(k-1);   
7 n% C9 Y/ U$ d# i: b5 R; G+ Q+ aend   + q/ m6 p* c  ]2 Z
V=randn(1,N);   
- O5 S- v( L3 p6 D; f0 s- ]! pq1=std(V);   
* i3 b# o$ C( b) @' T1 _( M+ ~Rvv=q1.^2;   
# O& v- g+ I2 ]: n" @q2=std(x);   
0 Q( u: _2 {3 W0 B1 g4 jRxx=q2.^2;   
5 V8 C* b7 j( h& t/ ?q3=std(w);   
9 \! N. c7 l1 f, uRww=q3.^2;   - t! |* g9 H. N4 n/ h
c=0.2;   
* Q& u$ \6 ^- }0 w  _Y=c*x+V;   
' S7 m$ ~! r+ Q( ip(1)=0;   $ s; O6 {% o# @
s(1)=0;
* p6 e) s+ P1 _for t=2:N;   
' p) t8 m% K5 _  ^$ cp1(t)=a.^2*p(t-1)+Rww;   
) ~2 H/ _9 F: _4 h- B; Vb(t)=c*p1(t)/(c.^2*p1(t)+Rvv);   ; A3 C5 b0 \/ `( @& L) k+ w; ?
s(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1));   * g- t5 I! P% b) i% w
p(t)=p1(t)-c*b(t)*p1(t);   
1 s% b# C0 j) a/ K. M* B& g7 tend   ' l5 r7 M7 y7 A( E7 o7 l4 w" q
t=1:N;   ( R2 o  l: c+ {
plot(t,s,'r',t,Y,'g',t,x,'b');   * ^" O  z6 k0 o) p7 V1 E- v
function [x, V, VV, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, varargin)   8 D* l' `) C5 |3 c. w) U
% Kalman filter.   6 B; ?/ k6 U* G6 Z" b
% [x, V, VV, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, ...)   
' d2 ]. M: M9 Q$ C8 @; \%   
9 N! R! H/ s9 s$ `% INPUTS:   : ^+ o4 F7 v. |& q+ H- j; q
% y(:,t) - the observation at time t   $ z5 u. H, E: W' m
% A - the system matrix   5 F: V0 F- C1 F% ~8 i
% C - the observation matrix   ' `, }" r$ S. [9 f6 E* x  P
% Q - the system covariance   7 g- U, Z4 ?5 c* r/ {+ g2 J7 n
% R - the observation covariance   * A( ~9 ?  ]$ v4 S# O
% init_x - the initial state (column) vector   
* K' X, r: I  y) J% W+ r% init_V - the initial state covariance   8 x: @) S0 k0 D* C) X4 _) f' b7 |6 k1 J
%   " Y" M: I8 F4 {7 h, @
% OPTIONAL INPUTS (string/value pairs [default in brackets])   % F* i; j, ?# P1 @/ F; w9 y
% 'model' - model(t)=m means use params from model m at time t [ones(1,T) ]   
$ ?/ \- J, C" P( z* h% In this case, all the above matrices take an additional final dimension,   
# v$ ]& D. m* O4 A) N% i.e., A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m).   
* s. _* p. T5 O6 t+ j% However, init_x and init_V are independent of model(1).   & f0 _) Q0 L# _& \# O" Z1 i* j" A
% 'u' - u(:,t) the control signal at time t [ [] ]   
7 y1 G# f' w, S7 J( W% 'B' - B(:,:,m) the input regression matrix for model m   $ B9 m# P: U* q$ V$ z
%   
' c$ S! L3 ~% x" Z7 l; t) b/ ^% OUTPUTS (where X is the hidden state being estimated)   
8 ?) I+ y' G1 [8 Q/ r% x(:,t) = E[X(:,t) | y(:,1:t)]   
, {* l5 `" O" e" e6 c& P( F% V(:,:,t) = Cov[X(:,t) | y(:,1:t)]   
+ ]6 ]: S' ?7 q% VV(:,:,t) = Cov[X(:,t), X(:,t-1) | y(:,1:t)] t >= 2   - I+ l" W+ Z" x6 ]3 `6 h5 _8 ]+ w' g
% loglik = sum{t=1}^T log P(y(:,t))   / p0 J9 O9 m3 V- L' [0 L
%   - y5 Q" B6 V3 |, p: [/ o
% If an input signal is specified, we also condition on it:   
: d, t7 f- H# R! e  V3 G# d( Z" Y* p% e.g., x(:,t) = E[X(:,t) | y(:,1:t), u(:, 1:t)]   0 B0 a$ }6 O7 u
% If a model sequence is specified, we also condition on it:   
/ `! K6 y$ \/ u5 s4 r; s% z: [0 I% e.g., x(:,t) = E[X(:,t) | y(:,1:t), u(:, 1:t), m(1:t)]   0 H- {7 Z# J6 v) s& b
[os T] = size(y);   
5 e; |/ U9 Q/ U9 K4 C' w  B& q3 O0 ?! pss = size(A,1); % size of state space   
7 i9 f, h" s; ?5 Y% set default params   
. O2 G0 v2 [8 Q7 gmodel = ones(1,T);   5 N, o+ v( S, K& B4 L9 z
u = [];   4 L2 L# u0 G( ~' D8 x; E# y
B = [];   
8 a- ^& A, E) c3 bndx = [];   : h8 D! A& D' U/ }; }# o5 |* V
args = varargin;   $ I5 B# R0 \) b8 n: k/ a. _! S- O2 w
nargs = length(args);   
$ g/ ]) `0 b( m( |1 bfor i=1:2:nargs   1 W  {* ]8 N9 _. T5 N; X5 X
switch args   
' x' x4 d1 N5 t0 H$ t; T4 dcase 'model', model = args{i+1};   * ~3 c! D& k9 l: ^
case 'u', u = args{i+1};   
* J: M$ v2 p; ^$ N( scase 'B', B = args{i+1};   2 p. W3 T# w7 T$ O
case 'ndx', ndx = args{i+1};     D$ L8 u! f  @3 |+ |
otherwise, error(['unrecognized argument ' args])   3 W2 w! v; e2 G9 g$ T* l7 E8 d
end   
  Z# S+ _( F% J  |7 ^9 o+ {end   ' T7 F. k6 e2 {! U, U
x = zeros(ss, T);   
3 `) `* [# T6 t' ?4 o; z# \7 AV = zeros(ss, ss, T);   2 \  P. `8 T7 C! M
VV = zeros(ss, ss, T);   8 k6 y& F; F. _0 Y. u, g' A3 V1 p
loglik = 0;   1 R3 L; V# w; z
for t=1:T   m = model(t);   
' s% G* ~5 D% o' _, |7 _if t==1   %prevx = init_x(:,m);   % M* i% k' B: g& o/ c) e/ a8 @
%prevV = init_V(:,:,m);   
2 U" A& `: ~* F: {0 l: B& Eprevx = init_x;   ( x$ _! p1 h5 a6 z- C/ t4 g- Y
prevV = init_V;   
, ~9 o6 q3 J" y4 _% dinitial = 1;   / j  k. W/ X$ a9 u; n/ @
else   prevx = x(:,t-1);   - t8 @* \8 m+ Z4 `
prevV = V(:,:,t-1);   7 n  _4 _- f7 l- Q- G: s+ c1 u
initial = 0;   0 S0 J1 n! p1 _- h
end   
5 P" F0 J/ i# Q$ Z* b- Qif isempty(u)   
7 w* p# a. `( E+ _8 [[x(:,t), V(:,:,t), LL, VV(:,:,t)] = ...   + T( f. _% p" c, u3 N
kalman_update(A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m), y(:,t), prevx, prevV, 'initial', initial);   else   
9 I% q/ V# P0 s/ }  if isempty(ndx)   [x(:,t), V(:,:,t), LL, VV(:,:,t)] = ...   6 R" K% z! f* q# _3 d2 ^
     kalman_update(A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m), y(:,t), prevx, prevV, ...   'initial', initial,     'u', u(:,t), 'B', B(:,:,m));   - s$ b. Z5 B9 O0 X# `% ]' t
else   9 [+ ]- F: y, ^1 }# c) y! H8 ~( g
i = ndx;   
7 m- v! V( s. l% copy over all elements; only some will get updated   x(:,t) = prevx;   * \+ b, w" r9 n. q/ l. S# {* {
prevP = inv(prevV);   ( R( g' L+ L5 i* a! U& w( E! P: V4 V( e" E
prevPsmall = prevP(i,i);   
) i) {. J$ t6 vprevVsmall = inv(prevPsmall);   
6 h5 ~& g3 i2 W5 P- z  K[x(i,t), smallV, LL, VV(i,i,t)] = ...   kalman_update(A(i,i,m), C(:,i,m), Q(i,i,m), R(:,:,m), y(:,t), prevx(i), prevVsmall, ...   'initial', initial, 'u', u(:,t), 'B', B(i,:,m));   
$ M% G( P( E8 e5 |+ |3 JsmallP = inv(smallV);   
  ~2 b4 O' d& X% A8 q/ L" `prevP(i,i) = smallP;   ( T* `% i/ D; x9 W: B2 D
V(:,:,t) = inv(prevP);   
8 U* Q% r8 K2 J9 o: ]end   
& p4 H1 j+ H  {) Oend   
) q: b/ n" ?/ A* x3 |/ Cloglik = loglik + LL;   ( K3 F1 h* D/ T7 E( ^
end
作者: 剑游九天    时间: 2011-12-19 23:38
厚积薄发 发表于 2011-11-28 10:48
0 f: I7 |3 G5 tmatlab下面的kalman滤波程序
! J7 x+ f, F1 X% |) X2 {7 Y( ]+ Zclear  N=200; w(1)=0;   4 p8 C; [* p, \. Z- l9 D% x; I2 }) B( l
w=randn(1,N)

% T: Y% l( Y& N+ O. x: R8 ~9 x强大呀,下了
作者: 工科男    时间: 2012-2-5 10:24
太强大了
作者: 244190977    时间: 2013-10-7 19:59
太令人。。。。
作者: yulun9988    时间: 2014-1-13 21:45
谢谢。。。。。。。。。。。。。。。。。。。。。。。。。。
作者: yulun9988    时间: 2014-1-13 21:45





欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5