%% Finding the Optimal Value
% Have you ever needed to solve an optimization problem where there were local minima?  What strategy do you use to solve it, trying to find the "best" answer?  Today I'm going to talk about a simple strategy, readily available in the Global Optimization Toolbox.
%
%% Solve a Simple Problem
% Or at least let's try.  I have some data and I want to fit a particular form of a curve to it.   First let's look at the pharmacokinetic data. Here's the reference: Parameter estimation in nonlinear algebraic models via global optimization. Computers & Chemical Engineering, Volume 22, Supplement 1, 15 March 1998, Pages S213-S220 William R. Esposito, Christodoulos A. Floudas.
% The data are time vs. concentration

t = [ 3.92,  7.93, 11.89, 23.90, 47.87, 71.91, 93.85, 117.84 ]
c = [0.163, 0.679, 0.679, 0.388, 0.183, 0.125, 0.086, 0.0624 ]

% I like to see the data, in part to be sure I have no entry mistakes, and in part to get a feel for the overall system.  In fact, let's visualize the data.

plot(t,c,'o')
xlabel('Time')
ylabel('Concentration')

%
%% 3 Compartment Model
% As in the reference, we fit a 3 compartment model, sum of 3 decaying exponentials.
% c=b_1 \;e^{\left(-b_4 \;t\right)} +b_2 \;e^{\left({-b}_5 \;t\right)} +b_3 \;e^{\left({-b}_6 \;t\right)}
% and we can express that model as an anonymous function of t (time) and the model parameters [b(1) b(2) ... b(6)].

model = @(b,t) b(1)*exp(-b(4)*t) + b(2)*exp(-b(5)*t) + b(3)*exp(-b(6)*t)

%
%% Define Optimization Problem
% We next define the optimization problem to solve using the problem-based formulation.  This allows us to choose the solver we want, supply the data, and naturally express constraints and options.

problem = createOptimProblem('lsqcurvefit', ...
    'objective', model, ...
    'xdata', t, 'ydata', c, ...
    'x0',ones(1,6),...
    'lb', [-10 -10 -10  0   0   0 ],...
    'ub', [ 10  10  10 0.5 0.5 0.5], ...
    'options',optimoptions('lsqcurvefit',...
    'OutputFcn', @curvefittingPlotIterates,...
    'Display','none'))

%
%% Solve the Problem
% First solve the problem directly once.

b = lsqcurvefit(problem)

% You'll notice that the model does not do a stellar job fitting the data or even following the shape of the data.
%
%% MultiStart
% Let's see if we can do better by starting at a bunch of different points.

ms = MultiStart;
ms.Display = 'iter';
rng default
figure
tic
[~,fval,exitflag,output,solutions] = run(ms, problem, 50)
serialTime = toc;

%
%% Visualize the Best Solution
% The 50th solution, which is what is plotted above, is not necessarily the best one.  Luckily for us, MultiStart orders the solutions from best to worst. So we need only look at the first one.

curvefittingPlotIterates(solutions)

% You can see now that the 50th was not the best solution as the mean squared error on this final one displayed is over a factor of 10 better.
%
%% MultiStart with Parallel Computing
% I will now see if we can improve the performance using all 4 of my cores as parallel workers locally.

ms.UseParallel = true;
gcp;
tic;
rng default
run(ms, problem, 50);
parallelTime = toc;

%
%% Calculate Speedup
% Speed up may not be evident until second run due to pool start up time.  Since I started mine earlier, I get to see decent speed up.

speedup = serialTime/parallelTime

%
%% Do You Have Problems Where the Solution is Sensitive to the Starting Point?
% Tell us about your exploration of the solution space to your problem.  If the solution is sensitive to where you start, you might consider using MultiStart and other techniques from the Global Optimization Toolbox.
%
% Copyright 2021 The MathWorks, Inc.
%
%% Appendix
% Here's the code for plotting the iterates.

dbtype curvefittingPlotIterates