% Class Example - car engine RSM calculations
% Barukyah Shaparenko and John M. Cimbala

clear all
close all

% Given:
% a = spark plug gap, in inches
% b = spark plug timing, before TDC, in degrees
% c = fuel:air mixture set screw, in number of turns of the screw
% P = power in horsepower from the engine

% To do: Calculate the direction of steepest ascent

% Solution: 

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

% Initial operating point: a=a_0, b=b_0, c=c_0
a_0 = 0.030;
b_0 = 10;
c_0 = 0;

del = [0.001; 1; 0.25]; % Initial deltas for "tweaking"

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

% 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));

P = [read_data{1,end}]; % From text file
% --OR--
%P = xlsread('sample_car_engine_RSM.xls',1,'F16:F21'); % From Excel file
% --OR--
%P = [156;162;159;154;158;160]; % From direct input

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('Physical variables\t\t\tResult\t\tCoded variables\n');
fprintf('a (in)\tb (deg)\tc (turns)\tP (hp)\t\tx_a\t\tx_b\t\tx_c\n');
for k = 1:length(a)
    fprintf('%-7g\t%-7g\t%-8g\t%-8g\t%-3g\t\t%-3g\t\t%-3g\n',...
        a(k),b(k),c(k),P(k),a_coded(k),b_coded(k),c_coded(k));
end

% Now perform a regression analysis on these normalized coded variables:
N = length(a_coded);
X = [ones(N,1) a_coded b_coded c_coded];
slopes = X\P; % Regression coefficients - vector of steepest ascent in coded variables

incr = zeros(length(slopes)-1,1);
incr_coded = zeros(length(slopes)-1,1);

[max_slope I] = max(abs(slopes(2:end))); % Find maximum magnitude coded variable

fprintf('\nSince coded variable x_%s has the maximum absolute value of slope,\n',parameters(I));
fprintf('we choose parameter x_%s as the one for which to choose the increment.\n',parameters(I));

incr_coded(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(I+1)
    incr_coded(I) = -incr_coded(I); % Sets correct direction of march: uphill
end

fprintf('Increment in x_%s = %g (unitless)\n',parameters(I),incr_coded(I));

fprintf('\nWe calculate the increments of the other two coded variables,\n');
fprintf('based on the direction of steepest ascent.\n');

for i=1:length(incr_coded)
    if i~=I
        incr_coded(i) = (incr_coded(I)/slopes(I+1))*slopes(i+1);
        fprintf('Increment in x_%s = %g\n',parameters(i),incr_coded(i));
    end
end

fprintf('\nNow convert back to the physical variables:\n');

incr(I) = incr_coded(I)*range(I)/2;

for i=1:length(incr)
    if i~=I
        incr(i) = incr_coded(i)*range(i)/2;
    end
end

fprintf('Increment in %s = %-10.3g (inches of spark plug gap)\n',parameters(1),incr(1));
fprintf('Increment in %s = %-10.3g (degrees before TDC)\n',parameters(2),incr(2));
fprintf('Increment in %s = %-10.3g (turns of the screw)\n',parameters(3),incr(3));

fprintf('\nFinal result:\n');
fprintf('March in direction of steepest ascent, i.e.\n');
fprintf('delta %s = %-12.5g (inches of spark plug gap)\n',parameters(1),incr(1));
fprintf('delta %s = %-12.5g (degrees before TDC)\n',parameters(2),incr(2));
fprintf('delta %s = %-12.5g (turns of the screw)\n',parameters(3),incr(3));
