%% Optimization Tips & Tricks
% This script contains the examples shown in the webinar titled
% _Optimization Tips and Tricks: Getting Started using Optimization with
% MATLAB_ presented live on 21 August 2008.  To view the webinar, please go
% <http://www.mathworks.com/webinars/ here> and click on recorded
% webinars.
%
% I've tried to make the demos self explanatory in this file.  But if
% something is broken or not clear, please let me know.
%
% To run the demos in this file, you need Optimization Toolbox and Genetic
% Algorithm and Direct Search Toolbox, which I will refer to as GADS.
% These files work with Release R2008a.  You will also need the Parallel
% Computing Toolbox to run hte parallel examples.  Most demos should work
% with older releases, but some will require R2008a products.  In
% particular, the parallel optimization examples require it.  You can
% obtain a demo copy from
% <http://www.mathworks.com/products/optimization/tryit.html> or 
% <http://www.mathworks.com/products/gads/tryit.html>
%% Introduction to Optimization - Mt. Washington Demo
% This first example shows you how optimization in general works.
% Optimization solvers are domain searching algorithms that differ in the
% types of problems (or domains) they can solve (search).  This example
% shows how the pattern search algorithm can be used to find the peak of
% the White Mountain Range.
%
% Note: this demo is available in the GADS Toolbox.

% Load data for white mountains
load mtWashington;   
% Load psproblem which have all required settings for pattern search
load psproblem;
% Run the solver and show displays
patternsearch(psproblem);
%%
% This command will load the example in optimization tool, a graphical user
% interface for setting up and running optimization problems.  You can
% access optimization tool from the *Start Menu --> Toolboxes -->
% Optimization --> Optimization Tool (optimtool)* or by typing |optimtool|
% at the command prompt.


% Remove the % symbol if you'd like to run this part of the code.
% mtwashingtondemo

%%
% You should see two plots when running the Mt. Washington demo.  One is a
% surface plot,  the other is a topology map.  On the topology map you can
% see that starting point and iterations (filled circles) and the tested
% points that were not selected (+ symbols).  There is also a slider bar on
% the right that you can use to speed up/slow down the process.
%
% Notice how the pattern search solver expands and contracts the search
% radius as it explores the domain for the maximum value.  This is an
% example of a patterned search, as the name implies, and is only one of
% many search patterns the pattern search solver supports.
%
% I showed this example to illustrate how effective the pattern search
% solver can be on a highly rough surface like this one, which gradient
% based solvers like the ones in Optimization Toolbox couldn't solve.
%
% This leads to the primary difference between the Optimization Toolbox
% solvers and the GADS solvers.  Optimization Toolbox solvers rely on
% gradients (well, most of them do anyway, |fminsearch| does not) to
% determine the search direction.  GADS solvers use gradient free methods
% to find search directions and optimal values.  This makes GADS solvers
% more useful on highly nonlinear, discontinuous, ill-defined, or
% stochastic problems. They are also better at finding global solutions,
% although they do not guarantee a global solution will be found.
% 
% The price you pay is that the GADS solvers are often slower and are not
% able to handle as large of problems as the Optimization Toolbox solvers.
%% Getting Started - MATLAB Backslash
% Let's start with a simple curve fitting example.  We have time (t) and a
% response (y).  Looking at the plot, you can see that the curve follows an
% exponential decay.
close all, clear all
t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
plot(t,y,'o')

%%
% Let's assume that an equation of this from can be used to describe the
% data
%
% $$yhat = c_1 + c_2e^{-t}$$
% 
% We can linearize this problem and solve it in MATLAB using the backslash
% operator.
E = [ones(size(t)) exp(-t)]
c = E\y

