Mercurial > hg > ltpda
view m-toolbox/classes/@ssm/modifyTimeStep.m @ 5:5a49956df427 database-connection-manager
LTPDAPreferences panel for new LTPDADatabaseConnectionManager
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% MODIFYTIMESTEP modifies the timestep of a ssm object %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % DESCRIPTION: MODIFYTIMESTEP modifies the timestep of a ssm object, and updates % the A and B matrices supposing there is no aliasing. % % CALL: sys = modifyTimeStep(sys,pl) % % INPUTS: % sys - (array of) ssm objects % pl - A plist or numeric value giving new timestep value (param name 'newtimestep') % % OUTPUTS: % sys - (array of) ssm % % <a href="matlab:utils.helper.displayMethodInfo('ssm', 'modifyTimeStep')">Parameters Description</a> % % VERSION: $Id: modifyTimeStep.m,v 1.11 2011/04/08 08:56:23 hewitson Exp $ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function varargout = modifyTimeStep(varargin) %% Check if this is a call for parameters if utils.helper.isinfocall(varargin{:}) varargout{1} = getInfo(varargin{3}); return end %% send starting message import utils.const.* utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename); %% collecting input in_names = cell(size(varargin)); for ii = 1:nargin,in_names{ii} = inputname(ii);end [sys, ssm_invars, rest] = utils.helper.collect_objects(varargin(:), 'ssm', in_names); [pl, invars2, rest] = utils.helper.collect_objects(rest(:), 'plist'); % override with a numeric input value if ~isempty(rest) % if the user inputs a single numerical value, we use it if isnumeric(rest{1}) && numel(rest{1}) == 1 if ~isa(pl, 'plist') % perhaps we dont' have a plist yet pl = plist(); end pl.pset('newtimestep', rest{1}); % set newtimestep else pl = combine(pl, plist(rest{:})); end end pl = combine(pl, getDefaultPlist()); %%% Internal call: Only one object + don't look for a plist internal = strcmp(varargin{end}, 'internal'); %% dealing with user custom inputs Nsys = numel(sys); %% begin function body %% building output depending on options sys = copy(sys,nargout); newtimestep = find(pl,'newtimestep'); outputAntiAlias = find(pl, 'outputAntiAlias'); if numel(newtimestep)~=1 error('parameter "newtimestep" should be of size 1x1') end for i = 1:Nsys if ~sys(i).isnumerical error(['error in modifTimeStep because system "',sys(i).name, '" should be numerical to be modified' ]); end %% retrieving input data timestep = sys(i).timestep; sssizes = sys(i).sssizes; inputsizes = sys(i).inputsizes; amat = ssm.blockMatFusion(sys(i).amats, sssizes, sssizes); bmat = ssm.blockMatFusion(sys(i).bmats, sssizes, inputsizes); if outputAntiAlias outputsizes = sys(i).outputsizes; cmat = ssm.blockMatFusion(sys(i).cmats, outputsizes, sssizes); dmat = ssm.blockMatFusion(sys(i).dmats, outputsizes, inputsizes); end %% take different actions depending on old and new timestep if timestep == newtimestep action = 'DoNothing'; elseif newtimestep == 0 action = 'MakeContinuous'; str = 'warning because system is sent back to continuous time'; utils.helper.msg(utils.const.msg.MNAME, str); elseif timestep == 0 action = 'Discretize'; elseif floor(newtimestep/timestep) == newtimestep/timestep action = 'TakeMultiple'; elseif newtimestep > timestep action = 'MakeLonger'; else action = 'MakeShorter'; str = 'warning because system is sent back to shorter time step'; utils.helper.msg(utils.const.msg.MNAME, str); end %% proceed with matrix modifications amat_input = zeros(size(amat)); if isequal(action,'DoNothing') %% same timestep : do nothing amat_2 = amat; bmat_2 = bmat; elseif isequal(action,'Discretize') %% from continuous to discrete amat_2 = expm(amat * newtimestep); amat_input = ExpInt(amat, newtimestep); bmat_2 = real(amat_input * bmat); if outputAntiAlias cmat_2 = cmat * 0.5*(eye(size(amat_2)) + amat_2); dmat_2 = dmat + cmat*0.5*bmat_2; sys(i).cmats = ssm.blockMatRecut(real(cmat_2), outputsizes, sssizes); sys(i).dmats = ssm.blockMatRecut(real(dmat_2), outputsizes, inputsizes); end elseif isequal(action,'MakeContinuous') %% for discrete to continuous [V, E] = eig(amat); amat_2 = real(V * (diag(log(diag(E)))/timestep) * V^(-1)); amat_input = ExpInt(amat_2, timestep); amat_input_inv = (amat_input)^(-1); bmat_2 = real(amat_input_inv * bmat); if outputAntiAlias cmat_2 = cmat * (0.5*(eye(size(amat)) + amat))^-1; dmat_2 = dmat - cmat_2*0.5*bmat; sys(i).cmats = ssm.blockMatRecut(real(cmat_2), outputsizes, sssizes); sys(i).dmats = ssm.blockMatRecut(real(dmat_2), outputsizes, inputsizes); end elseif isequal(action,'TakeMultiple') %% discrete to discrete with a multiple multiple = newtimestep/timestep; amat_2 = amat^multiple; for i_step = 1:multiple amat_input = amat_input + amat^(i_step-1); end bmat_2 = real(amat_input * bmat); if outputAntiAlias cmat_0 = cmat * (0.5*(eye(size(amat)) + amat))^-1; cmat_2 = cmat_0 * 0.5*(eye(size(amat_2)) + amat_2); dmat_2 = dmat + cmat_0*0.5*bmat_2 - cmat_0*0.5*bmat; sys(i).cmats = ssm.blockMatRecut(real(cmat_2), outputsizes, sssizes); sys(i).dmats = ssm.blockMatRecut(real(dmat_2), outputsizes, inputsizes); end elseif isequal(action,'MakeLonger')||isequal(action,'MakeShorter') %% discrete to discrete with no multiple relationship [V, E] = eig(amat); amat_c = real(V * (diag(log(diag(E)))/timestep) * V^(-1)); amat_2 = expm( amat * newtimestep/timestep); amat_input_1 = ExpInt(amat_c, timestep); amat_input_2 = ExpInt(amat_c, newtimestep); bmat_2 = real(amat_input_2 * (amat_input_1)^(-1) * bmat); if outputAntiAlias cmat_0 = cmat * (0.5*(eye(size(amat)) + amat))^-1; cmat_2 = cmat_0 * 0.5*(eye(size(amat_2)) + amat_2); dmat_2 = dmat + cmat_0*0.5*bmat_2 - cmat_0*0.5*bmat; sys(i).cmats = ssm.blockMatRecut(real(cmat_2), outputsizes, sssizes); sys(i).dmats = ssm.blockMatRecut(real(dmat_2), outputsizes, inputsizes); end end amat_2 = real(amat_2); bmat_2 = real(bmat_2); %% Save Matrix modifications sys(i).amats = ssm.blockMatRecut(amat_2, sssizes, sssizes); sys(i).bmats = ssm.blockMatRecut(bmat_2, sssizes, inputsizes); sys(i).timestep = newtimestep; if ~internal sys(i).addHistory(getInfo('None'), pl, ssm_invars(i), sys(i).hist ); end end if nargout == numel(sys) for ii = 1:numel(sys) varargout{ii} = sys(ii); end else varargout{1} = sys; end end function A_int = ExpInt(A, t) niter = 35; A1 = A*t; A2 = eye(size(A))*t; A_int = zeros(size(A)); for i=1:niter A_int = A_int + A2; A2 = A1*A2/(i+1); end end %-------------------------------------------------------------------------- % Get Info Object %-------------------------------------------------------------------------- function ii = getInfo(varargin) if nargin == 1 && strcmpi(varargin{1}, 'None') sets = {}; pl = []; else sets = {'Default'}; pl = getDefaultPlist; end % Build info object 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); end %-------------------------------------------------------------------------- % Get Default Plist %-------------------------------------------------------------------------- function pl = getDefaultPlist() pl = plist(); p = param({'newtimestep', 'Specify the desired new timestep.'}, paramValue.EMPTY_DOUBLE); pl.append(p); p = param({'outputAntiAlias', 'Uses a linear averaging method to compute the systems output.'}, paramValue.TRUE_FALSE); pl.append(p); end