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 8 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)
 
    ***     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 8 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 %
0033 %    ***     the plist CANNOT be arrays anymore.    ***
0034 %    Because of the very large number of output it is not desirable to have
0035 %    more than 1 simulation per function run.
0036 %
0037 %
0038 % M-FILE INFO: Get information about this methods by calling
0039 %              >> ssm.getInfo('simulate')
0040 %
0041 %              Get information about a specified set-plist by calling:
0042 %              >> ssm.getInfo('simulate', 'Default')
0043 %
0044 % ***** THERE ARE NO DEFAULT PARAMETERS *****
0045 %
0046 % VERSION: $Id:$
0047 %
0048 % HISTORY:
0049 % 04-08-2008 A Grynagier : new multiple state spaces, independant variables
0050 % in aos...
0051 % 12-06-2008 M Weyrich : change for multiple outputs, bug fixes
0052 % 17-04-2008 A Grynagier
0053 % 19-02-2008 A Grynagier
0054 % 17-02-2008 A Grynagier
0055 % Creation 30-01-2008 M Weyrich
0056 %
0057 % TO DO: Check input aos for the timestep, tsdata, and ssm.timestep
0058 % options to be defined (NL case)
0059 % add check if one input mach no ssm input variable
0060 % allow use of other LTPDA functions to generate white noise
0061 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0062 
0063 function varargout = simulate(varargin)
0064 
0065   %% starting initial checks
0066   utils.helper.msg(utils.const.msg.MNAME, ['running ', mfilename]);
0067 
0068   % Check if this is a call for parameters
0069   if utils.helper.isinfocall(varargin{:})
0070     varargout{1} = getInfo(varargin{3});
0071     return
0072   end
0073 
0074   % checking number of inputs work data
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   %% begin function body
0083 
0084   i_sys=1;
0085   % for i_sys = 1:Nsys
0086   %% retrieve system infos
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   %% finding input and variable positions of data in option plists
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       % getting position
0123       [positions, sum, value] = ssm.cellstrfind(noise_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0124       % checking against redundant input data
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       % getting position
0131       [positions, sum, value] = ssm.cellstrfind(aos_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0132       % checking against redundant input data
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       % getting position
0140       [positions, sum, value] = ssm.cellstrfind(constants_varnames,sys(i_sys).inputvarnames{i}{j},'all');
0141       % checking against redundant input data
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   %% getting correct number of samples
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     % maybe tdata should be retrieved and verified to be equal, rather than this.
0161   end
0162   if Nsamples == inf % case there is no input!
0163     Nsamples = 0;
0164   end
0165 
0166   %% generating random data
0167   [V,E]=eig(noise_in);
0168   noise_mat = (V*(E'.^0.5)*V');
0169   noise_array = noise_mat*randn(Nnoise, Nsamples);
0170   % How to generate initial position
0171   % Do one step with firts noise vector
0172   % muliply inputs by (I - Minv*A)^-1
0173   % add in the fashion    ssini = ssini.*(~isNan(ssini)) + ssini_moise.*isNaN(ssini)
0174   %
0175   %% evaluating data to return
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   %% simulation first step
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   %% simulation loop
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     %% resetting each field
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     %% evaluating differences
0215     for i=1:Nnoise %Noise  MUST BE MODIFIED IF IT IS REALLY SPECIFIED FOR TIME CONTINUOUS SYSTEM!!!!!!!
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     %% computing and storing outputs
0254     %differential equation
0255     lastX = ssm.cell_mult(Minv, ssm.cell_add(IncrementX, ssm.cell_mult(sys(i_sys).amats, lastX)));
0256     %output equations
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   %% saving in aos
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   %% construct output analysis object
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   % Build info object
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

Generated on Mon 25-Aug-2008 22:39:29 by m2html © 2003