MonteCarlo.m
From charlesreid1
Code
The following is the Matlab code for my Monte Carlo driver.
M file is available here: http://files.charlesmartinreid.com/ExperimentalDesign/MonteCarlo.m
% This is the Monte Carlo driver
% for response surface creation
keep_existing_data=true;
%% Using Matlab Parallel Toolbox:
%matlabpool open 8
%% (pools larger than 8 aren't supported by Matlab as of 2010a)
% Active parameters:
% (see http://wiki.charlesmartinreid.com:8888/wiki/Monte_Carlo_Experimental_Design)
%
% I U
% ---------
% m_dot 1.5 +/- 0.5
% k(T) 1.5 * 10^{ 0 +/- 2} (Log scale)
% L_mix 0.3 - 3.0 (Log scale)
% x_i 0.5 +/- 0.02
% 1.5 +/- 0.02
% 2.5 +/- 0.02
N_samples = 10000;
N_variables = 6;
% Set variable names
names{1} = 'Mass flowrate';
names{2} = 'Rxn rate k(T)';
names{3} = 'L_{mix}';
names{4} = 'z_1';
names{5} = 'z_2';
names{6} = 'z_3';
% Set minimum for variable ranges
alpha = zeros(N_variables,1);
alpha(1) = 0.0; % mdot (this is set later...)
alpha(2) = 1.0e-2; % k(T)
alpha(3) = 0.3; % Lmix
alpha(4) = 0.48; % z1 = 0.5 m
alpha(5) = 1.48; % z2 = 1.5 m
alpha(6) = 2.48; % z3 = 2.5 m
% Set maximum for variable ranges
beta = zeros(N_variables,1);
beta(1) = 0.0; % mdot (this is set later...)
beta(2) = 1.0e2; % k(T)
beta(3) = 3.0; % Lmix
beta(4) = 0.52; % z1 = 0.5 m
beta(5) = 1.52; % z2 = 1.5 m
beta(6) = 2.52; % z3 = 2.5 m
% Set log scale or not
logscale(1) = false; % mdot
logscale(2) = true; % k(T)
logscale(3) = true; % Lmix
logscale(4) = false; % z1
logscale(5) = false; % z2
logscale(6) = false; % z3
% Samples = actual parameter values for each of the N_samples random samples
% SamplesNormalized = random number values for each of the N_samples random samples
Samples = zeros(N_variables, N_samples);
SamplesNormalized = zeros(N_variables, N_samples);
% Samples for ONLY mdot = 1 response surface
Samples_mdot1 = zeros(N_variables, 1);
% Samples for ONLY mdot = 2 response surface
Samples_mdot2 = zeros(N_variables, 1);
% (both of these variables will expand to the proper size...)
% Create arrays for storing function outputs
% 3 = number of species concentrations, [A] [B] [P]
Yp_out = zeros(3, N_samples);
Yp_x1 = zeros(3, N_samples);
Yp_x2 = zeros(3, N_samples);
Yp_x3 = zeros(3, N_samples);
filename = sprintf('MonteCarloData%0d.mat',N_samples);
if( N_samples ~= 100 ...
&& N_samples ~= 1000 ...
&& N_samples ~= 10000 ...
&& N_samples ~= 100000 )
error('Use a number of samples that is a power of 10: (100, 1000, 10000, 100000). You used %0d. \n',N_samples);
end
if( exist(filename) && keep_existing_data )
fprintf('Found the data set for N_samples = %0d \n',N_samples);
fprintf('Loading... ');
load(filename);
fprintf('Done. \n\n');
else
% No MonteCarloData of same sample size N_samples was found.
if(keep_existing_data)
fprintf('No data set found for N_samples = %0d \n',N_samples);
else
fprintf('Ignoring data sets for N_samples = %0d \n',N_samples);
end
fprintf('Creating random samples... ');
count_m1 = 1;
count_m2 = 1;
for i=1:N_samples
% normalized values are between 0 and 1
r = rand(N_variables,1);
SamplesNormalized(:,i) = r;
% non-normalized values must be transformed
% from [0,1] to [alpha,beta]
% (and scale needs to be defined as linear or log)
mdot = 0; % dummy var
s = zeros(N_variables,1); % this particular sample
for j=1:N_variables
x_hat = r(j,1);
if( strcmp(names{j},'Mass flowrate') )
% mass flowrate is a multimodal variable
% and is therefore treated slightly differently
% (for more info see http://wiki.charlesmartinreid.com:8888/wiki/Example_Problem_for_Experimental_Design)
if( rand > 0.5 )
alpha(j) = 0.98;
beta(j) = 1.02;
mdot = 1;
else
alpha(j) = 1.98;
beta(j) = 2.02;
mdot = 2;
end
end
if( logscale(j) == true )
x = exp( ( log(beta(j)) - log(alpha(j)) )*x_hat + log(alpha(j)) );
else
x = (beta(j) - alpha(j))*x_hat + alpha(j);
end
s(j) = x;
end
% store the random sample in Samples array
Samples(:,i) = s;
% store the random sample in the appropriate mdot Samples array
if mdot == 1
Samples_mdot1(:,count_m1) = s;
count_m1 = count_m1 + 1;
else
Samples_mdot2(:,count_m2) = s;
count_m2 = count_m2 + 1;
end
end
fprintf('Done.\n\n');
fprintf('Evaluating function at random samples... \n');
loopTic = tic;
time_run_avg = 0;
times=zeros(N_samples,1);
% Don't use Matlab Parallel Computing Toolbox
%for i=1:N_samples
% Do use Matlab Parallel Computing Toolbox
% (http://www.mathworks.com/products/parallel-computing/demos.html)
count_m1 = 1;
count_m2 = 1;
parfor i=1:N_samples
ToyProblemTic = tic;
[Yp_x1(:,i) Yp_x2(:,i) Yp_x3(:,i) Yp_out(:,i)] = ...
ToyProblem_cmr( Samples(1,i), Samples(2,i), Samples(3,i), ...
Samples(4,i), Samples(5,i), Samples(6,i) );
fprintf('Evaluating sample %0d of %0d total... \n',i,N_samples);
ToyProblemToc = toc(ToyProblemTic);
times(i) = (time_run_avg)*((i-1)/i)+(ToyProblemToc/i);
if(Samples(1,i) <= 1.02 && Samples(1,i) >= 0.98)
Yp_x1_m1(:,count_m1) = Yp_x1(:,i);
count_m1 = count_m1 + 1;
end
if(Samples(1,i) <= 2.02 && Samples(1,i) >= 1.98)
Yp_x1_m2(:,count_m2) = Yp_x1(:,i);
count_m2 = count_m2 + 1;
end
end
loopTime = toc(loopTic);
fprintf('Done.\n\n');
fprintf('Total loop time: \t\t%0.4f s\n', loopTime);
fprintf('Function evaluation time average: \t%0.4f s\n', mean(times));
filename = sprintf('MonteCarloData%0d.mat',N_samples);
save(filename);
end
% ------------------------------------------------
% Postprocessing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make response surfaces
% (plots are created from MakeResponseSurface function)
% Perform variable transforms before calling MakeResponseSurface
% e.g.: Samples_mdot2(plotX,:) = Samples_mdot2(plotX,:).^(-1);
% Make quadratic response surface with all 6 input parameters
plotX = 2; plotY = 3;
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x1_mdot2(3,:), logscale,names,2, 'MCResponseSurface_Yp_x1_6dim_2deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x2_mdot2(3,:), logscale,names,2, 'MCResponseSurface_Yp_x2_6dim_2deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x3_mdot2(3,:), logscale,names,2, 'MCResponseSurface_Yp_x3_6dim_2deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_out_mdot2(3,:),logscale,names,2,'MCResponseSurface_Yp_out_6dim_2deg',plotX,plotY);
close all;
% Make custom cubic response surface with all 6 input parameters
model = [allVL1(6,2,'<='); eye(6)*3; 1 1 1 0 0 0; 0 1 2 0 0 0; 0 2 1 0 0 0 ];
plotX = 2; plotY = 3;
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x1_mdot2(3,:), logscale,names,model, 'MCResponseSurface_Yp_x1_6dim_3degmod',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x2_mdot2(3,:), logscale,names,model, 'MCResponseSurface_Yp_x2_6dim_3degmod',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x3_mdot2(3,:), logscale,names,model, 'MCResponseSurface_Yp_x3_6dim_3degmod',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_out_mdot2(3,:),logscale,names,model,'MCResponseSurface_Yp_out_6dim_3degmod',plotX,plotY);
close all;
% Make cubic resopnse surface
plotX = 2; plotY = 3;
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x1_mdot2(3,:), logscale,names,3, 'MCResponseSurface_Yp_x1_6dim_3deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x2_mdot2(3,:), logscale,names,3, 'MCResponseSurface_Yp_x2_6dim_3deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x3_mdot2(3,:), logscale,names,3, 'MCResponseSurface_Yp_x3_6dim_3deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_out_mdot2(3,:),logscale,names,3,'MCResponseSurface_Yp_out_6dim_3deg',plotX,plotY);
close all;
% Make quartic response surface
plotX = 2; plotY = 3;
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x1_mdot2(3,:), logscale,names,4,'MCResponseSurface_Yp_x1_6dim_4deg', plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x2_mdot2(3,:), logscale,names,4,'MCResponseSurface_Yp_x2_6dim_4deg', plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_x3_mdot2(3,:), logscale,names,4,'MCResponseSurface_Yp_x3_6dim_4deg', plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples_mdot2,Yp_out_mdot2(3,:),logscale,names,4,'MCResponseSurface_Yp_out_6dim_4deg',plotX,plotY);
close all;
% Make response surface with 2 input parameters:
% -> Lmix
% -> k(T)
% (this is much easier to plot...)
yy=2;
zz=3;
X = [Samples_mdot2(yy,:); Samples_mdot2(zz,:)];
logscale2(1) = logscale(yy);
logscale2(2) = logscale(zz);
names2 = cell(2,1);
names2{1} = names{yy};
names2{2} = names{zz};
response_surface_2 = MakeResponseSurface(X,Yp_x1_mdot2(3,:), logscale2,names2,2, 'MCResponseSurface_Yp_x1_2dim_2deg');
response_surface_2 = MakeResponseSurface(X,Yp_x2_mdot2(3,:), logscale2,names2,2, 'MCResponseSurface_Yp_x2_2dim_2deg');
response_surface_2 = MakeResponseSurface(X,Yp_x3_mdot2(3,:), logscale2,names2,2, 'MCResponseSurface_Yp_x3_2dim_2deg');
response_surface_2 = MakeResponseSurface(X,Yp_out_mdot2(3,:),logscale2,names2,2,'MCResponseSurface_Yp_out_2dim_2deg');
close all;
See also
| |||||||||||||||||||||||||||||||||||||||||
