% Response Surface Methodology Example, n = 3 parameters
% Barukyah Shaparenko and John M. Cimbala

clear all
close all

% COMMENT: A Taguchi design array (3 parameters and 3 levels) is used below.
% The original operating condition, although not part of the Taguchi array,
% is also included in the regression analysis, since that data point is available.

% NOTE: Coded (i.e. normalized) variables x_1, x_2, and x_3 are used for
% the regression analysis instead of a, b, and c.

parameters = ['a';'b';'c'];

% Initial operating point: a=a_0, b=b_0, c=c_0
a_0 = 20;
b_0 = 6;
c_0 = -2;

del = [2; 0.5; 1]; % Initial deltas for "tweaking"

% Read in test array of level numbers from text file
fid1 = fopen('Example_RSM.txt');
read_data = textscan(fid1,'%d %d %d',7);    
Taguchi = [read_data{1,:}];

fprintf('\n----------------------\n');

for j = 1:10 % Specify number of iterations
    % Low, medium, high values used in Taguchi array
    a_values = [a_0-del(1); a_0; a_0+del(1)];
    b_values = [b_0-del(2); b_0; b_0+del(2)];
    c_values = [c_0-del(3); c_0; c_0+del(3)];

    % Create Taguchi array with actual values
    a = a_values(Taguchi(:,1));
    b = b_values(Taguchi(:,2));
    c = c_values(Taguchi(:,3));

    % Calculate the response y(a,b,c)
    y = 300-0.8*(a-15).^2-0.7*(b-10).^2-0.6*(c+5).^2 + 0.05*(a-15).*(b-10).*(c+5)+rand(1)*0.05-0.025;

    range = [max(a)-min(a); max(b)-min(b); max(c)-min(c)];
    mid_value = [(max(a)+min(a))/2; (max(b)+min(b))/2; (max(c)+min(c))/2];

    % Transform to coded variables
    a_coded = 2*(a-mid_value(1))/range(1);
    b_coded = 2*(b-mid_value(2))/range(2);
    c_coded = 2*(c-mid_value(3))/range(3);
    
    fprintf('\n----------------------\n');
    fprintf('Iteration #%g\n',j);
    fprintf('original parameters\t\t\tresponse\t\tcoded variables\n');
    fprintf('a\t\tb\t\tc\t\t\ty\t\t\t\tx_1\t\tx_2\t\tx_3\n');
    for k = 1:length(a)
        fprintf('%-7g\t%-7g\t%-8g\t%-8g\t\t%-3g\t\t%-3g\t\t%-3g\n',...
            a(k),b(k),c(k),y(k),a_coded(k),b_coded(k),c_coded(k));
    end
 
    % Perform a regression analysis with a_coded, b_coded, and c_coded as
    % the independent variables, and y as the dependent variable:
    N = length(a_coded);
    X = [ones(N,1) a_coded b_coded c_coded];
    slopes = X\y; % Regression coefficients - vector of steepest ascent in coded variables
    slopes_unit = slopes(2:end)/sqrt(sum(slopes(2:end).^2)); % Unit vector of steepest ascent

    fprintf('\nVector of steepest ascent in coded variables = (%g, %g, %g)\n',slopes(2:end));
    fprintf('Magnitude of this vector = %g\n',sqrt(sum(slopes(2:end).^2)));
    
    fprintf('\nUnit vector of steepest ascent in coded variables = (%g, %g, %g)\n',slopes_unit);
    fprintf('Verify: Magnitude of this vector = %g\n',sqrt(sum(slopes_unit.^2)));

    incr = zeros(length(slopes_unit),1);
    incr_coded = zeros(length(slopes_unit),1);
    
    % COMMENT: We find the coded variable with the largest magnitude, and 
    % choose the increment of its dimensional (uncoded) value.
    % The increments of the other two parameters will be calculated, based
    % on the direction of steepest ascent.
    
    [max_slope I] = max(abs(slopes_unit)); % Find maximum magnitude coded variable
    incr(I) = del(I); % Choose an increment in its uncoded parameter
    % NOTE: This increment can be chosen at a different value, if desired.
    if max_slope == -slopes_unit(I)
        incr(I) = -incr(I); % Sets correct direction of march: uphill
    end

    incr_coded(I) = 2*incr(I)/range(I); % Increment in coded variable

    fprintf('\nChoose increment in parameter %s = %g\t',parameters(I),incr(I));
    fprintf('This corresponds to increment in %s_coded = %g\n\n',parameters(I),incr_coded(I));

    % Now, use ratios, based on the direction of steepest ascent to
    % calculate (coded and uncoded) increments in the remaining variables:
    for i=1:length(incr)
        if i~=I
            incr_coded(i) = (incr_coded(I)/slopes_unit(I))*slopes_unit(i);
            incr(i) = incr_coded(i)*range(i)/2;
            fprintf('Increment in %s_coded = %-10g\t\t',parameters(i),incr_coded(i));
            fprintf('This corresponds to increment in %s = %g\n',parameters(i),incr(i));
        end
    end
    
    % If desired, choose the final increments by rounding off, 
    % for simplicity in conducting the experiment.
    
    % Now march "uphill" from the previous middle point until y starts to decrease.
    march = [1:100]'; % Creates 100 steps in the direction of marching.
    a_march = a_0+incr(1)*march;
    b_march = b_0+incr(2)*march;
    c_march = c_0+incr(3)*march;

    y_march = 300-0.8*(a_march-15).^2-0.7*(b_march-10).^2-0.6*(c_march+5).^2 ...
        + 0.05*(a_march-15).*(b_march-10).*(c_march+5)+rand(1)*0.05-0.025;

    [max_y ind] = max(y_march); % Find maximum y value along path
    
    fprintf('\na\t\t\tb\t\t\tc\t\t\ty\n');
    fprintf('%-10g\t%-10g\t%-10g\t%-10g\tStarting point\n',a_0,b_0,c_0,y(1));
    for m = 1:ind+4
        if m~=ind
            fprintf('%-10g\t%-10g\t%-10g\t%-10g\n',a_march(m),b_march(m),c_march(m),y_march(m));
        else
            fprintf('%-10g\t%-10g\t%-10g\t%-10g\tMax value (new starting location)\n',...
                a_march(m),b_march(m),c_march(m),y_march(m));
        end
    end
    
    % Now start over, centering around the parameters that yieled the maximum value:
    % NOTE: This time use smaller increments around the operating point,
    % since we are closer to the optimum.  Also, it is not necessary to
    % exactly center around the operating point, for convenience.
    
    % Decrease deltas for next iteration
    del = del*0.75;

    % Specify new operating condition
    a_0 = a_march(ind);
    b_0 = b_march(ind);
    c_0 = c_march(ind);
end

% FINAL COMMENTS:
% We stop here since the maximum y is now very close to the known optimum of 300.
fprintf('\nThe optimum operating condition calculated here is\n');
fprintf('y = %g at (a,b,c) = (%g, %g, %g).\n',max_y,a_march(ind),b_march(ind),c_march(ind));
fprintf('This is very close to the real optimum condition of y = 300 at (a,b,c) = (15,10,-5).\n');
% The optimization analysis can continue if desired, but we may have
% reached a point of diminishing returns.