GEATbx: Tutorial 3 Writing Objective Functions

• 3.1 Parametric optimization functions
• 3.2 Defining default values of the objective function
• 3.3 Optimization of dynamic systems
• 3.4 Remark

3 Writing Objective Functions

When using the toolbox the implementation of an objective function consumes most of the work. Inside this function the calculation of the objective values depending on the variables is performed. Here you must implement your problem! The kind of implementation determines how good the evolutionary algorithm can work on and solve the problem.

Included in the distributed version of the toolbox are many examples of objective functions. All included objective functions follow the naming convention obj*.m (see also Naming Convention, Section 5.1). These functions implement a broad class of parameter optimization problems (and a few ordering/scheduling problems). When using these functions as a template it will be easier to implement your own functions/problems. For an overview of the mathematical description see Examples of Objective Functions.

Consider the following tasks:

1. The objective function is called with a matrix with as many rows as individuals. Every row corresponds to one individual. The number of columns determines the dimension/the number of variables of the objective function (see also Data Structures of the GEATbx, Chapter 8).
2. Because of being called with many individuals the objective function calculates the same mathematical expressions more than once. This could be done in a for-loop. However, here is a highly recommended place for vectorization.
3. Beside calculating the objective values the objective function could be used for defining default values (boundaries: lower and upper bounds of the variables, a default dimension of the problem; additionally, a descriptive string (name of the objective function) for labeling plots and, if known, the minimal objective value; for multi-objective problems the number of objective values).

3.1 Parametric optimization functions

Let's finish the theory. Here is an example.

