Genetic Algorithm and Direct Search Toolbox 2.4.2

Optimization of Non-smooth Objective Function

This is a demonstration of how to find a minimum of a non-smooth objective function using the GA and PATTERNSEARCH functions in the Genetic Algorithm and Direct Search Toolbox™. Traditional derivative-based optimization methods, like those found in the Optimization Toolbox™, 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.

Contents

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
reset(RandStream.getDefaultStream); % Reset the state of random number gener
ator

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

function [f, g] = nonSmoothFcn(x)
%NONSMOOTHFCN is a non-smooth objective function

%   Copyright 2005 The MathWorks, Inc.
%   $Revision: 1.1.6.1 $   $Date: 2005/06/21 19:22:11 $

for i = 1:size(x,1)
    if  x(i,1) < -7
        f(i) = (x(i,1))^2 + (x(i,2))^2 ;
    elseif x(i,1) < -3
        f(i) = -2*sin(x(i,1)) - (x(i,1)*x(i,2)^2)/10 + 15 ;
    elseif x(i,1) < 0
        f(i) = 0.5*x(i,1)^2 + 20 + abs(x(i,2))+ patho(x(i,:));
    elseif x(i,1) >= 0
        f(i) = .3*sqrt(x(i,1)) + 25 +abs(x(i,2)) + patho(x(i,:));
    end
end

%Calculate gradient
g = NaN;
if x(i,1) < -7
    g = 2*[x(i,1); x(i,2)];
elseif x(i,1) < -3
    g = [-2*cos(x(i,1))-(x(i,2)^2)/10; -x(i,1)*x(i,2)/5];
elseif x(i,1) < 0
    [fp,gp] = patho(x(i,:));
    if x(i,2) > 0
        g = [x(i,1)+gp(1); 1+gp(2)];
    elseif x(i,2) < 0
        g =  [x(i,1)+gp(1); -1+gp(2)];
    end
elseif x(i,1) >0
    [fp,gp] = patho(x(i,:));
    if x(i,2) > 0
        g = [.15/sqrt(x(i,1))+gp(1); 1+ gp(2)];
    elseif x(i,2) < 0
        g = [.15/sqrt(x(i,1))+gp(1); -1+ gp(2)];
    end
end

function [f,g] = patho(x)
Max = 500;
f = zeros(size(x,1),1);
g = zeros(size(x));
for k = 1:Max  %k
   arg = sin(pi*k^2*x)/(pi*k^2);
   f = f + sum(arg,2);
   g = g + cos(pi*k^2*x);
end

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].

showNonSmoothFcn(Objfcn,range);
set(gca,'CameraPosition',[-36.9991   62.6267  207.3622]);
set(gca,'CameraTarget',[0.1059   -1.8145   22.3668])
set(gca,'CameraViewAngle',6.0924)
% Plot of the starting point (not used by the GA solver)
plot3(X0(1),X0(2),Objfcn(X0),'or','MarkerSize',10,'MarkerFaceColor','r');
fig = gcf;

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 = gaoptimset('PlotFcns',@gaplotbestfun,'PlotInterval',5, ...
    'PopInitRange',[-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;
plot3(Xga(1),Xga(2),Fga,'vm','MarkerSize',10,'MarkerFaceColor','m');
hold off;
fig = gcf;

Optimization terminated: average change in the fitness value less than optio
ns.TolFun.

Xga =

   -4.5766    0.0216


Fga =

   13.0186

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 the 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 = gaoptimset(optionsGA,'Generations',15, 'PlotInterval',1,...
    'HybridFcn',{@fminunc,optimset('OutputFcn',@fminuncOut)});
[Xhybrid,Fhybrid] = ga(FitnessFcn,2,optHybrid);

Optimization terminated: maximum number of generations exceeded.

Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the function tolerance.

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 - Fhy
brid)]);

The norm of |Xga - Xhb| is  0.13746
The difference in function values Fga and Fhb is 0.018614

Plot the final solution

figure(fig);
hold on;
plot3(Xhybrid(1),Xhybrid(2),Fhybrid+1,'^c','MarkerSize',10,'MarkerFaceColor'
,'c');
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 = psoptimset('PlotFcns',@psplotbestf);
% Run pattern search solver
[Xps,Fps] = patternsearch(ObjectiveFunction,X0,[],[],[],[],[],[],optionsPS)
% Plot the final solution
figure(fig)
hold on;
plot3(Xps(1),Xps(2),Fps+1,'*y','MarkerSize',14,'MarkerFaceColor','y');
hold off;

Optimization terminated: mesh size less than options.TolMesh.

Xps =

   -4.7124         0


Fps =

   13.0000