%% 
% Notice that we have an over determined system.  The neat thing about
% backslasch (\) is that when you try to solve an over determined system
% (more equations than unknowns), it tries to minimize the error when
% estimating the solution.  It's an optimization solver (can be used as a
% least squares minimizer).
% 
% Let's plot the results.  Not a bad fit.   
yhat = @(c,t) c(1) + c(2)*exp(-t); % I'll describe function handles in a bit
ezplot(@(t)yhat(c,t),[0,2.5])
hold on
plot(t,y,'o')
hold off
%% Now Solve using |lsqnonneg| optimization solver
% Let's try solving this with an optimization solver.  Say lsqnonneg.
c = lsqnonneg(E,y)
%%
% We get the same result. Note that lsqnonneg solves a linear system of
% equations using least squares minimization.  If you haven't already, take
% a look at the optimization documentation for lsqnonneg.  And if you have
% R2008a or later, take a look at the new documentation (Optimization
% Toolbox --> Optimization Overview --> Choosing a solver).  This page
% gives a nice summary of the optimziation solvers, their problem types,
% and a decision matrix/table at the bottom of the page.  Look at the
% least-squares solvers and note the different formulations.
%
% That wasn't to hard, but let's now add a constraint to our curve fitting
% problem and make it a bit more challenging (i.e. something that Curve
% Fitting Toolbox doesn't currently support.
%% Constrain the y-intercept
% Let's say that we know the data value at time 0 precisely.  That is, we
% can set the starting condition precisely for this exponentially decaying
% system.  This is a constraint, and we can add it to our least squares minimization
% problem that we've developed so far.
%
% This constraint is a linear constraint on c(1) and c(2).
% At t = 0, c(1) + c(2)*1 = 0.8 --> Aeq*x = beq --> [1 1][c1 c2]' = 0.8
% I've written the constraint in the format accepted by the Optimization
% Toolbox Solvers.  And if we look at the doc for the solver we need, we'd
% find that the solver that matches our problems is |lsqlin|.

Aeq = [1 1];
beq = 0.82;
cc = lsqlin(E,y,[],[],Aeq,beq)
%%
% Now let's plot it against our data and the unconstrained fit earlier.  We
% can see that it has a slightly different shape, and if you were to zoom
% into the y-intercept, you'd see that it passes through our data point
% where the unconstrained fit did not.
ezplot(@(t)yhat(c,t),[0,2.5])
hold on
tt = 0:0.1:2.5;
plot(tt,yhat(cc,tt),'g')
plot(t,y,'o')
hold off
legend('Unconstrained Fit', 'Constrained Fit', 'Data')
%% Working with function handles
% The solver I showed earlier accepted MATLAB data types, arrays, to define
% the problem and constraints.  However, most of the optimization solvers
% require that a function be passed in as the objective function(also true
% for nonlinear constraints).  The mechanism that allows you to pass a
% function into another function is known as a function handle.
%
% A function handle keeps reference to your function and allows you to pass
% a function into another function.
%
% Let's redefine our equation as an anonymous function (an anonymous
% function is simply a function that exists only in memory - not in an
% M-fle).
yhat = @(c,t) c(1) + c(2)*exp(-t)
yhat(cc,t) % pass in the previous solution and the time vector
%%
% Note the @ symbol.  This declares a function hanlde, and you can see that
% it behaves like a function call.  Note that I can even redefine my
% function to only take in one input field (this will be important later).
% Let's redefine yhat to only be a function of c. 
yhat2 = @(c) yhat(c,t)
yhat2(cc)
%%
% Notice that when I redefined yhat as yhat2, the time vector t is not
% passed.  How does MATLAB know how to use t in this case?  When I
% redefined the function to accept only c, MATLAB stored the time vector in
% the function handle object, so it can use it.
%
% To prove this point, notice what happens if I redefine t, yhat2 does
% not use this new value since it is still referencing the t I had defined
% at the time I declared yhat2.  So t became a constant vector, and was
% saved with the function handle definition of yhat2.
%
% This is often a point of confusion when first starting to use function
% handles and it's good to be aware of how parameters get saved when
% creating a function handle.
t = 0:0.1:2.5; % has more points than before.
y2 = yhat2(cc)
y1 = yhat(cc,t')
%% Solve using |lsqcurvefit|
% Ok, so let's get back to the reason I introduced function handles.  You
% need them for most of the nonlinear optimization solvers, as they expect
% you to pass in a function as the objective or nonlinear constraint. 
%
% |lsqcurvefit| expects a function to be passed in.  Here I pass it yhat
% (note that it is already a function handle in the workspace, if I wanted
% to pass in an M-file function called yhat, I need to pass it as @yhat).
% Also take notice that |lsqcurvefit| allows you to pass in the data to fit
% to, in this case t and y.  You can see we get the same result as in the
% unconstrained case (this solver does not allow linear constraints, that's
% why we're using the unconstrained case). 
t = [0 .3 .8 1.1 1.6 2.3]';     % define again since we changed it earlier
clear yhat2;                    % no longer needed
c = lsqcurvefit(yhat,[1 1],t,y)

%% Solve using |lsqnonlin|
% |lsqcurvefit| is one of the few solvers that let you pass in an objective
% function with more than one field.  But most of our solvers only allow
% you to pass in one input field - the decision vector.  |lsqnonlin| is an
% example.  Notice in this case that I redefine yhat to be only a function
% of the decision variables, c.  It is also different from yhat2 I defined
% earlier in that it is subtracting out the y data - i.e. returning only
% the residuals of the fit.  So y and t are passed in as constant
% parameters to this optimization problem.
%
% This is an example of how you can use function handles to pass in data
% that your objective function may need to perform a calculation but is
% considered constant for the optimization problem.
%
% These results agree with prior ones.
c = lsqnonlin(@(c)yhat(c,t)-y,[1 1])
%% Passing data using function handles with an M-file objective function
% This example shows how to pass in the constant data t and y when using a
% M-file function called myCurveFcn (it's the same as yhat, but in an
% M-file).
c = lsqnonlin( @(c)myCurveFcn(c,t,y), [1 1])
%% Nested functions
% While function handles are very useful for passing constant parameters,
% there is another method that works well.  This example shows how to use
% nested function to pass data around in functions.  Note that the
% objective function doesn't get passed t and y into it from the function
% calling syntax.  The nested function has access to it's parent's function
% workspace, this is where it gets values of t and y that it needs to use
% in the body of the objective function.
type nestedOptimExample
c = nestedOptimExample()
%% Nonlinear Optimization and Topology Considerations
% Let's switch gears now, and talk about nonlinear optimization.  PEAKS is
% a function that ships with MATLAB and is shown below.  It is a highly
% nonlinear function and note how close the peaks and valleys are together.
% This is a very challenging surface to minimize using a gradient based
% solver, as we'll soon see.
peaks
%%
% Let's define our optimization problem.  I want to find the minimum value
% of the PEAKS function.  My objective function is the PEAKS function,
% which I've redefined in a M-file function called peaksObj.  It's been
% redefined here so that it can accept data in the format that the
% optimization solvers prefer.  I'd also like to add a nonlinear constraint
% to this problem, so you can see how to add them for your problem.  I've
% defined my nonlinear constraint in a function called peaksCon.
clear all, close all
type peaksObj
type peaksCon
%%
% We'll us the constrained nonlinear minimization solver, |fmincon|, to
% solve this problem.  Let's define the starting point and solve the
% problem.  (Note that I've created an output plotting function so we can
% see the solver's iteration history and defined this using solver
% options).
x0 = [0 1.5];
options = optimset('Disp','iter','OutputFcn',@peaksOutputFcn);
x = fmincon(@peaksObj,x0,[],[],[],[],[],[],@peaksCon,options)
%% 
% The top plot is a surface plot of the peaks function.  The lower plot is
% a contour plot of the upper plot with the constraint boundary shown in
% the shaded blue region.  This is where we violate the constraint.
% 
% Notice how fast it converged to a solution.  Looks good right?  Let's
% see what happens when we change the staring point.
x0 = [-3 -3];
x = fmincon(@peaksObj,x0,[],[],[],[],[],[],@peaksCon,options)
%%
% Bummer!  We didn't get the right solution.  We did get a local minima,
% but not the best or global minima.  Is this a problem with |fmincon|?
% No, this solver was designed to be fast and efficient at finding local
% minima.  It wasn't designed to be a global solver.  Most optimization
% solvers are only local solvers.  This is because it is very difficult to
% guarantee a global solution for nonlinear optimization problems like this
% one.
%
% This problem also illustrates that the choice of starting point is often
% critical to finding the right solution.
%% Try a random starting point (uniform distribution grid of 4 points)
% For most practical nonlinear problems, you probably will not run into
% this issue.  But if you suspect that your problem may be this nonlinear,
% we can get around this starting point issue using multiple staring points
% (commonly referred to as multistart optimization).  Let's do a simple
% example.  I'll run the same optimization problem from four different
% staring points.  How should I choose them?  There are a variety of
% methods you can choose from (such as design of experiments, space filling
% designs, grid patterns, or any other you can think of).  I'm simply going
% to generate four random starting points ands see how good this works.
% I'll use a uniform random number generator (all points are equally likely
% to be selected) and bound the starting points in the range of +/- 3 on
% each axis.
a = -3;
b = 3;
x1 = a + (b-a) * rand(2);
x2 = a + (b-a) * rand(2);
x0r = [x1; x2];

for i = 1:length(x0r)
    xopt(i,:) = fmincon(@peaksObj,x0r(i,:),[],[],[],[],[],[],@peaksCon,options);
end
xopt
%%
% Now this result will be different each time you run it, but you should
% see that in this case we found the global minima at least once (possibly
% more).  So multistart can be an effective method for finding the global
% minima.  Also note that I get additional minima when using this method.
% This may also be useful, depending upon your problem.  For example, the
% global minima may not be a stable condition when you take uncertainty in
% X1/X2 into consideration and you migh prefer the other minima in this
% case.
%
% I used only 4 points, but depending upon your problem and the size of
% your domain you would like to search you may need many more points.  So
% choose your multistart strategy carefully and incorporate any knowledge
% you may know about your problem to guide you on deciding how many point
% you need.
%% Solve using genetic algorithms
% One of the benefits of the solvers in the Genetic Algorithm and Direct
% Search Toolbox is that they are useful for problems that are highly
% nonlinear, such as this one.  They also are useful for problems that are
% stochastic, contain custom data types such as integer sets, problems
% that have undefined derivates, holes in the search domain, or are
% ill-defined.  They also are better at finding global solutions than the
% gradient based solver in Optimization Toolbox.
%
% But they are more expensive and are limited to solving problems of
% smaller size.
%
% Let's solve this problem using genetic algorithms (GA).
options = gaoptimset('Disp','iter','OutputFcns',@peaksOutputFcn);
x = ga(@peaksObj,2,[],[],[],[],[],[],@peaksCon,options)
%%
% The plot shows how GA progresses towards a solution.  It uses a
% stochastic search, based upon an initial population of points (green),
% and the selects a subpopulation (black) for creating a new population
% to evaluate.  GA get's it's name from it's search method, which
% selectively generate new candidate points to evaluate based upon a method
% that is similar to breeding among two 'parents' to generate a 'child'.  I
% won't get into the details here, but if you want to know more about GA
% solvers, I suggest you take a look at our documentation, or watch the
% Genetic Algorithms for Financial Application Webinar (has a nice
% description of GA algorithms in it).
%
% You can see that the final population (red) contains the global minima.
%% Solve using simulated annealing
% Another popular solver is Simulated Annealing.  This solver is a little
% more cost effective then GA, but note that it doesn't allow you to solve a
% problem with nonlinear constraints.  But for our problem, we can
% approximate the nonlinear constraint as a rectangular region, vs. a
% circular one.  So we'll use the lower/upper bounds to add this
% constraint.  Notice how this solver works.  It starts out be randomly
% sampling the domain and gradually reduces the search radius around a
% mimima that it's found.  Periodically, the solver will reset and start
% searching a wider radius.  This is the mechanism that allows this solver
% to jump out of a local minima and find the global minima.
lb = [-3 -3];
ub = -lb;
options = saoptimset('Display','iter','OutputFcns',@peaksOutputFcn);
x = simulannealbnd(@peaksObj,x0,lb,ub,options)

%% Solve using pattern search
% Let's now try out pattern search and see how well it does.  I explained
% |patternsearch| earlier so I won't get into details here.  I will however
% note that this solver's default settings are set to find the local
% optimum.  We'll add a search option that will make it more global.
options = psoptimset('Display','iter','OutputFcns',@peaksOutputFcn,...
                    'CompleteSearch','on','SearchMethod',@searchlhs);
x = patternsearch(@peaksObj,x0,[],[],[],[],[],[],@peaksCon,options)
%% Which solver should I choose?
% Let's time these approaches and see how they fare.
compareApproaches
%%
% From this chart it’s clear that using a gradient based solver, such as
% fmincon is the best approach.  Notice that even if run with a random
% start, it is still faster than the other methods and often results in
% fewer function calls. Among the gradient-free methods, we can see that
% the order of preference in terms of function calls is PS, SA, then GA.
% And if we considered time to solution, it would be SA, PS, then GA. (Note
% that I haven't used the vectorize option in these solvers).
% 
% Notice that we also can’t necessarily relate function evaluation with
% time to find the solution when using different solvers.  This becomes
% even more of a consideration when you use parallel computing, which I’ll
% discuss in a minute.  In general, my recommendation is that if you think
% you have a nonlinear problem like this, or one that you can’t use a
% gradient based solver on, I’d start with PS as it is typically the most
% efficient.  I’d then use simulated annealing, as it is a compromise
% between PS and GA and can also handle customized data types, such as
% integers.
%
% However, GA is likely your best choice if you have custom data types or
% you have knowledge of the problem that can lend itself to modifying GA’s
% population selection and modification  functions.  GA can be efficient if
% you take the time to set up these functions.  I've glossed over the
% complexity of GA in this demo, and you can tune it to perform better.
%% Using parallel computing with optimization
% Both Optimization Toolbox and Genetic Algorithm and Direct Search Toolbox
% have built-in support for parallel computing (requires Parallel Computing
% Toolbox).  The built-in support is in selected constrained nonlinear
% solvers in Optimization Toolbox and Genetic algorithms and Pattern Search
% solvers in GADS Toolbox.
%
% Let's see how it works.  First lets time our objective
tic
peaksObj(x0)
toc
%%
% Let's now run our multistart |fmincon| example and time it.
clear all, close all, clc
a = -3;
b = 3;
x1 = a + (b-a) * rand(2);
x2 = a + (b-a) * rand(2);
x0r = [x1; x2];

options = optimset('Algorithm','active-set','Display','none');
tic
for i = 1:length(x0r)
    xopt(i,:) = fmincon(@peaksObj,x0r(i,:),[],[],[],[],[],[],@peaksCon,options);
end
toc
%%
% And let's rerun it using the build in parallel computing capabilities.
% Note that to use parallel computing, I need to set an option and open up
% the matlabpool (start the additional labs). I'm running this on a dual
% core laptop, so I can use two processing cores (labs) at a time.
options = optimset(options, 'UseParallel','Always');
matlabpool open local 2
tic
for i = 1:length(x0r)
    xopt(i,:) = fmincon(@peaksObj,x0r(i,:),[],[],[],[],[],[],@peaksCon,options);
end
toc
%%
% Hey, what gives?  This doesn't provide me any benefit!
%
% True, the built-in support is only useful for objective function and/or
% nonlinear constraint functions that are computationally intensive.  What
% constitutes an expensive function relies a lot on your local network and
% hardware, but a general rule of thumb is that if your objective function
% executes on the order of a few hundred milliseconds, you'll probably
% benefit from using parallel computing.
%
% Let's tweak my problem so it is expensive.  I've added a line of code to
% peaksOjb that simply takes time to compute (if you're interested, it's in
% the file peaksObjExp).
%
% Time it again.
tic
peaksObjExp(x0r(1,:));
toc
%%
% We need to rerun the base case since we are using a different objective
% function.
options = optimset(options, 'UseParallel','Never');
tic
for i = 1:length(x0r)
    x = fmincon(@peaksObjExp,x0r(i,:),[],[],[],[],[],[],@peaksCon,options);
end
toc
%%
% And let's see if we get an improvement using my two cores.
options = optimset(options, 'UseParallel','Always');
tic
for i = 1:length(x0r)
    x = fmincon(@peaksObjExp,x0r(i,:),[],[],[],[],[],[],@peaksCon,options);
end
toc

%% 
% You can see that I was able to get roughly a 15% improvement.  This
% improvement wasn't bad considering I have one machine with two cores that
% share the same memory.  You'd probably be able to get a better result
% using more cores or more computers in a cluster.
%
% Let's now take a look at how you could do even better by exploiting the
% nature of this problem using parallel computing on the for loop.
%% Do-it-yourself parallel optimization
% The built-in parallel computing support in |fmincon| accelerates the
% estimation of gradients.  While this helped to speed up the solution
% time, keep in mind that each iteration is dependent upon the previous
% iteration's results.  In other words, there are limits to how much
% parallel computing can speed up a gradient based solver since it contains
% serial operations that can't be parallelized.
%
% This problem can benefit form using parallel computing in a different
% way.  Notice that the optimization is being performed four times (the for
% loop).  Each loop is independent of the others, so we can make use of
% parallel computing in this case, and the Parallel Computing Toolbox makes
% this easy with the |parfor| language construct.
options = optimset(options, 'UseParallel','Never');
tic
parfor i = 1:length(x0r)% Notice the for loop is changed to parfor
    x = fmincon(@peaksObjExp,x0r(i,:),[],[],[],[],[],[],@peaksCon,options);
end
toc

%%
% You can see that the solution time was reduced by more than 30%.  I've
% showed this method to illustrate that you may benefit from using
% parallel computing directly in your problems, rather than relying on our
% built-in solution in some cases.  If you do have an expensive objective
% function, you might be able to use parallel computing in the definition
% of your objective function to get better performance than using the
% built-in support.  Often the expense arises in the objective
% function evaluation, and if you can speed that up you'll see the
% benefit in the optimization solver performance.
clear all, close all, clc
matlabpool close

%% Scaling considerations
% I'd like to show you some additional considerations when working with
% optimzation.  When solving problems, you should pay attention to the
% magnitude of the decision variables and the objective functions.  If there
% is a wide range, this can impact how the solver proceeds to a
% solution (numerical issues).  Let's look at the case where I have two
% decision variables of different magnitudes.  The plot below shows an
% engine performance map.  The axes show engine speed in RPM and the intake
% pressure ratio.  The color contours are a measure of engine efficiency.
% Note how the magnitude of the engine speed is nearly 10,000 times that of
% the pressure ratio.  This can present a problem for an optimization
% solver.
load VEdata
VEPlot([3500 0.5]);
%%
% If you think about this magnitude difference in terms of how a gradient
% based solver searches using derivatives of the objective function with
% respect to the decision variables, we get a large difference in the
% change seen due to the magnitude of decision variables (i.e. df/dx1
% ~1/1000, df/dx2 ~1/0.1).  This large difference leads to numerical
% problems, and as we see in this example, the optimization solver doesn't
% progress towards the solution.
options = optimset('Disp','Iter','Algorithm','interior-point',...
                   'OutputFcn',@VEPlot);
objFun = @(x)-VEMap(x,RPM,Pratio,VE);
x0 = [3500 0.5];
lb = [1000 0 ];
ub = [6000 1.0];
tic
z = fmincon(objFun,x0,[],[],[],[],lb,ub,[],options)
toc

%% Redefine RPM to have same scale as Pratio
% It is good practice to define your optimization problems such that the
% magnitudes of the objective/constraint functions and decision variables
% are of similar magnitude.  In this example, I rescaled the decision
% variables to be of a similar magnitude.  I've also changed the objective
% function definition to accept the scaled decision variables.
%
% Notice that when the problem is scaled correctly, we get the right
% solution.
close(gcf)
options = optimset('Disp','Iter','Algorithm','interior-point',...
                   'OutputFcn',@(x,o,s)VEPlot([x(1)*10000 x(2)],o,s));
objFun = @(x)-VEMap([x(1)*10000 x(2)],RPM,Pratio,VE);
x0 = [3500/10000 0.5];
lb = [1000/10000 0 ];
ub = [6000/10000 1.0];
tic
z = fmincon(objFun,x0,[],[],[],[],lb,ub,[],options)
toc

%% Tolerances (TolX, TolFun, etc)
% Since we're on the topic of magnitudes, let's return to the Mt.
% Washington example.  Notice the scale of the axes, 1,000s of feet.  Also
% note that we are seeking the peak, which is also in 1,000s of feet.  If
% we run this problem with default settings, we'll get the right solution
% but it will take 68 iterations to do so.
clear all, close all
% Load data for white mountains
load mtWashington;   
% Load psproblem which have all required settings for pattern search
load psproblem;
% show iteration history
psproblem.options.Display = 'iter';
% Run the solver and show displays
patternsearch(psproblem);
%%
% If you take a look at the default settings for the tolerances on the
% objective function (TolFun) and on the decision variables (TolX) you'll
% note that it is 1e-6.  I don't need tolerances anywhere near this tight
% since my problem is in 1,000s of feet.  I doubt my elevation data is even
% this accurate.  This is unnecessary and you can significantly speed up
% the optimization solver performance by adjusting these parameters to be
% closer to the accuracy required for the problem.  In this case, I'll set
% the TolFun to 10 (feet), and rerun the solver.  We can see that the
% iterations have been reduced to 43, nearly a 40% reduction.
psproblem.options.TolFun
psproblem.options.TolX

% Adjust the tolerance
psproblem.options.TolFun = 10;

% Run the solver and show displays
patternsearch(psproblem);

%% End game
% This concludes the examples shown in the webinar (plus some additional
% ones).  I hope you enjoyed them and learned something new along the way.
