%% Run Rapid Simulations Over Range of Parameter Values
% This example shows how to use the RSim system target file to run simulations over a range 
% of parameter values. The example uses the Van
% der Pol oscillator and performs a parameter sweep over a range of
% initial state values to obtain the phase diagram of a nonlinear system.
%
% You can customize this example for your own application by
% modifying the MATLAB(R) script used to build this example. Click the link in
% the top left corner of this page to edit the MATLAB(R) script. Click the link
% in the top right corner to run this example from MATLAB(R). Run the example 
% from writable location.  The example creates files that you may want to investigate later. 
%
% For this example, set model configuration parameter *Default parameter
% behavior* to |Inlined|. In the example, you create |Simulink.Parameter| objects as tunable
% parameters. To use RSim with *Default parameter behavior* set to
% |Tunable|, and without explicitly declaring tunable parameters, see 
% <docid:rtw_ug#mw_4b11232a-7e69-4f70-a54b-e251b63f1b75 Run Batch Simulations Without Recompiling Generated Code>.
%
% To run multiple simulations in the Simulink environment, consider using rapid accelerator
% instead of RSim. For information about rapid accelerator, see
% <docid:simulink_ug.brc689t>. To sweep parameter values, see
% <docid:simulink_ug.bu1wrbi Optimize, Estimate, and Sweep Block Parameter Values>. 
%% Preparation Model
% Make sure the current folder is writable. The example 
% creates files.

[stat, fa] = fileattrib(pwd);
if ~fa.UserWrite
    disp('This script must be run in a writable directory');
    return
end

%%
% Open the model. Then, open the *Simulink Coder* app. Set model configuration 
% parameter *System target file* to |rsim.tlc|. 
% 
% For more information on doing this graphically and setting up other RSim target
% related options,
% <matlab:helpview(fullfile(docroot,'toolbox','rtw','helptargets.map'),'config_target') look here>.

mdlName = 'rtwdemo_rsim_vdp';
open_system(mdlName);
cs = getActiveConfigSet(mdlName);
cs.switchTarget('rsim.tlc',[]);

%%
% Specify variables |INIT_X1| (the initial
% condition for state x1), |INIT_X2| (the initial condition for state |x2|),
% and |MU| (the gain value) as tunable.  To create tunable parameters, convert the variables to 
% |Simulink.Parameter| objects, and configure the objects to use a storage class 
% other than |Auto|. In this example you investigate how the
% state trajectories evolve from different initial values for states |x1| 
% and |x2| in the model.
INIT_X1 = Simulink.Parameter(INIT_X1);
INIT_X1.StorageClass = 'Model default';

INIT_X2 = Simulink.Parameter(INIT_X2);
INIT_X2.StorageClass = 'Model default';

MU = Simulink.Parameter(MU);
MU.StorageClass = 'Model default';
%%
% Define the names of files that will be created.
prmFileName = [mdlName, '_prm_sets.mat'];
logFileName = [mdlName, '_run_scr.log'];
batFileName = [mdlName, '_run_scr'];
exeFileName = mdlName;
if ispc
    exeFileName = [exeFileName, '.exe'];
    batFileName = [batFileName, '.bat'];
end
aggDataFile = [mdlName, '_results'];
startTime = cputime;


%% Build Model
% Build the RSim executable program for the model. During the build process, a
% structural checksum is calculated for the model and embedded into the
% generated executable program. This checksum is used to check that a parameter set
% passed to the executable program is compatible with the program. 

rtwbuild(mdlName);


%% Get Default Parameter Set for Model
% Get the default |rtP| structure (parameter set) for the model. The
% |modelChecksum|
% field in the |rtP| structure is the structural checksum of the model. This must
% match the checksum embedded in the RSim executable program.  If the two checksums 
% do not match, the executable program produces an
% error. |rsimgetrtp| generates an |rtP| structure with entries for the named tunable variables
% |INIT_X1|, |INIT_X2|, and |MU|.

rtp = rsimgetrtp(mdlName)


%% Create Parameter Sets
% Using the |rtp| structure, build a structure array with different
% values for the tunable variables in the model. As mentioned earlier, in this
% example you want see how the state trajectories evolve for different initial 
% values for states |x1| and |x2| in the model. To do this, generate different 
% parameter sets with different values for |INIT_X1| and |INIT_X2| and leave  
% tunable variable |MU| set to the default value.

INIT_X1_vals = -5:1:5;
INIT_X2_vals = -5:1:5;
MU_vals = 1;
nPrmSets = length(INIT_X1_vals)*length(INIT_X2_vals)*length(MU_vals)

%%
% This example has |nPrmSets| parameter sets. You need to run
% that many simulations.  Initialize |aggData|, which is a structure array
% that holds the parameter set and corresponding results.

aggData = struct('tout', [], 'yout', [], ...
                'prms', struct('INIT_X1',[],'INIT_X2',[], 'MU', []))
aggData = repmat(aggData, nPrmSets, 1);

%%
% The utility function |rsimsetrtpparam| builds the |rtP| 
% structure by adding parameter sets one at a time with different parameters
% values.

