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
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063 function varargout = simulate(varargin)
0064
0065
0066 utils.helper.msg(utils.const.msg.MNAME, ['running ', mfilename]);
0067
0068
0069 if utils.helper.isinfocall(varargin{:})
0070 varargout{1} = getInfo(varargin{3});
0071 return
0072 end
0073
0074
0075 if ~nargin==2, error('wrong number of inputs!'), end
0076
0077 if ~isequal(class(varargin{1}),'ssm'), error(['argument is not a ssm but a ', class(varargin{1})]),end
0078 if ~isequal(class(varargin{2}),'plist'), error(['argument is not a plist but a ', class(varargin{2})]),end
0079
0080 sys = varargin{1};
0081 plist_inputs = combine(varargin{2}, getDefaultPlist('Default'));
0082
0083
0084 i_sys=1;
0085
0086
0087 if ~sys(i_sys).isnumerical
0088 error(['error because system ',sys(i_sys).name,' is not numerical']);
0089 end
0090 M2 = ssm.cell_fusion(sys(i_sys).mmats, sys(i_sys).sssizes, sys(i_sys).sssizes);
0091 Minv = ssm.cell_recut(M2^-1, sys(i_sys).sssizes, sys(i_sys).sssizes);
0092 if ~isequal(M2, eye(size(M2)))
0093 utils.helper.msg(utils.const.msg.MNAME, 'warning, M is not the identity!!! It is:');
0094 end
0095 sssizes = sys(i_sys).sssizes;
0096 ssini = sys(i_sys).ssini;
0097 timestep = sys(i_sys).timestep;
0098 Ninputs = sys(i_sys).Ninputs;
0099 Noutputs = sys(i_sys).Noutputs;
0100 Nss = sys(i_sys).Nss;
0101 outputsizes = sys(i_sys).outputsizes;
0102 inputsizes = sys(i_sys).inputsizes;
0103
0104
0105 noise_varnames = find(plist_inputs, 'noise variable names');
0106 Nnoise = length(noise_varnames);
0107 noise_in = find(plist_inputs, 'variance');
0108 noise_pos = zeros(Nnoise,2);
0109
0110 aos_varnames = find(plist_inputs, 'aos variable names');
0111 Naos_in = length(aos_varnames);
0112 aos_in = find(plist_inputs, 'aos');
0113 aos_pos = zeros(size(aos_varnames),2);
0114
0115 constants_varnames = find(plist_inputs, 'constant variable names');
0116 Nconstants = length(constants_varnames);
0117 constants_in = find(plist_inputs, 'constants');
0118 constants_pos = zeros(size(constants_varnames),2);
0119
0120 for i=1:Ninputs
0121 for j=1:inputsizes(i)
0122
0123 [positions, sum, value] = ssm.cellstrfind(noise_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0124
0125 if sum==1
0126 noise_pos(positions,:) = [i,j];
0127 elseif sum>1
0128 error(['more than one noise variable is named ', value]);
0129 end
0130
0131 [positions, sum, value] = ssm.cellstrfind(aos_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0132
0133 if sum==1
0134 aos_pos(positions,:) = [i,j];
0135 aos_data{positions} = aos_in(positions).data.getY ;
0136 elseif sum>1
0137 error(['more than one so variable is named ', value]);
0138 end
0139
0140 [positions, sum, value] = ssm.cellstrfind(constants_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0141
0142 if sum==1
0143 constants_pos(positions,:) = [i,j];
0144 elseif sum>1
0145 error(['more than one constant variable is named ', value]);
0146 end
0147 end
0148 end
0149
0150 Nsamples = find(plist_inputs, 'Nsamples');
0151 f0 = 1/timestep;
0152 for i=1:Naos_in
0153 Nsamples = min(Nsamples,length(aos_in(i).data.y));
0154 if ~(f0==aos_in(i).data.fs)
0155 str = ['WARNING : ssm frequency is ',num2str(f0),...
0156 ' but sampling frequency of ao named ',...
0157 aos_in(i).name, ' is ', num2str(aos_in(i).data.fs) ];
0158 utils.helper.msg(utils.const.msg.MNAME, str);
0159 end
0160
0161 end
0162 if Nsamples == inf
0163 Nsamples = 0;
0164 end
0165
0166
0167 [V,E]=eig(noise_in);
0168 noise_mat = (V*(E'.^0.5)*V');
0169 noise_array = noise_mat*randn(Nnoise, Nsamples);
0170
0171
0172
0173
0174
0175
0176 return_states = find(plist_inputs, 'return states');
0177 returnStates = false(1,Nss);
0178 for i=1:Nss
0179 returnStates(i) = ~isempty(find( ismember(return_states, i)));
0180 end
0181 return_outputs = find(plist_inputs, 'return outputs');
0182 returnOutputs = false(1, Noutputs);
0183 for i=1:Noutputs
0184 returnOutputs(i) = ~isempty(find( ismember(return_outputs, i)));
0185 end
0186
0187 x = cell(Nss, 1);
0188 lastX = ssini;
0189 for i=1:Nss
0190 if returnStates(i)
0191 x{i} = zeros(sssizes(i), Nsamples);
0192 end
0193 end
0194 y = cell(Noutputs, 1);
0195 for i=1:Noutputs
0196 if returnOutputs(i)
0197 y{i} = zeros(outputsizes(i), Nsamples);
0198 end
0199 end
0200 IncrementX = cell(Nss,1);
0201 IncrementY = cell(Noutputs,1);
0202
0203 for k = 1:Nsamples
0204 if rem(k,500)==0
0205 utils.helper.msg(utils.const.msg.MNAME, ['simulation time : ',num2str(k*timestep) ]);
0206 end
0207
0208 for j = 1:Nss
0209 IncrementX{j} = zeros(sssizes(j),1);
0210 end
0211 for j = 1:Noutputs
0212 IncrementY{j} = zeros(outputsizes(j),1);
0213 end
0214
0215 for i=1:Nnoise
0216 for s=1:Nss
0217 if ~ isempty(sys(i_sys).bmats{s,noise_pos(i,1)})
0218 IncrementX{s} = IncrementX{s} + sys(i_sys).bmats{s,noise_pos(i,1)}(:,noise_pos(i,2))*noise_array(i,k);
0219 end
0220 end
0221 for j = 1:Noutputs
0222 if ~ isempty(sys(i_sys).dmats{j,noise_pos(i,1)})
0223 IncrementY{j} = IncrementY{j} + sys(i_sys).dmats{j,noise_pos(i,1)}(:,noise_pos(i,2))*noise_array(i,k);
0224 end
0225 end
0226 end
0227
0228 for i=1:Naos_in
0229 for s=1:Nss
0230 if ~ isempty(sys(i_sys).bmats{s,aos_pos(i,1)})
0231 IncrementX{s} = IncrementX{s} + sys(i_sys).bmats{s,aos_pos(i,1)}(:,aos_pos(i,2))*aos_data{i}(k);
0232 end
0233 end
0234 for j = 1:Noutputs
0235 if ~ isempty(sys(i_sys).dmats{j,aos_pos(i,1)})
0236 IncrementY{j} = IncrementY{j} + sys(i_sys).dmats{j,aos_pos(i,1)}(:,aos_pos(i,2))*aos_data{i}(k);
0237 end
0238 end
0239 end
0240
0241 for i=1:Nconstants
0242 for s=1:Nss
0243 if ~ isempty(sys(i_sys).bmats{s,constants_pos(i,1)})
0244 IncrementX{s} = IncrementX{s} + sys(i_sys).bmats{s,constants_pos(i,1)}(:,constants_pos(i,2))*constants_in(i);
0245 end
0246 end
0247 for j = 1:Noutputs
0248 if ~ isempty(sys(i_sys).dmats{j,constants_pos(i,1)})
0249 IncrementY{j} = IncrementY{j} + sys(i_sys).dmats{j,constants_pos(i,1)}(:,constants_pos(i,2))*constants_in(i);
0250 end
0251 end
0252 end
0253
0254
0255 lastX = ssm.cell_mult(Minv, ssm.cell_add(IncrementX, ssm.cell_mult(sys(i_sys).amats, lastX)));
0256
0257 lastY = ssm.cell_add(ssm.cell_mult(sys(i_sys).cmats, lastX) ,IncrementY);
0258 for s=1:Nss
0259 if returnStates(s)
0260 x{s}(:,k) = lastX{s};
0261 end
0262 end
0263 for o=1:Noutputs
0264 if returnOutputs(o)
0265 y{o}(:,k) = lastY{o};
0266 end
0267 end
0268 end
0269
0270 Naos = 0;
0271 for i=1:Nss
0272 if returnStates(i)
0273 Naos = Naos +sssizes(i);
0274 end
0275 end
0276 for i=1:Noutputs
0277 if returnOutputs(i)
0278 Naos = Naos +outputsizes(i);
0279 end
0280 end
0281 if Naos>0
0282 ao_out(Naos) = ao;
0283 if ~exist('t','var')
0284 t = (1:Nsamples)*timestep;
0285 end
0286 count = 1;
0287 for i=1:Nss
0288 if returnStates(i)
0289 for j=1:sssizes(i)
0290 ao_out(count) = ao(tsdata(t,x{i}(j,:)));
0291 ao_out(count).setName(['sim_X_',num2str(i),'_',num2str(j)]);
0292 ao_out(count).setDescription(...
0293 ['ssm#: ',num2str(i_sys),...
0294 ', state: ',sys(i_sys).ssnames{i},' ',sys(i_sys).ssvarnames{i}{j}]);
0295 ao_out(count) = ao_out(count).addHistory(ssm.getInfo(mfilename), plist_inputs , {''}, sys(i_sys).hist );
0296 count = count+1;
0297 end
0298 end
0299 end
0300 for i=1:Noutputs
0301 if returnOutputs(i)
0302 for j=1:outputsizes(i)
0303 ao_out(count) = ao(tsdata(t,y{i}(j,:)));
0304 ao_out(count).setName(['sim_Y_',num2str(i),'_',num2str(j)]);
0305 ao_out(count).setDescription(...
0306 ['ssm: ',num2str(i_sys),...
0307 ', output: ',sys(i_sys).outputnames{i},' variable: ',sys(i_sys).outputvarnames{i}{j}]);
0308 ao_out(count) = ao_out(count).addHistory(ssm.getInfo(mfilename), plist_inputs , {''}, sys(i_sys).hist );
0309 count = count+1;
0310 end
0311 end
0312 end
0313 end
0314
0315 varargout = {ao_out};
0316 end
0317
0318 function ii = getInfo(varargin)
0319 if nargin == 1 && strcmpi(varargin{1}, 'None')
0320 sets = {};
0321 pls = [];
0322 elseif nargin == 1 && ~isempty(varargin{1}) && ischar(varargin{1})
0323 sets{1} = varargin{1};
0324 pls = getDefaultPlist(sets{1});
0325 else
0326 sets = {'Default'};
0327 pls = [];
0328 for kk=1:numel(sets)
0329 pls = [pls getDefaultPlist(sets{kk})];
0330 end
0331 end
0332
0333 ii = minfo(mfilename, 'ssm', '', 'STATESPACE', '$Id: resp.m,v 1.17 2008/07/22 10:22:38 ingo Exp $', sets, pls);
0334 end
0335
0336 function plo = getDefaultPlist(set)
0337 switch set
0338 case 'Default'
0339 plo = plist(...
0340 'noise variable names', {}, ...
0341 'variance', [], ...
0342 'aos variable names', {}, ...
0343 'aos', [], ...
0344 'constant variable names',{}, ...
0345 'constant', [], ...
0346 'return states', [], ...
0347 'return outputs', 1, ...
0348 'Nsamples', 1e4 ...
0349 );
0350 otherwise
0351 error('### Unknown parameter set [%s].', set);
0352 end
0353 end