% Class 19 Design of Experiments
% Coded by NFR on 10.29.2019
%
% Characterization Designs (after you have narrowed
% down to significant variables,
% use this type of DoE to determine which
% variables to focus on for response surface
%
% Full Factorial
% 2 level (each parameter has a high and low)
Runs = ff2n(5);
% Also specify the levels of each (not 
% constrained to a two level design)
Runs2 = fullfact([2 3 4 2 3])
Nruns = length(Runs)
% This results in too many experiments
% so statisticians determined the minumum
% number of experiments needed to still
% have statistical confidence
% Fractional factorial designs
% We have 32 runs for full factorial design
% could we reduce this number?
% two level design that only looks at linear terms
GenD = fracfactgen('a b c d e');
Run3 = fracfact(GenD)
NumRuns = length(Run3)
% two level fractional design that has
% linear terms AND some interaction terms
% (you can decide which to include based
% on physcial understanding of system)
GenD2 = fracfactgen('a b c d e ab ac ae bc')
Run4 = fracfact(GenD2)
NumRuns = length(Run4)
% At this phase, this is the list of 
% experiments we need to run, and you would 
% fit to determine which variables to 
% focus on in response surface
% for now, we will jump to a response surface
% (main focus of our two days on DoE)
%
% Response Surface DoE problems (optimization)
% Step 1) Generate design (list of experiments)
%       a) use matlab to make a coded list
%       b) randomize order and convert to real levels
%       c) place centerpoints throughout the exp
%          (at least one at start, middle and end)
% Step 2) Run these experiments
% Step 3) Fit the results to a model (quadratic) 
%      a) which coeffcients are important
%      b) Use reduced model to make predictions
% Step 4) Test that optimal point -- if confirmed, celebrate
%
% Generate a design, we can use Box-Behnken or central composite
Runs5 = bbdesign(3)
NumRuns = length(Runs5)
% Note: multiple centerpoints allow us to determine
% if there is experimental drift
Runs6 = ccdesign(3,'type','faced')
NumRuns = length(Runs6)
% Example with car radiator fan
% P1 = distance to radiator
% P2 = pith of the blade
% P3 = clearance of the blade
bounds = [1 1.5; 15 35; 1 2];
% Create experimental list
CodedVal = bbdesign(3);
% Randomize the order of experiments
ReorderVec = randperm(length(CodedVal));
% We also want to convert coded units to real units
RealValues = zeros(size(CodedVal));
for i = 1:3
    zmin = min(CodedVal(:,i));
    zmax = max(CodedVal(:,i));
    RealValues(:,i) = interp1([zmin zmax],...
        bounds(i,:),CodedVal(:,i));
end
% Next step is to randomize the list to hand 
% off for experiments
disp({'Run Number','Distance','Pitch','Clearance'})
disp(sortrows([ReorderVec' RealValues]))
% Manually check that your center points
% are distributed throuought
%
% Now you run the tests, get data and optimize!
% Results in non-randomized order
Outcomes = [837 864 829 856 880 879 872 874 834 833 860 859 874 876 875];
% Typically fit these response surfaces
% initially with a full quadratic model (except 
% for the three parameter term)
% (see paper for equation in yellow)
% Before we used fittype to specify such models
% 1) this is big equation and 2) fittype only likes up to 2 independent
% parameters
%
% We will use fitlm (fit linear model)
mdl = fitlm(RealValues,Outcomes,'quadratic')
% Upon inspection we see some coefficients
% have a very high p-value, we should remove these
% NOTE: you want to remove them one by one, or in small numbers
mdl2 = removeTerms(mdl,'x1:x2+x2:x3')
% NOTE: you can keep removing terms until
% you are left with just signficant ones
%
% Let's show how you would optimize
% (full model)
P = mdl.Coefficients.Estimate;
f = @(x) -1*(P(1) + P(2)*x(1) + P(3)*x(2) + P(4)*x(3)...
    + P(5)*x(1)*x(2) + P(6)*x(1)*x(3) + P(7)*x(2)*x(3)...
    + P(8)*x(1)^2 + P(9)*x(2)^2 + P(10)*x(3)^2);
% Find optimum
xo = [1 25 1];
fminsearch(f,xo)
% We need bound the search within the model space
lb = bounds(:,1);
ub = bounds(:,2);
[x,feval] = fmincon(f,xo,[],[],[],[],lb,ub)
% What happens when the model has poor coeffcient p-values?
mdl2 = removeTerms(mdl,'x1:x2+x2:x3');
mdl3 = removeTerms(mdl2,'x1:x1+x3:x3');
mdl4 = removeTerms(mdl3,'x1:x3+x3');
% Plot the reduced model and compare to the full model!
P = mdl4.Coefficients.Estimate;
f2 = @(x) -1*(P(1) + P(2)*x(1) + P(3)*x(2)+ P(4)*x(2)^2);
[x2,feval2] = fmincon(f2,xo,[],[],[],[],lb,ub)