idx = 1;
for iX1 = INIT_X1_vals
    for iX2 = INIT_X2_vals
        for iMU = MU_vals
            rtp = rsimsetrtpparam(rtp,idx,'INIT_X1',iX1,'INIT_X2',iX2,'MU',iMU);
            aggData(idx).prms.INIT_X1 = iX1;
            aggData(idx).prms.INIT_X2 = iX2;
            aggData(idx).prms.MU      = iMU;
            idx = idx + 1;
        end
    end
end

%%
% Save the |rtP| structure array with the parameter sets to a MAT-file.

save(prmFileName,'rtp');


%% Create Batch File
% Create a batch script file to run the RSim executable program over the parameter
% sets. Each run reads the specified parameter set from the parameter MAT-file
% and writes the results to the specified output MAT-file.  Use the
% time out option so that if a run hangs (for example, because the model
% has a singularity for that parameter set), abort the run
% after the specified time limit and proceed to the next run.
% 
% For example, this command (on Windows(R)) specifies using the third parameter 
% set from the |rtP| structure in |prm.mat|, writes the results to |run3.mat|, 
% and aborts execution if a run takes longer than 3600 seconds of CPU time.
% While running, the script pipes messages from |model.exe| to |run.log|, which you 
% can use to debug issues.
%
%    model.exe -p prm.mat@3 -o run3.mat -L 3600 2>&1>> run.log 

fid = fopen(batFileName, 'w');
idx = 1;
for iX1 = INIT_X1_vals
    for iX2 = INIT_X2_vals
        for iMU = MU_vals
            outMatFile = [mdlName, '_run',num2str(idx),'.mat'];
            cmd = [exeFileName, ...
                   ' -p ', prmFileName, '@', int2str(idx), ...
                   ' -o ', outMatFile, ...
                   ' -L 3600'];
            if ispc
                cmd  = [cmd, ' 2>&1>> ', logFileName];
            else % (unix)
                cmd  = ['.' filesep cmd, ' 1> ', logFileName, ' 2>&1'];
            end
            fprintf(fid, ['echo "', cmd, '"\n']);
            fprintf(fid, [cmd, '\n']);
            idx = idx + 1;
        end
    end
end
if isunix,
    system(['touch ', logFileName]);
    system(['chmod +x ', batFileName]);
end
fclose(fid);

%%
% Creating a batch/script file to run simulations enables you to call the system
% command once to run the simulations (or even run the batch script outside
% MATLAB(R)) instead of calling the system command in a loop for each
% simulation. This improves performance because the
% system command has significant overhead.


%% Execute Batch File to Run Simulations
% Run the batch/script file, which runs the RSim executable program once for each
% parameter set and saves the results to a different MAT-file each time.
% You can run the batch file from outside MATLAB(R).

[stat, res] = system(['.' filesep batFileName]);
if stat ~= 0
    error(['Error running batch file ''', batFileName, ''' :', res]);
end

%%
% This example logs simulation run results in one batch file as it runs the batch
% file to sequentially run |n| simulations over |n| parameter sets. You can 
% modify the script to generate multiple batch files and run them in
% parallel by distributing them across multiple
% computers. Also, you can run the batch files without launching MATLAB(R).


%% Load Output MAT-Files and Collate Results
% Collect the simulation results from the output MAT-files into the
% structure |aggData|. If the output MAT-file corresponding to a particular run is
% not found, set the results corresponding to that run to |NaN| (not a
% number). This situation occurs if a simulation run with a particular set
% of parameters encounters singularities in the model.

idx = 1;
for iX1 = INIT_X1_vals
    for iX2 = INIT_X2_vals
        for iMU = MU_vals
            outMatFile = [mdlName, '_run',num2str(idx),'.mat'];
            if exist(outMatFile,'file')
                load(outMatFile);
                aggData(idx).tout = rt_tout;
                aggData(idx).yout = rt_yout;
            else
                aggData(idx).tout = nan;
                aggData(idx).yout = nan;
            end
            idx = idx + 1;
        end
    end
end

%%
% Save the |aggData| structure to the results MAT-file. At this point, you can
% delete the other MAT-files, as the |aggData| data structure contains the aggregation
% of input (parameters sets) and output data (simulation results).

save(aggDataFile,'aggData');
disp(['Took ', num2str(cputime-startTime), ...
      ' seconds to generate results from ', ...
      num2str(nPrmSets), ' simulation runs (Steps 2 to 7).']);

  
%% Analyze Simulation Results
% Plot the phase diagram (|X2| versus |X1|) with
% different initial values for |X1| and |X2|. The diagram shows that, irrespective
% of the initial condition, the Van der Pol oscillator converges to its natural
% oscillator mode.

colors = {'b','g','r','c','m','y','k'}; nColors = length(colors);
for idx = 1:nPrmSets
    col = colors{idx - nColors*floor(idx/nColors) + 1};
    plot(aggData(idx).prms.INIT_X1, aggData(idx).prms.INIT_X2, [col,'x'], ...
         aggData(idx).yout(:,1), aggData(idx).yout(:,2),col);
    hold on
end
grid on
xlabel('X1'); 
ylabel('X2'); 
axis('square');
title('Phase diagram for Van der Pol equation');


%% 
% Copyright 2012-2019 The MathWorks, Inc.