1 function sim_data = lin_ds_detailed_simulation(model,model_data)
2 %
function sim_data = lin_ds_detailed_simulation(model,model_data)
4 %
function performing a general time evolution of a linear dynamical
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]
11 % @f[ x(0) = x_0(\mu) @f]
13 % by theta-scheme time discretization.
16 % required fields of model:
18 % nt : number of time-intervals until T, i.e. nt+1
19 % solution slices are computed
20 % theta : real value in [0,1]
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
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:
49 % Generated fields of sim_data:
50 % X: matrix containing state vectors of solution trajectory as
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
57 % sim_data: structure with simulation data results
59 % Bernard Haasdonk 2.4.2009
62 disp(
'entered detailed simulation ');
65 %model.affine_decomp_mode =
'complete';
66 model.decomp_mode = 0;
67 model = model.set_time(model,0);
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;
75 save_time_index = ones(1,model.nt+1);
76 save_time_indices = 1:(model.nt+1);
79 % initial values by midpoint evaluation
81 if ~model.affinely_decomposed
82 xt = model.x0_function_ptr(model,model_data);
83 C0 = model.C_function_ptr(model,model_data);
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;
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));
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
104 if save_time_index(t_ind+1)
105 % disp('\n OK, saving initial time')
108 time_indices(t_column) = t_ind;
110 u = model.u_function_ptr(model);
111 Y(:,t_column) = C0 * xt + D * u;
112 t_column = t_column + 1;
115 model.dt = model.T/model.nt;
118 operators_required = 1;
119 previous_implicit_operators_available = 0;
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
125 disp(['entered time-loop step ',num2str(t_ind)]);
130 % get matrices and bias-vector
131 model = model.set_time(model,t);
133 if operators_required
134 if ~model.affinely_decomposed
135 if previous_implicit_operators_available
139 A0 = model.A_function_ptr(model,model_data);
140 B0 = model.B_function_ptr(model,model_data);
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);
148 A1 = model.A_function_ptr(model,model_data);
149 B1 = model.B_function_ptr(model,model_data);
150 previous_implicit_operators_available = 1;
152 model = model.set_time(model,t);
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,[]);
160 model = model.set_time(model,t+model.dt);
162 A1_coeff = model.A_function_ptr(model,[]);
163 B1_coeff = model.B_function_ptr(model,[]);
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);
171 if previous_implicit_operators_available
175 A0 = lincomb_sequence(model_data.A_components,A0_coeff);
176 B0 = lincomb_sequence(model_data.B_components,B0_coeff);
179 A1 = lincomb_sequence(model_data.A_components,A1_coeff);
180 B1 = lincomb_sequence(model_data.B_components,B1_coeff);
182 C1 = lincomb_sequence(model_data.C_components,C1_coeff);
185 if model.data_const_in_time
186 operators_required = 0;
189 u0 = model.u_function_ptr(model);
191 xt = xt + dt * A0 * xt + dt * B0 * u0;
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);
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;
208 if save_time_index(t_ind+1)
209 % disp('\n OK, saving next time');
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;
219 t_column = t_column + 1;
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...