Home > classes > @ssm > simulate.m

simulate

PURPOSE ^

simulate simulates a discrete ssm with given inputs

SYNOPSIS ^

function varargout = simulate(varargin)

DESCRIPTION ^

simulate simulates a discrete ssm with given inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: %simulate simulates a discrete ssm with given inputs

 CALL: [ao_out] = simulate(sys, plist_inputs)

 INPUTS:
         - sys, (array of) ssm object
         - plist_inputs contains 9 fields :
            - 2 noise fields
               - 'noise variable names' with a cell vector of strings of
               the desired input variable names
               - 'covariance' which gives the covariance of this noise between
               the different corresponding inputs given for the *TIME
               CONTINUOUS* noise model
            - 2 fields for input aos :
               - 'aos variable names' with a cell vector of strings of
               the desired input variable names
               - 'mean' which gives the aos for the different
               corresponding inputs
            - 2 fields to plug in constants :
               - 'constant variable names' with a cell vector of strings of
               the desired input variable names
               - 'constants' which gives the DC value for the different
               corresponding inputs
            - 2 arrays telling which blocks should be stored and returned
               - 'return states' : for blocks of the state space (ex [8])
               - 'return outputs' : for blocks of the output (ex [1 5])
            - 1 field 'Nsamples' telling max number of samples to use
                (it is smaller if aos are shorter than Nsamples)
            - 1 field 'ssini' is a cell array of vectors that give the
               initial position for simulation. 0 is the default.
 
    ***     the plist CANNOT be arrays anymore.    ***
    Because of the very large number of output it is not desirable to have
    more than 1 simulation per function run.
 
 
 M-FILE INFO: Get information about this methods by calling
              >> ssm.getInfo('simulate')

              Get information about a specified set-plist by calling:
              >> ssm.getInfo('simulate', 'Default')

 ***** THERE ARE NO DEFAULT PARAMETERS *****

 VERSION: $Id:$

 HISTORY:
 04-08-2008 A Grynagier : new multiple state spaces, independant variables
 in aos...
 12-06-2008 M Weyrich : change for multiple outputs, bug fixes
 17-04-2008 A Grynagier
 19-02-2008 A Grynagier
 17-02-2008 A Grynagier
 Creation 30-01-2008 M Weyrich

 TO DO: Check input aos for the timestep, tsdata, and ssm.timestep
 options to be defined (NL case)
 add check if one input mach no ssm input variable
 allow use of other LTPDA functions to generate white noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 %simulate simulates a discrete ssm with given inputs
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: %simulate simulates a discrete ssm with given inputs
0005 %
0006 % CALL: [ao_out] = simulate(sys, plist_inputs)
0007 %
0008 % INPUTS:
0009 %         - sys, (array of) ssm object
0010 %         - plist_inputs contains 9 fields :
0011 %            - 2 noise fields
0012 %               - 'noise variable names' with a cell vector of strings of
0013 %               the desired input variable names
0014 %               - 'covariance' which gives the covariance of this noise between
0015 %               the different corresponding inputs given for the *TIME
0016 %               CONTINUOUS* noise model
0017 %            - 2 fields for input aos :
0018 %               - 'aos variable names' with a cell vector of strings of
0019 %               the desired input variable names
0020 %               - 'mean' which gives the aos for the different
0021 %               corresponding inputs
0022 %            - 2 fields to plug in constants :
0023 %               - 'constant variable names' with a cell vector of strings of
0024 %               the desired input variable names
0025 %               - 'constants' which gives the DC value for the different
0026 %               corresponding inputs
0027 %            - 2 arrays telling which blocks should be stored and returned
0028 %               - 'return states' : for blocks of the state space (ex [8])
0029 %               - 'return outputs' : for blocks of the output (ex [1 5])
0030 %            - 1 field 'Nsamples' telling max number of samples to use
0031 %                (it is smaller if aos are shorter than Nsamples)
0032 %            - 1 field 'ssini' is a cell array of vectors that give the
0033 %               initial position for simulation. 0 is the default.
0034 %
0035 %    ***     the plist CANNOT be arrays anymore.    ***
0036 %    Because of the very large number of output it is not desirable to have
0037 %    more than 1 simulation per function run.
0038 %
0039 %
0040 % M-FILE INFO: Get information about this methods by calling
0041 %              >> ssm.getInfo('simulate')
0042 %
0043 %              Get information about a specified set-plist by calling:
0044 %              >> ssm.getInfo('simulate', 'Default')
0045 %
0046 % ***** THERE ARE NO DEFAULT PARAMETERS *****
0047 %
0048 % VERSION: $Id:$
0049 %
0050 % HISTORY:
0051 % 04-08-2008 A Grynagier : new multiple state spaces, independant variables
0052 % in aos...
0053 % 12-06-2008 M Weyrich : change for multiple outputs, bug fixes
0054 % 17-04-2008 A Grynagier
0055 % 19-02-2008 A Grynagier
0056 % 17-02-2008 A Grynagier
0057 % Creation 30-01-2008 M Weyrich
0058 %
0059 % TO DO: Check input aos for the timestep, tsdata, and ssm.timestep
0060 % options to be defined (NL case)
0061 % add check if one input mach no ssm input variable
0062 % allow use of other LTPDA functions to generate white noise
0063 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0064 
0065 function varargout = simulate(varargin)
0066 
0067   %% starting initial checks
0068   utils.helper.msg(utils.const.msg.MNAME, ['running ', mfilename]);
0069 
0070   % Check if this is a call for parameters
0071   if utils.helper.isinfocall(varargin{:})
0072     varargout{1} = getInfo(varargin{3});
0073     return
0074   end
0075 
0076   % checking number of inputs work data
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   %% begin function body
0085 
0086   i_sys=1;
0087   % for i_sys = 1:Nsys
0088   %% retrieve system infos
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   %% finding input and variable positions of data in option plists
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       % getting position
0132       [positions, sum, value] = ssm.cellstrfind(noise_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0133       % checking against redundant input data
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       % getting position
0140       [positions, sum, value] = ssm.cellstrfind(aos_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0141       % checking against redundant input data
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       % getting position
0149       [positions, sum, value] = ssm.cellstrfind(constants_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0150       % checking against redundant input data
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   %% getting correct number of samples
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     % maybe tdata should be retrieved and verified to be equal, rather than this.
0170   end
0171   if Nsamples == inf % case there is no input!
0172     Nsamples = 0;
0173   end
0174 
0175   %% generating random data
0176   [V,E]=eig(noise_in);
0177   noise_mat = (V*(E'.^0.5)*V');
0178   noise_array = noise_mat*randn(Nnoise, Nsamples);
0179   % How to generate initial position
0180   % Do one step with firts noise vector
0181   % muliply inputs by (I - Minv*A)^-1
0182   % add in the fashion
0183   %
0184   %% evaluating data to return
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   %% simulation first step
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   %% simulation loop
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     %% resetting each field
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     %% evaluating differences
0224     for i=1:Nnoise %Noise  MUST BE MODIFIED IF IT IS REALLY SPECIFIED FOR TIME CONTINUOUS SYSTEM!!!!!!!
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     %% computing and storing outputs
0263     %differential equation
0264     lastX = ssm.cell_mult(Minv, ssm.cell_add(IncrementX, ssm.cell_mult(sys(i_sys).amats, lastX)));
0265     %output equations
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   %% saving in aos
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   %% construct output analysis object
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   % Build info object
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

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003