-
Notifications
You must be signed in to change notification settings - Fork 3
/
EnKF_l63.m
247 lines (217 loc) · 6.12 KB
/
EnKF_l63.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
function A = enkf_l63(E)
%% EnKF_l63.m
%
% Run the Lorenz 1963 model and
% assimilate observations using the Ensemble Kalman Filter
% Lisa Neef; started 17 June 2013
%
% INPUT:
% E: a Matlab structure that holds all the model and assimilation
% parameters. It is generated using E = set_enkf_inputs
%
% OUTPUT:
% A: a Matlab structure that holds the truth,observed, and analysis values
% for the model variables x,y,and z
% A.xt
% A.yt
% A.zt
% A.xo
% A.yo
% A.zo
% A.xa
% A.ya
% A.za
% as well as the average state error (i.e. the true error)
% A.ETave
% and the average analysis (or estimated) error
% A.EAave
%-----------------------------------------------------------
%% extract the individual assimilation parameters from E
sigma = E.sigma;
rho = E.rho;
beta = E.beta;
dt = E.dt;
xt0 = E.xt0;
xf0 = E.xf0;
sig0 = E.sig0;
sig_obs = E.sig_obs;
obsx = E.obsx;
obsy = E.obsy;
obsz = E.obsz;
N = E.N;
Tend = E.Tend;
tobs = E.tobs;
obs_meanxy = E.obs_meanxy;
obs_meanyz = E.obs_meanyz;
obs_meanxz = E.obs_meanxz;
%% Initialize arrays to hold everything
t = 1:dt:Tend;
nT = length(t);
XT = zeros(3,nT)+NaN; % array to hold the true state in time
XENS = zeros(N,3,nT)+NaN; % array to hold the ensemble
S = zeros(3,nT)+NaN; % array to hold analysis error variance in time
YOBS = zeros(3,nT)+NaN;
%% initial conditions
XT(:,1) = xt0; % initial truth
xf = xt0+sig0*rand(3,1); % initial forecast
for iens = 1:N
XENS(iens,:,1) = xf0 + sig0*randn(3,1);
end
%% define the observation error covariance matrix
R = sig_obs*eye(3);
%% define the observation operator, based on the desired obs configuration
H = zeros(3);
if obsx, H(1,1) = 1; end
if obsy, H(2,2) = 1; end
if obsz, H(3,3) = 1; end
% note that this code is not (yet) equipped to handle cases were we observe single obs AND
% averages -- throw an error in this case
observe_averages = obs_meanxy+obs_meanyz+obs_meanxz;
observe_single_variables = obsx+obsy+obsz;
if (observe_averages > 0) && (observe_single_variables > 0)
error('Cannot observe both variable averages and single variables...please change the input file.')
end
if observe_averages > 0
if obs_meanxy
H(1,1) = 0.5;
H(1,2) = 0.5;
end
if obs_meanyz
H(2,2) = 0.5;
H(2,3) = 0.5;
end
if obs_meanxz
H(3,1) = 0.5;
H(3,3) = 0.5;
end
end
%% Loop in time
xfens = squeeze(XENS(:,:,1));
for k = 1:nT-1
if E.run_filter
% are there observations?
if (mod(t(k),tobs) == 0)
% create observations of the true state
YOBS(:,k) = H*XT(:,k)+H*sig_obs*rand(3,1);
% get the forecast error covariance matrix from the ensemble
D = zeros(3,N);
for iens = 1:N
D(:,iens) = xfens(iens,:) - mean(xfens,1);
end
Pf1 = (1/N)*D*D';
% localize the covariance matrix?
if E.localize
Pf = eye(3).*Pf1;
else
Pf = Pf1;
end
% update the entire ensemble with the observations
K_bottom = H*Pf*H' + R;
K_top = H*Pf';
K = K_top * inv(K_bottom);
for iens = 1:N
% for each ensemble member, create a perturbed observation vector
yens = YOBS(:,k)+H*sig_obs*rand(3,1);
XENS(iens,:,k) = xfens(iens,:)' + K*(yens - H*xfens(iens,:)');
end
else
% if no observation, then the forecast becomes the analysis
XENS(:,:,k) = xfens;
end
else
% if not running the filter, then the forecast is the analysis
XENS(:,:,k) = xfens;
end
% regardless of whether there's been an observation, the analysis error covariance matrix comes from the ensemble
D = zeros(3,N);
for iens = 1:N
D(:,iens) = XENS(iens,:,k) - mean(XENS(:,:,k),1);
end
Pa = (1/N)*D*D';
% save the diagonals of the analysis error covariance matrix -- these are the variances
S(:,k) = diag(Pa);
% evolve the truth forward
XT(:,k+1) = lorenz63(XT(:,k), sigma, rho, beta, dt);
% evolve the analysis ensemble forward to become the next forecast
xfens = zeros(N,3);
for iens = 1:N
xfens(iens,:) = lorenz63(XENS(iens,:,k), sigma, rho, beta, dt);
end
end
%---------------PLOTTING----------------------------------
%% Compute a few other output quantities
XA = squeeze(mean(XENS,1));
%% Produce Plots!
YL = {'x','y','z'};
%% definte some plot settings
LW = 2; % line width
tcol = [0,0,0];
acol = [217,95,2]/256.0;
acol = [102,166,30]/256.0;
ocol = [241,41,138]/256.0;
ecol = .7*ones(1,3);
% plot the state analysis versis truth
figure(1),clf
h = zeros(1,4); % legend handle
T = ones(N,1)*t;
for ic = 1:3
subplot(3,1,ic)
ens = transpose(squeeze(XENS(:,ic,:)));
dum = plot(transpose(T),ens,'Color',ecol,'LineWidth',1);
hold on
h(1) = dum(1);
h(2) = plot(t,XT(ic,:),'Color',tcol,'LineWidth',LW);
h(3) = plot(t,XA(ic,:),'Color',acol,'LineWidth',LW);
h(4) = plot(t,YOBS(ic,:),'o','Color',ocol,'MarkerSize',5,'LineWidth',LW);
if ic == 1
title('Lorenz 1963 Model - State Variables')
end
if ic == 3
xlabel('time')
legend(h, 'ensemble','truth','analysis','obs','Location','SouthOutside','Orientation','Horizontal')
end
ylabel(YL(ic))
end
% plot the estimated and true errors
ET = sqrt((XT-XA).^2);
EA = sqrt(S);
figure(2),clf
for ic = 1:3
subplot(3,1,ic)
h = zeros(1,2);
h(1) = semilogy(t,ET(ic,:),'Color',tcol,'LineWidth',LW);
hold on
h(2) = semilogy(t,EA(ic,:),'Color',acol,'LineWidth',LW);
ylabel(YL(ic))
if ic == 1
title('Lorenz 1963 Model - RMSE')
end
if ic == 3
xlabel('time')
legend(h, 'true RMSE','estimated RMSE','Location','SouthOutside','Orientation','Horizontal')
end
end
% average the errors for each variable over the integration time, and print to the screent
ETave = nanmean(ET,2);
EAave = nanmean(EA,2);
names = {'x','y','z'};
disp('+++++++Assimilation Run Average Errors++++++++++')
for ic = 1:3
str_out = strcat(names(ic),': True Error = ',num2str(ET(ic),3),' Estimated = ',num2str(EA(ic),3));
disp(str_out)
end
disp('average state error')
str_out = strcat('True Error = ',num2str(mean(ETave),3),' Estimated = ',num2str(mean(EAave),3));
disp(str_out)
% pack everything into the output structure
A = struct('xt',XT(1,:),...
'yt',XT(2,:),...
'zt',XT(3,:),...
'xa',XA(1,:),...
'ya',XA(2,:),...
'za',XA(3,:),...
'xo',YOBS(1,:),...
'yo',YOBS(2,:),...
'zo',YOBS(3,:),...
'ETave',ETave,...
'EAave',EAave);