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