Mercurial > hg > ltpda

comparison m-toolbox/classes/@ssm/modifyTimeStep.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % MODIFYTIMESTEP modifies the timestep of a ssm object
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: MODIFYTIMESTEP modifies the timestep of a ssm object, and updates
5 % the A and B matrices supposing there is no aliasing.
6 %
7 % CALL: sys = modifyTimeStep(sys,pl)
8 %
9 % INPUTS:
10 % sys - (array of) ssm objects
11 % pl - A plist or numeric value giving new timestep value (param name 'newtimestep')
12 %
13 % OUTPUTS:
14 % sys - (array of) ssm
15 %
16 % <a href="matlab:utils.helper.displayMethodInfo('ssm', 'modifyTimeStep')">Parameters Description</a>
17 %
18 % VERSION: $Id: modifyTimeStep.m,v 1.11 2011/04/08 08:56:23 hewitson Exp $
19 %
20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21
22 function varargout = modifyTimeStep(varargin)
23
24 %% Check if this is a call for parameters
25 if utils.helper.isinfocall(varargin{:})
26 varargout{1} = getInfo(varargin{3});
27 return
28 end
29
30 %% send starting message
31 import utils.const.*
32 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
33
34 %% collecting input
35 in_names = cell(size(varargin));
36 for ii = 1:nargin,in_names{ii} = inputname(ii);end
37
38 [sys, ssm_invars, rest] = utils.helper.collect_objects(varargin(:), 'ssm', in_names);
39 [pl, invars2, rest] = utils.helper.collect_objects(rest(:), 'plist');
40 % override with a numeric input value
41 if ~isempty(rest)
42 % if the user inputs a single numerical value, we use it
43 if isnumeric(rest{1}) && numel(rest{1}) == 1
44 if ~isa(pl, 'plist') % perhaps we dont' have a plist yet
45 pl = plist();
46 end
47 pl.pset('newtimestep', rest{1}); % set newtimestep
48 else
49 pl = combine(pl, plist(rest{:}));
50 end
51 end
52 pl = combine(pl, getDefaultPlist());
53
54 %%% Internal call: Only one object + don't look for a plist
55 internal = strcmp(varargin{end}, 'internal');
56
57 %% dealing with user custom inputs
58
59 Nsys = numel(sys);
60
61 %% begin function body
62
63 %% building output depending on options
64 sys = copy(sys,nargout);
65 newtimestep = find(pl,'newtimestep');
66 outputAntiAlias = find(pl, 'outputAntiAlias');
67 if numel(newtimestep)~=1
68 error('parameter "newtimestep" should be of size 1x1')
69 end
70
71 for i = 1:Nsys
72 if ~sys(i).isnumerical
73 error(['error in modifTimeStep because system "',sys(i).name, '" should be numerical to be modified' ]);
74 end
75
76 %% retrieving input data
77 timestep = sys(i).timestep;
78 sssizes = sys(i).sssizes;
79 inputsizes = sys(i).inputsizes;
80
81 amat = ssm.blockMatFusion(sys(i).amats, sssizes, sssizes);
82 bmat = ssm.blockMatFusion(sys(i).bmats, sssizes, inputsizes);
83 if outputAntiAlias
84 outputsizes = sys(i).outputsizes;
85 cmat = ssm.blockMatFusion(sys(i).cmats, outputsizes, sssizes);
86 dmat = ssm.blockMatFusion(sys(i).dmats, outputsizes, inputsizes);
87 end
88
89 %% take different actions depending on old and new timestep
90 if timestep == newtimestep
91 action = 'DoNothing';
92 elseif newtimestep == 0
93 action = 'MakeContinuous';
94 str = 'warning because system is sent back to continuous time';
95 utils.helper.msg(utils.const.msg.MNAME, str);
96 elseif timestep == 0
97 action = 'Discretize';
98 elseif floor(newtimestep/timestep) == newtimestep/timestep
99 action = 'TakeMultiple';
100 elseif newtimestep > timestep
101 action = 'MakeLonger';
102 else
103 action = 'MakeShorter';
104 str = 'warning because system is sent back to shorter time step';
105 utils.helper.msg(utils.const.msg.MNAME, str);
106 end
107
108 %% proceed with matrix modifications
109
110 amat_input = zeros(size(amat));
111 if isequal(action,'DoNothing')
112 %% same timestep : do nothing
113 amat_2 = amat;
114 bmat_2 = bmat;
115 elseif isequal(action,'Discretize')
116 %% from continuous to discrete
117 amat_2 = expm(amat * newtimestep);
118 amat_input = ExpInt(amat, newtimestep);
119 bmat_2 = real(amat_input * bmat);
120 if outputAntiAlias
121 cmat_2 = cmat * 0.5*(eye(size(amat_2)) + amat_2);
122 dmat_2 = dmat + cmat*0.5*bmat_2;
123 sys(i).cmats = ssm.blockMatRecut(real(cmat_2), outputsizes, sssizes);
124 sys(i).dmats = ssm.blockMatRecut(real(dmat_2), outputsizes, inputsizes);
125 end
126 elseif isequal(action,'MakeContinuous')
127 %% for discrete to continuous
128 [V, E] = eig(amat);
129 amat_2 = real(V * (diag(log(diag(E)))/timestep) * V^(-1));
130 amat_input = ExpInt(amat_2, timestep);
131 amat_input_inv = (amat_input)^(-1);
132 bmat_2 = real(amat_input_inv * bmat);
133 if outputAntiAlias
134 cmat_2 = cmat * (0.5*(eye(size(amat)) + amat))^-1;
135 dmat_2 = dmat - cmat_2*0.5*bmat;
136 sys(i).cmats = ssm.blockMatRecut(real(cmat_2), outputsizes, sssizes);
137 sys(i).dmats = ssm.blockMatRecut(real(dmat_2), outputsizes, inputsizes);
138 end
139 elseif isequal(action,'TakeMultiple')
140 %% discrete to discrete with a multiple
141 multiple = newtimestep/timestep;
142 amat_2 = amat^multiple;
143 for i_step = 1:multiple
144 amat_input = amat_input + amat^(i_step-1);
145 end
146 bmat_2 = real(amat_input * bmat);
147 if outputAntiAlias
148 cmat_0 = cmat * (0.5*(eye(size(amat)) + amat))^-1;
149 cmat_2 = cmat_0 * 0.5*(eye(size(amat_2)) + amat_2);
150 dmat_2 = dmat + cmat_0*0.5*bmat_2 - cmat_0*0.5*bmat;
151 sys(i).cmats = ssm.blockMatRecut(real(cmat_2), outputsizes, sssizes);
152 sys(i).dmats = ssm.blockMatRecut(real(dmat_2), outputsizes, inputsizes);
153 end
154 elseif isequal(action,'MakeLonger')||isequal(action,'MakeShorter')
155 %% discrete to discrete with no multiple relationship
156 [V, E] = eig(amat);
157 amat_c = real(V * (diag(log(diag(E)))/timestep) * V^(-1));
158 amat_2 = expm( amat * newtimestep/timestep);
159 amat_input_1 = ExpInt(amat_c, timestep);
160 amat_input_2 = ExpInt(amat_c, newtimestep);
161 bmat_2 = real(amat_input_2 * (amat_input_1)^(-1) * bmat);
162 if outputAntiAlias
163 cmat_0 = cmat * (0.5*(eye(size(amat)) + amat))^-1;
164 cmat_2 = cmat_0 * 0.5*(eye(size(amat_2)) + amat_2);
165 dmat_2 = dmat + cmat_0*0.5*bmat_2 - cmat_0*0.5*bmat;
166 sys(i).cmats = ssm.blockMatRecut(real(cmat_2), outputsizes, sssizes);
167 sys(i).dmats = ssm.blockMatRecut(real(dmat_2), outputsizes, inputsizes);
168 end
169 end
170 amat_2 = real(amat_2);
171 bmat_2 = real(bmat_2);
172 %% Save Matrix modifications
173 sys(i).amats = ssm.blockMatRecut(amat_2, sssizes, sssizes);
174 sys(i).bmats = ssm.blockMatRecut(bmat_2, sssizes, inputsizes);
175 sys(i).timestep = newtimestep;
176 if ~internal
177 sys(i).addHistory(getInfo('None'), pl, ssm_invars(i), sys(i).hist );
178 end
179 end
180 if nargout == numel(sys)
181 for ii = 1:numel(sys)
182 varargout{ii} = sys(ii);
183 end
184 else
185 varargout{1} = sys;
186 end
187 end
188
189 function A_int = ExpInt(A, t)
190 niter = 35;
191 A1 = A*t;
192 A2 = eye(size(A))*t;
193 A_int = zeros(size(A));
194 for i=1:niter
195 A_int = A_int + A2;
196 A2 = A1*A2/(i+1);
197 end
198 end
199
200
201 %--------------------------------------------------------------------------
202 % Get Info Object
203 %--------------------------------------------------------------------------
204 function ii = getInfo(varargin)
205
206 if nargin == 1 && strcmpi(varargin{1}, 'None')
207 sets = {};
208 pl = [];
209 else
210 sets = {'Default'};
211 pl = getDefaultPlist;
212 end
213 % Build info object
214 ii = minfo(mfilename, 'ssm', 'ltpda', utils.const.categories.op, '$Id: modifyTimeStep.m,v 1.11 2011/04/08 08:56:23 hewitson Exp $', sets, pl);
215 end
216
217 %--------------------------------------------------------------------------
218 % Get Default Plist
219 %--------------------------------------------------------------------------
220 function pl = getDefaultPlist()
221 pl = plist();
222
223 p = param({'newtimestep', 'Specify the desired new timestep.'}, paramValue.EMPTY_DOUBLE);
224 pl.append(p);
225
226 p = param({'outputAntiAlias', 'Uses a linear averaging method to compute the systems output.'}, paramValue.TRUE_FALSE);
227 pl.append(p);
228
229 end
230
231