# 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™, 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 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
```

function [f, g] = nonSmoothFcn(x) %NONSMOOTHFCN is a non-smooth objective function % Copyright 2005 The MathWorks, Inc. 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].

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;

Optimization terminated: average change in the fitness value less than options.FunctionTolerance. Xga = -4.6555 -0.0164 Fga = 13.0034

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);

Optimization terminated: maximum number of generations exceeded. Local minimum found. Optimization completed because the size of the gradient is less than the selected value of the optimality 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 - Fhybrid)]);

The norm of |Xga - Xhb| is 0.059157 The difference in function values Fga and Fhb is 0.0033556

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;

Optimization terminated: mesh size less than options.MeshTolerance. Xps = -4.7124 0 Fps = 13.0000