rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
lin_ds_detailed_simulation.m
1 function sim_data = lin_ds_detailed_simulation(model,model_data)
2 %function sim_data = lin_ds_detailed_simulation(model,model_data)
3 %
4 % function performing a general time evolution of a linear dynamical
5 % system
6 %
7 % @f[ \frac{d}{dt} x = A(t,\mu) x + B(t,\mu) u @f]
8 % @f[ y = C(t,\mu) x + D(t,\mu) u \quad \mbox{ for }
9 % \quad t \in [0,T] @f]
10 %
11 % @f[ x(0) = x_0(\mu) @f]
12 %
13 % by theta-scheme time discretization.
14 %
15 %
16 % required fields of model:
17 % T : final time
18 % nt : number of time-intervals until T, i.e. nt+1
19 % solution slices are computed
20 % theta : real value in [0,1]
21 % - 0 = explicit
22 % - 0.5 = Crank-Nicolson
23 % - 1 = Implicit discretization
24 % x0_function_ptr: pointer to function for computing the
25 % initvalues-DOF with arguments (model)
26 % A_function_ptr: pointer to function for computing the matrix A
27 % B_function_ptr: pointer to function for computing the matrix B
28 % C_function_ptr: pointer to function for computing the matrix C
29 % D_function_ptr: pointer to function for computing the matrix D
30 % u_function_ptr: pointer to input function u(t)
31 % data_const_in_time : if this optional field is 1, the time
32 % evolution is performed with constant operators,
33 % i.e. only the initial-time-operators are computed
34 % and used throughout the time simulation.
35 % affinely_decomposed: if this optional field is set, model_data is
36 % assumed to contain the components of the matrices and online
37 % assembly is performed by coefficient-component linear combinaton
38 % save_time_indices: list of integer time indices between 0 (!) and nt,
39 % which are to be computed, stored and returned
40 % verbose: integer verbosity level
41 %
42 % Required fields of model_data:
43 % A_components: cell array of affine decomposed A-components
44 % B_components: cell array of affine decomposed B-components
45 % C_components: cell array of affine decomposed C-components
46 % x0_components: cell array of affine decomposed x0-components
47 % Required fields of model_data:
48 %
49 % Generated fields of sim_data:
50 % X: matrix containing state vectors of solution trajectory as
51 % columns
52 % Y: matrix containing output vectors of trajectory as columns
53 % time: vector of time-values in [0,T] for columns in X,Y
54 % time_indices: integer time step numbers for columns in X,Y
55 %
56 % Return values:
57 % sim_data: structure with simulation data results
58 
59 % Bernard Haasdonk 2.4.2009
60 
61 if model.verbose >= 9
62  disp('entered detailed simulation ');
63 end;
64 
65 %model.affine_decomp_mode = 'complete';
66 model.decomp_mode = 0;
67 model = model.set_time(model,0);
68 %model.t = 0;
69 
70 if isfield(model,'save_time_indices')
71  save_time_index = zeros(1,model.nt+1);
72  save_time_index(model.save_time_indices(:)+1) = 1;
73  save_time_indices = find(save_time_index==1)-1;
74 else
75  save_time_index = ones(1,model.nt+1);
76  save_time_indices = 1:(model.nt+1);
77 end;
78 
79 % initial values by midpoint evaluation
80 
81 if ~model.affinely_decomposed
82  xt = model.x0_function_ptr(model,model_data);
83  C0 = model.C_function_ptr(model,model_data);
84 else
85  old_decomp_mode = model.decomp_mode;
86  model.decomp_mode = 2;
87  x0_coeff = model.x0_function_ptr(model,[]);
88  xt = lincomb_sequence(model_data.x0_components,x0_coeff);
89  C0_coeff = model.C_function_ptr(model,[]);
90  C0 = lincomb_sequence(model_data.C_components,C0_coeff);
91  model.decomp_mode = old_decomp_mode;
92 end;
93 D = model.D_function_ptr(model,model_data); % no decomposition required
94 X = zeros(length(xt(:)),length(save_time_indices));
95 Y = zeros(size(C0,1),length(save_time_indices));
96 time_indices = zeros(1,length(save_time_indices));
97 time = zeros(1,length(save_time_indices));
98 
99 t_column = 1; % next column index to be filled in output
100 t_ind = 0; % current t_ind between 0 and nt
101 t = 0; % absolute time between 0 and T
102 
103 % store init data
104 if save_time_index(t_ind+1)
105 % disp('\n OK, saving initial time')
106 % keyboard;
107  X(:,t_column) = xt;
108  time_indices(t_column) = t_ind;
109  time(t_column) = t;
110  u = model.u_function_ptr(model);
111  Y(:,t_column) = C0 * xt + D * u;
112  t_column = t_column + 1;
113 end;
114 
115 model.dt = model.T/model.nt;
116 dt = model.dt;
117 theta = model.theta;
118 operators_required = 1;
119 previous_implicit_operators_available = 0;
120 
121 % loop over time-steps
122 for t_ind = 1:model.nt
123  t = model.dt * (t_ind-1); % compute new x and y at time t
124  if model.verbose >= 10
125  disp(['entered time-loop step ',num2str(t_ind)]);
126  end;
127  if model.verbose >= 9
128  fprintf('.');
129  end;
130  % get matrices and bias-vector
131  model = model.set_time(model,t);
132 % model.t = t;
133  if operators_required
134  if ~model.affinely_decomposed
135  if previous_implicit_operators_available
136  A0 = A1;
137  B0 = B1;
138  else
139  A0 = model.A_function_ptr(model,model_data);
140  B0 = model.B_function_ptr(model,model_data);
141  end;
142  % C and D required at time t+dt;
143  model = model.set_time(model,t+model.dt);
144 % model.t = t+model.dt;
145  C1 = model.C_function_ptr(model,model_data);
146  D1 = model.D_function_ptr(model,model_data);
147  if theta>0
148  A1 = model.A_function_ptr(model,model_data);
149  B1 = model.B_function_ptr(model,model_data);
150  previous_implicit_operators_available = 1;
151  end;
152  model = model.set_time(model,t);
153  % model.t = t;
154  else
155  model.decomp_mode = 2;
156  if ~previous_implicit_operators_available
157  A0_coeff = model.A_function_ptr(model,[]);
158  B0_coeff = model.B_function_ptr(model,[]);
159  end;
160  model = model.set_time(model,t+model.dt);
161  if theta>0
162  A1_coeff = model.A_function_ptr(model,[]);
163  B1_coeff = model.B_function_ptr(model,[]);
164  end;
165  %model.t = t+model.dt;
166  C1_coeff = model.C_function_ptr(model,[]);
167  model.decomp_mode = old_decomp_mode;
168  D1 = model.D_function_ptr(model,[]); % no decomposition of D
169  model = model.set_time(model,t);
170  %model.t = t;
171  if previous_implicit_operators_available
172  A0 = A1;
173  B0 = B1;
174  else
175  A0 = lincomb_sequence(model_data.A_components,A0_coeff);
176  B0 = lincomb_sequence(model_data.B_components,B0_coeff);
177  end;
178  if theta>0
179  A1 = lincomb_sequence(model_data.A_components,A1_coeff);
180  B1 = lincomb_sequence(model_data.B_components,B1_coeff);
181  end;
182  C1 = lincomb_sequence(model_data.C_components,C1_coeff);
183  end;
184 
185  if model.data_const_in_time
186  operators_required = 0;
187  end;
188  end;
189  u0 = model.u_function_ptr(model);
190  if theta == 0
191  xt = xt + dt * A0 * xt + dt * B0 * u0;
192  else
193  % get input at next time
194  model = model.set_time(model,t+model.dt);
195  u1 = model.u_function_ptr(model);
196  model = model.set_time(model,t);
197  % solve M1 x1= rhs := M0 x0 + b
198  M1 = speye(size(A0))-dt * theta * A1;
199  M0 = speye(size(A0))+dt * (1-theta) * A0;
200  b = dt * ( theta * B1 * u1 + (1-theta)*B0 *u0);
201  rhs = M0 * xt + b;
202 % keyboard;
203  xt = M1 \ rhs;
204  end;
205  % X(:,t_ind+1) = X(:,t_ind) + model.dt * (A * X(:,t_ind) + B * u);
206  %Y(:,t_ind+1) = C * X(:,t_ind) + D * u;
207 
208  if save_time_index(t_ind+1)
209 % disp('\n OK, saving next time');
210 % keyboard;
211  X(:,t_column) = xt;
212  time_indices(t_column) = t_ind;
213  % y evaluation requires + dt !!
214  model = model.set_time(model,t+model.dt);
215  %model.t = t + model.dt;
216  u = model.u_function_ptr(model);
217  Y(:,t_column) = C1 * xt + D1 * u;
218  time(t_column) = t;
219  t_column = t_column + 1;
220  fprintf('.');
221  end;
222 
223 end;
224 
225 if model.verbose >= 9
226  fprintf('\n');
227 end;
228 
229 sim_data = [];
230 sim_data.X = X;
231 sim_data.Y = Y;
232 sim_data.time = time;
233 sim_data.time_indices = time_indices;
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
Definition: verbose.m:17