www.gusucode.com > globaloptim 案例源码程序 matlab代码 > globaloptim/nonSmoothOpt.m
%% Optimization of Non-smooth Objective Function
% This example shows how to find a minimum of a non-smooth
% objective function using the |ga| and |patternsearch| functions in the
% Global Optimization Toolbox. Traditional derivative-based optimization
% methods, like those found in the Optimization Toolbox(TM), are fast
% and accurate for many types of optimization problems. These methods
% are designed to solve 'smooth', i.e., twice continuously differentiable,
% minimization problems, as they use derivatives to determine the direction
% of descent. While using derivatives makes these methods fast and
% accurate, they often are not effective when problems lack smoothness,
% e.g., problems with discontinuous, non-differentiable, or stochastic
% objective functions. When faced with solving such non-smooth problems,
% methods like the genetic algorithm or the more recently developed pattern
% search methods are effective alternatives.
% Copyright 2005-2015 The MathWorks, Inc.
%% Initialization
Objfcn = @nonSmoothFcn; % Handle to the objective function
X0 = [2 -2]; % Starting point
range = [-6 6;-6 6]; % Range used to plot the objective function
rng(0,'twister') % Reset the global random number generator
%% Non-smooth Objective Function
% To help visualize the problem and results, we have chosen a problem with
% only two variables, but the pattern search and genetic algorithms are
% certainly not limited to such small problems. We can view the code for
% this objective function.
type nonSmoothFcn.m
%%
% The objective function of our sample problem is a piece-wise continuous
% function, that is, it has smooth regions separated by discontinuities,
% with one exceptional region that is non-differentiable almost everywhere.
% We plot the objective function over range = [-6 6;-6 6].
fig = figure('Color','w');
[ah,sh] = showNonSmoothFcn(Objfcn,range);
ah.CameraPosition = [-36.9991 62.6267 207.3622];
ah.CameraTarget = [0.1059 -1.8145 22.3668];
ah.CameraViewAngle = 6.0924;
% Plot of the starting point (not used by the |ga| solver)
ph = [];
ph(1) = plot3(X0(1),X0(2),Objfcn(X0),'ob','MarkerSize',10,'MarkerFaceColor','b');
% Filter out FaceColor warning
ws = warning('off','MATLAB:legend:UnsupportedFaceColor');
oc = onCleanup(@()warning(ws));
% Add legend to plot
legendLabels = {'Non-differentiable regions','Smooth regions','Start point'};
lh = legend([sh,ph],legendLabels,'Location','SouthWest');
lp = lh.Position;
lh.Position = [0.005 0.005 lp(3) lp(4)];
%% Minimization Using the Genetic Algorithm
% The motivation for the genetic algorithm is evolutionary biology and
% genetics, mainly Darwin's theory of survival of the fittest. The Genetic
% Algorithm (|ga|) works on a population using a set of operators that are
% applied to the population. A population is a set of points in the design
% space. The initial population is generated randomly by default. The next
% generation of the population is determined using the fitness of the
% individuals in the current generation. The genetic algorithm does not use
% derivatives to detect descent in its minimization steps. Therefore, it is
% a good choice for problems such as this non-differentiable problem, as
% well as for discontinuous or stochastic problems.
%%
% To start, we will use the Genetic Algorithm, |ga|, alone to find the
% minimum of the objective function (also called a fitness function). We
% need to supply |ga| with a function handle to the fitness function
% nonSmoothFcn.m. Also, |ga| needs to know the how many variables are in the
% problem, which is two for this function.
FitnessFcn = @nonSmoothFcn;
numberOfVariables = 2;
%%
% Some plot functions are selected to monitor the performance of the solver.
optionsGA = optimoptions(@ga,'PlotFcn',@gaplotbestfun, ...
'InitialPopulationRange',[-5;5]);
% We run |ga| with the options 'optionsGA' as the tenth argument.
[Xga,Fga] = ga(FitnessFcn,numberOfVariables,[],[],[],[],[],[],[],optionsGA)
% Plot the final solution
figure(fig)
hold on;
ph(2) = plot3(Xga(1),Xga(2),Fga,'vm','MarkerSize',10,'MarkerFaceColor','m');
legendLabels = [legendLabels, '|ga| solution'];
lh = legend([sh,ph],legendLabels,'Location','SouthWest');
lp = lh.Position;
lh.Position = [0.005 0.005 lp(3) lp(4)];
hold off;
%%
% The optimum is at x* = (-4.7124, 0.0). |ga| found the point
% (-4.7775,0.0481) near the optimum, but could not get closer with the
% default stopping criteria. By changing the stopping criteria, we might
% find a more accurate solution, but it may take many more function
% evaluations to reach x* = (-4.7124, 0.0). Instead, we can use a more
% efficient local search that starts where |ga| left off. The hybrid function
% field in |ga| provides this feature automatically.
%% Minimization Using a Hybrid Function
% We will use a hybrid function to solve the optimization problem, i.e.,
% when |ga| stops this hybrid function will start from the final point
% returned by |ga|. Our choices are |fminsearch|, |patternsearch|, or |fminunc|.
% Since this optimization example is smooth near the optimizer, we can use
% the |fminunc| function from Optimization Toolbox as our hybrid function
% as this is likely to be the most efficient. Since |fminunc| has its own
% options structure, we provide it as an additional argument when
% specifying the hybrid function.
%%
% Run |ga|-|fminunc| hybrid
optHybrid = optimoptions(optionsGA,'MaxGenerations',15, ...
'HybridFcn',{@fminunc,optimoptions(@fminunc,'OutputFcn',@fminuncOut)});
[Xhybrid,Fhybrid] = ga(FitnessFcn,2,optHybrid);
%%
% The plot shows the best value of the population in every generation. The
% best value found by |ga| when it terminated is also shown in the plot
% title. When |ga| terminated, |fminunc| (the hybrid function) was
% automatically called with the best point found by |ga| so far. The solution
% |Xhybrid| with function value |Fhybrid| is the result of using |ga| and
% |fminunc| together. As shown here, using the hybrid function can
% efficiently improve the accuracy of the solution. Also, the total number
% of generations was reduced from 100 to 18 by using |fminunc| solver as
% hybrid a function. The improvement in |X| and function value is calculated
% below.
disp(['The norm of |Xga - Xhb| is ', num2str(norm(Xga-Xhybrid))]);
disp(['The difference in function values Fga and Fhb is ', num2str(Fga - Fhybrid)]);
%%
% Plot the final solution
figure(fig);
hold on;
ph(3) = plot3(Xhybrid(1),Xhybrid(2),Fhybrid+1,'^g','MarkerSize',10,'MarkerFaceColor','g');
legendLabels = [legendLabels, '|ga|-|fminunc| solution'];
lh = legend([sh,ph],legendLabels,'Location','SouthWest');
lp = lh.Position;
lh.Position = [0.005 0.005 lp(3) lp(4)];
hold off;
%% Minimization Using the Pattern Search Algorithm
% Pattern search, although much less well known, is an attractive
% alternative to the genetic algorithm as it is often computationally less
% expensive and can minimize the same types of functions.
%%
% Pattern search operates by searching a set of points called a pattern,
% which expands or shrinks depending on whether any point within the
% pattern has a lower objective function value than the current point. The
% search stops after a minimum pattern size is reached. Like the genetic
% algorithm, the pattern search algorithm does not use derivatives to
% determine descent, and so works well on non-differentiable, stochastic,
% and discontinuous objective functions. And similar to the genetic
% algorithm, pattern search is often very effective at finding a global
% minimum because of the nature of its search.
%%
% To minimize our objective function using the |patternsearch| function,
% we need to pass in a function handle to the objective function as well
% as specifying a start point as the second argument.
ObjectiveFunction = @nonSmoothFcn;
X0 = [2 -2]; % Starting point
% Some plot functions are selected to monitor the performance of the solver.
optionsPS = optimoptions(@patternsearch,'PlotFcn',@psplotbestf);
% Run pattern search solver
[Xps,Fps] = patternsearch(ObjectiveFunction,X0,[],[],[],[],[],[],optionsPS)
% Plot the final solution
figure(fig)
hold on;
ph(4) = plot3(Xps(1),Xps(2),Fps+1,'or','MarkerSize',4,'MarkerFaceColor','r');
legendLabels = [legendLabels, 'Pattern Search solution'];
lh = legend([sh,ph],legendLabels,'Location','SouthWest');
lp = lh.Position;
lh.Position = [0.005 0.005 lp(3) lp(4)];
hold off;