0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036 function varargout = modiftimestep(varargin)
0037
0038
0039 utils.helper.msg(utils.const.msg.MNAME, ['running ', mfilename]);
0040
0041
0042 if utils.helper.isinfocall(varargin{:})
0043 varargout{1} = getInfo(varargin{3});
0044 return
0045 end
0046
0047
0048 if ~nargin==2, error('wrong number of inputs!'), end
0049
0050 if ~isequal(class(varargin{1}),'ssm'), error(['argument is not a ssm but a ', class(varargin{1})]),end
0051 if ~isequal(class(varargin{2}),'plist'), error(['argument is not a plist but a ', class(varargin{2})]),end
0052
0053 sys_in = varargin{1};
0054 options = combine(varargin{2}, getDefaultPlist('Default'));
0055 Nsys = length(sys_in);
0056 Noptions = length(options);
0057
0058
0059
0060 if length(options)==1
0061 sys_out(:,1) = reshape(sys_in,length(sys_in),1);
0062 else
0063 for j=1:Noptions
0064 sys_out(:,j)=copy(reshape(sys_in,length(sys_in),1),true);
0065 end
0066 end
0067
0068 if ~(Noptions==1)
0069 sys_out(Nsys,Noptions) = ssm();
0070 end
0071
0072 for i = 1:Nsys
0073 for j = 1:Noptions
0074
0075 options(j) = combine(options(j), getDefaultPlist('Default'));
0076 newtimestep = find(options(j),'newtimestep');
0077 timestep = sys_in(i).timestep;
0078 sssizes = sys_in(i).sssizes;
0079 inputsizes = sys_in(i).inputsizes;
0080 amat = ssm.cell_fusion(sys_in(i).amats, sssizes, sssizes);
0081 mmat = ssm.cell_fusion(sys_in(i).mmats, sssizes, sssizes);
0082 bmat = ssm.cell_fusion(sys_in(i).bmats, sssizes, inputsizes);
0083
0084 minv = mmat^(-1);
0085
0086 amat = minv*amat;
0087 bmat = minv*bmat;
0088 mmat_2 = eye(size(amat));
0089
0090 if timestep == newtimestep
0091 action = 'DoNothing';
0092 elseif newtimestep == 0
0093 action = 'MakeContinuous';
0094 str = 'warning because system is sent back to continuous time';
0095 utils.helper.msg(utils.const.msg.MNAME, str);
0096 elseif timestep == 0
0097 action = 'Discretize';
0098 elseif floor(newtimestep/timestep) == newtimestep/timestep
0099 action = 'TakeMultiple';
0100 elseif newtimestep > timestep
0101 action = 'MakeLonger';
0102 else
0103 action = 'MakeShorter';
0104 str = 'warning because system is sent back to shorter time step';
0105 utils.helper.msg(utils.const.msg.MNAME, str);
0106 end
0107
0108
0109 NSubstep = 35;
0110 amat_input = zeros(size(amat));
0111 if isequal(action,'DoNothing')
0112
0113 amat_2 = amat;
0114 bmat_2 = bmat;
0115 elseif isequal(action,'Discretize')
0116
0117 amat_2 = expm(amat * newtimestep);
0118 for i_step = 1:NSubstep
0119 position = newtimestep*(2*i_step-1)/(2*NSubstep);
0120 amat_input = amat_input + expm(amat*position);
0121 end
0122 amat_input = amat_input*newtimestep/NSubstep;
0123 bmat_2 = real(amat_input * bmat);
0124 elseif isequal(action,'MakeContinuous')
0125
0126 [V, E] = eig(amat);
0127 amat_2 = V * (diag(diag(log(E)))/timestep) * V^(-1);
0128 for i_step = 1:NSubstep
0129 position = timestep*(2*i_step-1)/(2*NSubstep);
0130 amat_input = amat_input + expm( amat_2*position );
0131 end
0132 amat_input = amat_input*timestep/NSubstep;
0133 amat_input_inv = (amat_input)^(-1);
0134 bmat_2 = real(amat_input_inv * bmat);
0135 elseif isequal(action,'TakeMultiple')
0136
0137 multiple = newtimestep/timestep;
0138 amat_2 = amat^multiple;
0139 for i_step = 1:multiple
0140 amat_input = amat_input + amat^(i_step-1);
0141 end
0142 bmat_2 = real(amat_input * bmat);
0143
0144 elseif isequal(action,'MakeLonger')||isequal(action,'MakeShorter')
0145
0146 [V, E] = eig(amat);
0147 amat_1 = V * (diag(diag(log(E)))/timestep) * V^(-1);
0148 amat_2 = expm( amat_1 * newtimestep);
0149 amat_input_1 = amat_input;
0150 amat_input_2 = amat_input;
0151 for i_step = 1:NSubstep
0152 position_1 = timestep*(2*i_step-1)/(2*NSubstep);
0153 position_2 = newtimestep*(2*i_step-1)/(2*NSubstep);
0154 amat_input_1 = amat_input_1 + expm( amat_1*position_1 );
0155 amat_input_2 = amat_input_2 + expm( amat_1*position_2 );
0156 end
0157 amat_input_1 = amat_input_1*timestep/NSubstep;
0158 amat_input_2 = amat_input_2*newtimestep/NSubstep;
0159 amat_input_1inv = (amat_input_1)^(-1);
0160 bmat_2 = real(amat_input_2 * amat_input_1inv * bmat);
0161 end
0162 amat_2 = real(amat_2);
0163 bmat_2 = real(bmat_2);
0164
0165 sys_out(i,j).amats = ssm.cell_recut(amat_2, sssizes, sssizes);
0166 sys_out(i,j).mmats = ssm.cell_recut(mmat_2, sssizes, sssizes);
0167 sys_out(i,j).bmats = ssm.cell_recut(bmat_2, sssizes, inputsizes);
0168 sys_out(i,j).timestep = newtimestep;
0169 sys_out(i,j).addHistory(ssm.getInfo(mfilename), options(j) , {''}, sys_in(i).hist );
0170 end
0171 end
0172 if ~(Noptions==1)
0173 varargout = {sys_out};
0174 else
0175 varargout = {sys_in};
0176 end
0177 end
0178
0179
0180 function ii = getInfo(varargin)
0181 if nargin == 1 && strcmpi(varargin{1}, 'None')
0182 sets = {};
0183 pls = [];
0184 elseif nargin == 1 && ~isempty(varargin{1}) && ischar(varargin{1})
0185 sets{1} = varargin{1};
0186 pls = getDefaultPlist(sets{1});
0187 else
0188 sets = {'Default'};
0189 pls = [];
0190 for kk=1:numel(sets)
0191 pls = [pls getDefaultPlist(sets{kk})];
0192 end
0193 end
0194
0195 ii = minfo(mfilename, 'ssm', '', 'STATESPACE', '$Id: resp.m,v 1.17 2008/07/22 10:22:38 ingo Exp $', sets, pls);
0196 end
0197
0198 function plo = getDefaultPlist(set)
0199 switch set
0200 case 'Default'
0201 plo = plist('newtimestep',0.1);
0202 otherwise
0203 error('### Unknown parameter set [%s].', set);
0204 end
0205 end
0206