%% ME 2016
%% Fall 2007
%%
%% AUTHOR:  Chris Paredis
%%
%% PURPOSE: systematically explore and solve an optimization problem for
%% finding the differential gear ratio at which a car drives as fast as
%% possible over a quarter-mile drag race
%%
%% ASSUMPTIONS: 
%%         - the same assumptions as are being made in the car simulation
%%         - the gear ratio can take on values between 1 and 20

clear all;
warning off;

% these two lines set up a global variable that can be used to determine
% the total number of evaluations of the function car_model_w_parameters,
% as asked for in question 4
global num_model_evals;
num_model_evals = 0;

%% Step 1: create a plot to get an approximate view of what the function
%% looks like for a slope of zero
% pick as an initial range the values from 1 to 20 in steps of 1
differential_ratios = 1:20;  
finish_times = zeros(length(differential_ratios),1);
for i=1:length(differential_ratios)
    finish_times(i)=evaluate_car(differential_ratios(i));
end

% plot the results
figure(1)
plot(differential_ratios,finish_times,'b');
xlabel('differential ratio []')
ylabel('finish time [s]')
title('finish time for quarter-mile drag race');
grid on

%% Step 2: perform the optimization using fminbnd
% reset the number of model evaluations for question 4
num_model_evals = 0;

% perform the actual optimization -- minimization of the finish time
[opt_ratio,opt_time, flags, output] = fminbnd(@evaluate_car, ...
                                      1, 20, optimset('Display','Iter'));
      
% print the results of the optimization
fprintf('the optimum occurs at a differential ratio of %g\n',opt_ratio);
fprintf('the optimum finish time is %g\n',opt_time);
fprintf('the solution required %d iterations and %d function evaluations\n',...
        output.iterations, output.funcCount);
fprintf('overall, the function car_model_w_parameters was called %d times\n',...
         num_model_evals);

%% Step 3: investigate the vicinity of the optimum
% use a range of +/- 5% with 100 equally spaced points
range = 0.05;
differential_ratios = linspace((1-range)*opt_ratio,(1+range)*opt_ratio, 100);
finish_times = zeros(length(differential_ratios),1);
for i=1:length(differential_ratios)
    finish_times(i)=evaluate_car(differential_ratios(i));
end

%% plot the results
figure(2)
plot(differential_ratios,finish_times,'b');
xlabel('differential ratio []')
ylabel('finish time [s]')
title('finish time for quarter-mile drag race');
grid on

%% Answers to questions:

%% Question 1: At which differential ratio does the minimum appear to
%% occur?
% the minimum appears to be somewhere between 7 and 8

%% Question 2: Does the problem satisfy all the conditions that are
%% necessary for golden section search (which is used by Matlab) to be
%% applicable? List at least three characteristics. 
% yes, at least appears so from this graph.  The characteristics are:
% 1) the function is uni-modal -- only has one minimum
% 2) the initial interval is such that a minimum occurs in it
% 3) the problem is one-dimensional -- only one independent variable
%
% The derivative of the function is not provided, but the derivative is not
% needed for golden section search.

%% Question 3: How many function evaluations of evaluate_car are necessary
%% to reach the optimum?
% Here is the output:
% Func-count     x          f(x)         Procedure
%    3        8.25735      14.8885        initial
%    4        12.7426      18.5824        golden
%    5        5.48529      14.9746        golden
%    6        7.00319      14.7723        parabolic
%    7        7.06186      14.7619        parabolic
%    8        7.42486      14.7291        parabolic
%    9        7.74284      14.7583        golden
%   10        7.41218      14.7281        parabolic
%   11        7.27837      14.7384        golden
%   12        7.38144      14.7297        parabolic
%   13        7.40537      14.7281        parabolic
%   14        7.40814       14.728        parabolic
%   15         7.4084       14.728        parabolic
%   16        7.40984       14.728        golden
%   17         7.4104       14.728        parabolic
%   18        7.41035       14.728        parabolic
%   19        7.41032       14.728        parabolic
% 
% Optimization terminated:
% the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 
%
% the optimum occurs at a differential ratio of 7.41035
% the optimum finish time is 14.728
% the solution required 16 iterations and 19 function evaluations
% overall, the function car_model_w_parameters was called 2929 times

% so the answer to the question is 19 -- NOTE if you started from a
% different interval, your numbers may be different

%% Question 5: At which differential ratio does the minimum appear to
%% occur? Include your graph in your report. Explain the shape of your
%% graph.
% the minimum occurs at 7.4066 which is different from the solution found
% previously!  The shape of the graph is very noisy because of the
% round-off errors that occur in the simulation.  Matlab uses an adaptive
% step-size in ode45 with a default desired accuracy of 0.1%.  In this case
% that corresponds to a a potential error in the result of approximately
% 14.7*0.001 = 0.015 seconds -- notice that this corresponds more or less
% to the size of the noise peaks that seem to be occuring.

%% Question 6: Does the problem satisfy all the conditions that are
%% necessary for golden section search (which is used by Matlab) to be
%% applicable? 
% oops... we were wrong in our initial assessment.  The function is not
% uni-modal at least not when investigated at a level where the noise is
% visible.

%% Question 7: Discuss the accuracy with which the optimum can be
%% determined. What are some of the factors that play a role in determining
%% the accuracy of the optimum? 
% given that the noise is approximately 0.015 seconds, optimization
% routines cannot distinguish between function values that differ by less
% than 0.015.  That implies that the true minimum (according to the model)
% could lie almost anywhere in the range between 7.1 and 7.8.  When
% considering that our model is only an approximation (All models are
% wrong, some are useful), our ability to find the "true" optimum is
% limited even further.  If we want to know the optimum differential ratio
% more accurately, we definitely will need to increase the accuracy of our
% model and take such phenomena as tire slip into account.
% conclusion:  be very careful when interpreting the results of an
% optimization problem:  given the inaccuracy in the model and in the
% numerical algorithms, you are virtually guaranteed that the "optimum"
% found by your algorithm is not the optimum.  Even worse -- given that all
% models are wrong, the optimum can never be found.