function PracDOE
% Coded by NFR on 10.23.17
% This code how MATLAB can be used for DOE problems
%
% Factors
% - Distance from radiator
% - Pitch angle
% - Blade tip clearance
%
CodedValue = bbdesign(3);
m = 10;
%Randomize run order and set up real world limits
runorder = randperm(15);
bounds = [1 1.5; 15 35; 1 2];
RealValue = zeros(size(CodedValue));
for i = 1:3
    zmax = max(CodedValue(:,i));
    zmin = min(CodedValue(:,i));
    RealValue(:,i) = interp1([zmin zmax],bounds(i,:),CodedValue(:,i));
end
TestResult = [837 864 829 856 880 879 872 874 834 833 860 859 874 876 875]';
disp({'Run Number','Distance','Pitch','Clearance','Airflow'})
disp(sortrows([runorder' RealValue TestResult]))
% SKIP the table saving bit
% How to optimize? Fit quadratic model
mdl = fitlm(RealValue,TestResult,'quadratic') % Don't suppress and you get model parameters
% Look at p values to determine dominant terms
% plotSlice(mdl) % use for interactive response surface plot...
% You could use this to try and find the max, but let's use our existing
% tools...we first need to define our fit model as an anonymous function.
% You could code out each of the model parameter fit estimates
%
P = mdl.Coefficients.Estimate;
%f = @(x) -x2fx(x,'quadratic')*P % This will be confusing to all!
f = @(x) -(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);
xo = [1.25 27 1.5];
lb = bounds(:,1);
ub = bounds(:,2);
[xout, feval] = fmincon(f,xo,[],[],[],[],lb,ub)
% Perfect!  However the 'optimal' point exists at the boundary...is this
% good? No...would need to redo experiments and recenter the bounds
%
% You can use monte carlo to look at product spec tolerance
%
% Use the rsmdemo for second example
% Hydrogen 100 to 470
% nPentane 80 to 300
% isoPentane 10 to 120
CodedValue = ccdesign(3,'type','inscribed');
runorder = randperm(length(CodedValue));
bounds = [100 470; 80 300; 10 120];
RealValue = zeros(size(CodedValue));
for i = 1:3
    zmax = max(CodedValue(:,i));
    zmin = min(CodedValue(:,i));
    RealValue(:,i) = interp1([zmin zmax],bounds(i,:),CodedValue(:,i));
end
% from the rsmdemo
TestResult = [6.0358 2.4873 12.0536 6.7753 3.4668 1.8512 7.7262 4.8287 7.3329 3.8373 1.4935 8.3146 8.1546 3.3491 3.3407 4.9223 5.0017 5.0560 5.1173 5.1664 5.2539 4.4553 5.1306 5.5768]';
disp({'Run Number','Hydrogen','n-Pentane','iso Pentane','Reaction Rate'})
disp(sortrows([runorder' RealValue TestResult]))
% SKIP the table saving bit
% How to optimize? Fit quadratic model
mdl = fitlm(RealValue,TestResult,'quadratic')
mdl.VariableNames
%plotEffects(mdl)
%plotInteraction(mdl,'x1','x2')
%plotSlice(mdl)
P = mdl.Coefficients.Estimate;
%f = @(x) -x2fx(x,'quadratic')*P % This will be confusing to all!
f = @(x) abs(15-(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));
xo = [1.25 27 1.5];
lb = bounds(:,1);
ub = bounds(:,2);
[xout, feval] = fmincon(f,xo,[],[],[],[],lb,ub)
m = 10;