Code

The following is the Matlab code for my composite design driver.
M file is available here: http://files.charlesmartinreid.com/ExperimentalDesign/CompositeDesign.m
% This is the Composite Design driver
% for response surface creation

keep_existing_data=true;

% Active parameters:
% (see http://wiki.charlesmartinreid.com:8888/wiki/Monte_Carlo_Experimental_Design)
% I     U        
% x_i   0.5 +/- 0.02
%       1.5 +/- 0.02
%       2.5 +/- 0.02
% m_dot 1.5 +/- 0.5 
% k(T)  1.5 * 10^{ 0 +/- 2}     (Log scale)
% L_mix 0.3 - 3.0               (Log scale)

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
alph = zeros(N_variables,1);
alph(1) = 0.98;   % mdot
alph(2) = 1.0e-2; % k(T)
alph(3) = 0.3;    % Lmix
alph(4) = 0.48;   % z1 = 0.5 m
alph(5) = 1.48;   % z2 = 1.5 m
alph(6) = 2.48;   % z3 = 2.5 m

% Set maximum for variable ranges
beta = zeros(N_variables,1);
beta(1) = 1.02;  % mdot = 1
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

% Create composite design matrix
SamplesNormalized = ccdesign(N_variables,'center',1,'type','inscribed');
SamplesNormalized = SamplesNormalized';
Samples = zeros(size(SamplesNormalized));

% Create arrays for storing function outputs
% 3 = number of species concentrations, [A] [B] [P]
Yp_out = zeros(3, size(SamplesNormalized,2));
Yp_x1  = zeros(3, size(SamplesNormalized,2));
Yp_x2  = zeros(3, size(SamplesNormalized,2));
Yp_x3  = zeros(3, size(SamplesNormalized,2));


filename = sprintf('CompositeDesignData.mat');



% If a data set already exists, load it
if( exist(filename) && keep_existing_data )

    fprintf('Found an existing data set. \n');
    fprintf('Loading... ');
    load(filename);
    fprintf('Done. \n\n');



% Otherwise, make a new data set
else

    fprintf('No composite design data set found. \n');
    fprintf('Creating samples... ');

    for jj=1:size(SamplesNormalized,1) %which variable
        for kk=1:size(SamplesNormalized,2) %which variable value
            if( logscale(jj) == true )
                x_hat = SamplesNormalized(jj,kk);
                m2 = (log(beta(jj)) - log(alph(jj)))/2;
                x = exp( m2*x_hat + m2 + log(alph(jj)) );
            else
                x_hat = SamplesNormalized(jj,kk);
                m2 = (beta(jj) - alph(jj))/2;
                x = m2*x_hat + m2 + alph(jj);
            end
    
            Samples(jj,kk) = x;
        end
    end

    fprintf('Done. \n');

    fprintf('Evaluating function at random samples... \n');
    loopTic = tic;
    time_run_avg = 0;

    % Loop over all samples
    parfor i=1:size(Samples,2)
        
        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,size(Samples,2));

        ToyProblemToc = toc(ToyProblemTic);
        times(i) = (time_run_avg)*((i-1)/i)+(ToyProblemToc/i);

    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('CompositeDesignData.mat');

    save(filename);

end



% ------------------------------------------------



% Postprocessing



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make response surfaces
% (plots are created from MakeResponseSurface function)

% Make quadratic response surface with all 6 input parameters
plotX = 2; plotY = 3;
response_surface_6 = MakeResponseSurface(Samples,Yp_x1(3,:), logscale,names,2, 'CDResponseSurface_Yp_x1_6dim_2deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples,Yp_x2(3,:), logscale,names,2, 'CDResponseSurface_Yp_x2_6dim_2deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples,Yp_x3(3,:), logscale,names,2, 'CDResponseSurface_Yp_x3_6dim_2deg',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples,Yp_out(3,:),logscale,names,2,'CDResponseSurface_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,Yp_x1(3,:), logscale,names,model, 'CDResponseSurface_Yp_x1_6dim_3degmod',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples,Yp_x2(3,:), logscale,names,model, 'CDResponseSurface_Yp_x2_6dim_3degmod',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples,Yp_x3(3,:), logscale,names,model, 'CDResponseSurface_Yp_x3_6dim_3degmod',plotX,plotY);
response_surface_6 = MakeResponseSurface(Samples,Yp_out(3,:),logscale,names,model,'CDResponseSurface_Yp_out_6dim_3degmod',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(yy,:); Samples(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(3,:), logscale2,names2,2, 'CDResponseSurface_Yp_x1_2dim_2deg');
response_surface_2 = MakeResponseSurface(X,Yp_x2(3,:), logscale2,names2,2, 'CDResponseSurface_Yp_x2_2dim_2deg');
response_surface_2 = MakeResponseSurface(X,Yp_x3(3,:), logscale2,names2,2, 'CDResponseSurface_Yp_x3_2dim_2deg');
response_surface_2 = MakeResponseSurface(X,Yp_out(3,:),logscale2,names2,2,'CDResponseSurface_Yp_out_2dim_2deg');
close all;
CompositeDesign.m

From charlesreid1

Code

See also