Consider implementing the simple quadratic function (sum of quadrate or bowl function; known as De Jong's function 1 - objfun1). This functions works on real value variables.

Fig. 3-1: Definition of an objective function

function objval = objfun1(x)
    objval=sum((x.*x),2);   % vectorized, thus fast
    % for i=1:size(x,1), objval(i)=sum(x(i,:).^2); end

objfun1 returns the objective values of all individuals and because of the vectorization it will be fast. The second line (commented) shows the implementation of this function unvectorized. Every individual will be computed separately. The result is the same, however, Matlab takes a considerably longer time.

A fully documented version of the above objective function is implemented in objfun1.

Other examples of objective functions (using real value variables, vectorized implementation and a definable number of dimensions) are:

• objfun1a (axis parallel hyper-ellipsoid)
• objfun1b (rotated hyper-ellipsoid)
• objfun2 (Rosenbrock's function)
• objfun6 (Rastrigin's function)
• objfun7 (Schwefel's function)
• objfun8 (Griewangk's function)
• objfun9 (sum of different power)
• objfun10 (Ackley's path)
• objfun11 (Langermann's function)
• objfun12 (Michalewicz's function)

These functions are very often used as a standard set of test functions for evaluation of the performance of different evolutionary algorithms.

3.2 Defining default values of the objective function

The next step is defining default values for boundaries (domain of variables), best objective value and a description for the function. During the work on the toolbox the following style developed and is implemented inside all provided objective functions: Call the objective function with one special individual only (first variable is NaN, the optional second defines the type of info requested, for example x = [NaN, 1]).

Fig. 3-2: Definition of special return values of an objective function

function objval = objfun1(x, option)
    [Nind,Nvar]=size(x);
    % if Chrom is [NaN xxx] define size of boundary-matrix and others    
    if all([Nind == 1, isnan(Chrom(1))]),    
       % If only NaN is provided    
       if length(Chrom) == 1, option = 1; else option = Chrom(2); end    
       % Default dimension of objective function    
       Dim = 20;    
       % return text of title for graphic output    
       if option == 2, ObjVal = ['DE JONGs function 1'];    
       % return value of global minimum    
       elseif option == 3, ObjVal = 0;    
       % define size of boundary-matrix and values    
       else    
          % lower and upper bound, identical for all n variables    
          ObjVal = repmat([-512; 512], [1 Dim]);    
       end    
    else    
        % Compute objective function    
    end

The line objfun1([NaN,1]) or objfun1([NaN]) will return the default boundaries (includes implicitly the dimension) of the objective function. This is used inside a number of functions to retrieve the default boundaries and implicitly the number of variables for a given objective function.

objfun='objfun1';
 VLUB = feval(objfun, [NaN, 1]);
 geamain2(objfun, [], VLUB);

Getting the descriptive name of the objective function (objfun1([NaN,2])) or the minimal objective value (objfun1([NaN,3])) is not used in the distributed version of the toolbox (because I can not guaranty that you write all your objective functions according to this scheme). However, when comparing different algorithms and defining a termination criterion against the minimal objective value it is very useful (and I use it quite often for my own work).

A fully documented version of the above objective function is implemented in objfun1.

3.3 Optimization of dynamic systems

Often the solution of a problem involves the simulation of a system or the call of other functions. Consider the optimization of the control vector of a double integrator (push cart system). An overview of this system is given in Examples of Objective Functions.

function objval=objdopi(Chrom,option)
    [Nind,Nvar]=size(Chrom);
    XINIT=[0;-1];XEND=[0;0];
    TSTART=0; TEND=1;TIMEVEC=linspace(TSTART,TEND,Nvar)';
    if Nind==0,
    % see above example or objdopi
    else
       STEP=abs((TEND-TSTART)/(Nvar-1)));
       for indrun=1:Nind
          control=[TIMEVEC [Chrom(indrun,:)]'];
          [t x]=rk23('simdopiv', [TSTART TEND], XINIT, ...
                     [1e-3,STEP,STEP],control);
          % Calculate objective function
          objval(indrun)=sum(abs(x(size(x,1),:)'-XEND))+ ...
                         trapz(t,Chrom(indrun,:).^2));
       end
    end

At the beginning of the calculation of the objective values problem specific parameters are defined (XINIT, XEND, TSTART, TEND, TIMEVEC). The direct calculation is done separately for every individual (rk23 is written for only one system at one time). Every individual is converted in the form required by rk23 (vector of time values in the first column, next column(s) contain control values at specified time). The simulation function is called with appropriate parameters, the state values are returned. (rk23 is a Matlab4 specific function. The calling and setup of the integration routine changed under Matlab5 - see help for sim and simset.)

The objective value is a combination of two terms:

• sum of all values of the individual (the needed energy/force to change the state of the system) and
• difference between reached and needed end value of the states.

The objective function is finished. However, there is still another function needed - the s-function for the simulation routine simdopiv. (for an introduction to writing s-functions see the Simulink reference guide or look at the provided s-functions of the GEATbx sim*.m).

function [sys,x0]=simdopiv(t,x,u,flag);
    if abs(flag) == 1
       sys(1,:) = u(1,:);
       sys(2,:) = x(1,:);
    elseif abs(flag)==0
       sys=[2,0,0,1,0,0]; x0 = [0; -1];
    end

Let's go one step further. The toolbox provides the possibility to pass up to 10 parameters to the objective function. In the above example it could be useful to have the chance of changing TSTART, TEND, XINIT, XEND from outside and thus optimizing different situations of the double integrator system. The beginning of the function would be changed to:

function objval=objdopi(Chrom,option,TSTART,TEND,XINIT,XEND)
    [Nind,Nvar]=size(Chrom);
    if nargin<3, TSTART=0; end
    if nargin<4, TEND=1; end
    if nargin<5, XINIT=[0;-1]; end
    if nargin<6, XEND=[0;0]; end

The problem specific parameters are checked and if not provided set to default values. Thus, if necessary the parameters can be passed to the function or default values will be used. For optimizing objdopi the script would be (implemented in demodopi):

objfun='objdopi';
 TSTART=0; TEND=2; XINIT=[0;-2]; XEND=[0;0];
 tbxmpga(objfun, [], [], 0, TSTART, TEND, XINIT, XEND)

An even more advanced parameter checking could be done with (example):

   if nargin<3, TSTART=[]; end
    if isempty(TSTART), TSTART=0; end
    if nargin<4, TEND=[]; end
    if isempty(TEND), TEND=1; end

The parameter passing mechanism offers the possibility of getting multiple solutions at once without editing the objective function:

   for i=1:10,
       XINIT=[0;-i];
       geamain2(objfun,GeaOpt,VLUB,[],0,[],[],XINIT)
    end

Undefined parameters (TSTART, TEND, XEND) will be set to default values - really useful in day to day work.

One problem remains - vectorization. Many simulation problems are not vectorizable and thus slow. The evolutionary algorithm spends most of the time calculating the objective values. For simulation functions there are actually two problems: The Matlab-provided integration functions (rk23, rk45 or sim) are not vectorized. And it's often difficult to vectorize the s-function (see sim*.m). However, for this example both of these problems are solved. The s-function simdopiv is vectorized and together with the GEATbx a vectorized integration routine, intrk4, is provided. For more information see the documentation of intrk4 and method=12 inside objdopi.

If the objective function computes quite a few temporary results that are not part of the objective value it is good style to return them as additional output parameters.

function [objval, t, x]=objdopi(...)

What is the benefit? This opens the possibility of writing problem specific result plotting or special result computing routines. An example is implemented for the double integrator, which serves at the same time as an example for these features.

1. special initialization function (initdopi), [provides problem specific knowledge during initialization of the population].
2. special state plot function (plotdopi), [plot special results for best individual during optimization, for instance state and output variables of simulation].

For using these features the appropriate GeaOpt parameters (Special.InitDo, Special.InitFunction, Output.StatePlotInterval, Output.StatePlotFunction) must be defined. When calling the evolutionary algorithm Special.InitDo should be set to 1 and Output.StatePlotInterval should be set to 1 or greater. The corresponding *Function options contain the manes of the special initialization or state plot functions. The concept of setting the parameter in the options structure for using (or not) the special feature and defining the name of the special function is very flexible. For every objective function a special function for initialization or state plot can be defined (every problem needs it's own initialization or state plot function). Via parameter the use can be switched on or off. An example of all this is provided with the demo function demodopi, the objective function objdopi, the initialization function initdopi, and the state plot function plotdopi.

3.4 Remark

Two things should be stated at the end:

• Before starting to write own functions have a look to the provided example functions. It's much easier to use one of them as a template than starting from scratch. It is not necessary that everybody solves the same problems again and again.
• If good examples are known, which could be used for introducing a whole class of problems, please send them to the me (support@geatbx.com). If appropriate they will be included into the documentation of the toolbox and this tutorial as an example function.