Home > m > timetools > statespacefunctions > ltpda_ss_simulate.m

ltpda_ss_simulate

PURPOSE ^

simulates a discrete state space system

SYNOPSIS ^

function varargout = ltpda_ss_simulate(varargin)

DESCRIPTION ^

simulates a discrete state space system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: simulates a discrete state space system

 CALL: ao_out = simulate_discrete(ao_in, syst)

 INPUTS: ao_in_list - (array of) plists of :
                         - ao containing concatenated data stream for each
                         input
                         - plist describing an input noise with a mean 
                         'MEAN' (vector) and a sigma 'SIGMA' (Matrix)
 syst - (array of) system data in plist format

 OUTPUTS: ao_out - (array of) output analysis object in two dimensions
 (first for the AOs, second for the systems) 

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

 &&VERSION: $Id: ltpda_ss_simulate.m,v 1.1 2008/03/11 16:52:56 adrien Exp $

 HISTORY: 19-02-2008 A Grynagier
 17-02-2008 A Grynagier
 Creation 30-01-2008 M Weyrich

 TO DO: options to be defined (NL case)
 use LTPDA functions to generate white noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function varargout = ltpda_ss_simulate(varargin)
0002 %simulates a discrete state space system
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % DESCRIPTION: simulates a discrete state space system
0006 %
0007 % CALL: ao_out = simulate_discrete(ao_in, syst)
0008 %
0009 % INPUTS: ao_in_list - (array of) plists of :
0010 %                         - ao containing concatenated data stream for each
0011 %                         input
0012 %                         - plist describing an input noise with a mean
0013 %                         'MEAN' (vector) and a sigma 'SIGMA' (Matrix)
0014 % syst - (array of) system data in plist format
0015 %
0016 % OUTPUTS: ao_out - (array of) output analysis object in two dimensions
0017 % (first for the AOs, second for the systems)
0018 %
0019 % ***** THERE ARE NO DEFAULT PARAMETERS *****
0020 %
0021 % &&VERSION: $Id: ltpda_ss_simulate.m,v 1.1 2008/03/11 16:52:56 adrien Exp $
0022 %
0023 % HISTORY: 19-02-2008 A Grynagier
0024 % 17-02-2008 A Grynagier
0025 % Creation 30-01-2008 M Weyrich
0026 %
0027 % TO DO: options to be defined (NL case)
0028 % use LTPDA functions to generate white noise
0029 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0030 
0031 ALGONAME = mfilename;
0032 VERSION = '$Id: ltpda_ss_simulate.m,v 1.1 2008/03/11 16:52:56 adrien Exp $';
0033 CATEGORY = 'STATESPACE';
0034 display(['starting ' ALGONAME]);
0035 
0036 if not(isempty(varargin))
0037     if isequal( varargin{1}, 'Version')
0038         varargout = VERSION;
0039         return;
0040     elseif isequal(varargin{1}, 'Params')
0041         varargout = plist();
0042         return;
0043     elseif isequal(varargin{1}, 'Category')
0044         varargout = CATEGORY;
0045         return;
0046     end
0047 end
0048 systs = varargin{2};
0049 ao_in_lists = varargin{1};
0050 ao_out = ao();
0051 
0052 for i_syst = 1:length(systs)
0053     for i_ao_in = 1:length(ao_in_lists)
0054         syst =systs(i_syst);
0055         ao_in_list = ao_in_lists(i_ao_in);
0056 
0057         %% retrieve system infos
0058         A = find(syst, 'AMAT');
0059         B = find(syst, 'BMATS');
0060         C = find(syst, 'CMAT');
0061         D = find(syst, 'DMATS');
0062         A = A{1};
0063         C = C{1};
0064         nss = size(A,1);
0065         %check iniside A, Bs, c, Ds there is no symbolic variable
0066         Xini = find(syst, 'XINI');
0067         timestep = find(syst, 'TIMESTEP');
0068         NBINPUTS = find(syst, 'NBINPUTS');
0069         INPUTISUSED = find(syst, 'INPUTISUSED');
0070         INPUTNAMES = find(syst, 'INPUTNAMES');
0071         INPUTSIZES = find(syst, 'INPUTSIZES');
0072         XISOUTPUT = find(syst, 'XISOUTPUT');
0073         YISOUTPUT = find(syst, 'YISOUTPUT');
0074 
0075 
0076 
0077         %% read out AO (time-series data)
0078         ao_in = cell(NBINPUTS,1);
0079         for i=1:NBINPUTS
0080             if INPUTISUSED(i)
0081                 ao_in{i} = find(ao_in_list, INPUTNAMES{i});
0082                 if isequal(class(ao_in{i}), 'ao')
0083                     if not(isequal(ao_in{i}.data.x,{})) %making sure ao is not empty
0084                         t = ao_in{i}.data.x;
0085                         dt = diff(t);
0086                         if sum(not(isequal(dt,timestep)))>1
0087                             msg = ['error in input number ', i , ' called ', INPUTNAMES{i}, 'because wrong timestep for simulation'];
0088                             error(msg);
0089                         end
0090                         if not( isequal( size(B{i},1), size(D{i},1), size(ao_in{i}.data.y,1) ) )
0091                             msg = ['error in input number ', i , ' called ', INPUTNAMES{i}, 'because wrong input size'];
0092                             error(msg);
0093                         end
0094                     end
0095                 end
0096             end
0097         end
0098         %% getting correct number of samples
0099         nsamples = inf;
0100         for i=1:NBINPUTS
0101             if isequal(class(ao_in{i}), 'ao')
0102                 nsamples = min(nsamples,size(ao_in{i}.data.y,1));
0103             end
0104         end
0105         if nsamples == inf % case there is no input!
0106             nsamples = 0;
0107         end
0108 
0109         %% generating random data
0110         for i=1:NBINPUTS
0111             if isequal(class(ao_in{i}), 'plist')
0112                 sigma = find(ao_in{i}, 'SIGMA');
0113                 mean = find(ao_in{i}, 'MEAN');
0114                 if sigma==0
0115                     ao_in{i} = ao(ones(nsamples,INPUTSIZES(i))*diag(mean));
0116                 else
0117                     ao_in{i} = ao(ones(nsamples,INPUTSIZES(i))*diag(mean) + randn(nsamples,INPUTSIZES(i))*sigma);
0118                 end
0119             end
0120         end
0121 
0122         %% simulation first step
0123         x = zeros(nss, nsamples);
0124         y = zeros(size(C,1), nsamples);
0125         x(:,1) = Xini;
0126         IncrementY = zeros(size(C,1),1);
0127         for i=1:NBINPUTS
0128             if INPUTISUSED(i)
0129                 IncrementY = IncrementY +  D{i}* ao_in{i}.data.y(1,:)' ;
0130             end
0131         end
0132         y(:,1) = C* x(:,1) +IncrementY;
0133 
0134         %% simulation loop
0135         for k = 2:nsamples
0136             IncrementX = zeros(nss,1);
0137             IncrementY = zeros(size(C,1),1);
0138             for i=1:NBINPUTS
0139                 if INPUTISUSED(i)
0140                     IncrementX = IncrementX +  B{i}*ao_in{i}.data.y(k,:)' ;
0141                     IncrementY = IncrementY +  D{i}*ao_in{i}.data.y(k,:)' ;
0142                 end
0143             end
0144             x(:,k) = A* x(:,k-1) + IncrementX ;
0145             y(:,k) = C* x(:,k) + IncrementY ;
0146         end
0147 
0148         %% data object
0149         out_dat = zeros(0, nsamples);
0150         if XISOUTPUT
0151             out_dat = [out_dat; x]; %#ok<AGROW>
0152         end
0153         if YISOUTPUT
0154             out_dat = [out_dat; y]; %#ok<AGROW>
0155         end
0156 
0157         if exist('t','var')
0158             ts = tsdata(t,out_dat) ;
0159         else
0160             ts = tsdata(1:nsamples,out_dat) ;
0161         end
0162 
0163         %% history object, old history from ao_in plus new
0164         % hi_old = ao_in.hist ;
0165         hi_new = history('simulate_discrete','version', syst);
0166         % hi = [hi_old hi_new];
0167         % h = history(hi);
0168         h = history(hi_new);
0169         ao_out_loc = ao(ts,h);
0170         ao_out_loc.name = ALGONAME;
0171         ao_out(i_syst , i_ao_in) = ao_out_loc;
0172     end
0173 end
0174 %% construct output analysis object
0175 varargout  = {ao_out} ;
0176 end

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003