GEATbx: Tutorial 4 Writing Objective Functions

• 4.1 Parametric optimization functions
• 4.2 Optimization of dynamic systems
• 4.3 Remark

4 Writing Objective Functions

When using the toolbox the implementation of an objective function consumes most of the work. Inside this function all the problem specific variables are defined. The kind of implementation determines how good the evolutionary algorithm can work on and solve the problem.

One possible way for implementing the objective function is to provide the function as matlab expression directly (see Chapter 2). However, this is only applicable for small problems.

Included in the distributed version of the toolbox are many examples of objective functions. All included objective functions follow the naming convention obj*.m. These functions implement a broad class of parameter optimization problems. When using these functions as a template it will be easier to implement 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.
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 for lower and upper bound and a default dimension of the problem. Additionally, the example functions provide a descriptive string for labeling plots and, if known, the minimal objective value.

4.1 Parametric optimization functions

Let's finish the theory. Here is a first 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.

function objval = objfun1(x)
    objval=sum((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.

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 no individuals (the matrix with the individuals is empty). A second parameter determines which default parameter is to return.

function objval = objfun1(x, option)
    [Nind,Nvar]=size(x);
    if Nind== 0,
       if option==2, objval='DE JONG function 1';
       elseif option==3, objval=0;
       else Dim=20; objval=repmat([-512;+512],[1 Dim]);
    end
    else
       objval=sum((x.^2)')';
    end

The line objfun1([],1) will return the default boundaries (including the dimension) of the objective function. This is used inside, for instance, scrfun1 for retrieving default boundaries.

objfun='objfun1';
 bounds=feval(objfun, [], 1);
 VLB=bounds(1,:); VUB=bounds(2,:);
 tbxmpga(objfun, [], VLB, VUB);

Getting the descriptive name of the objective function (objfun1([],2)) or the minimal objective value (objfun1([],3)) is not used in the distributed version of the toolbox. However, when comparing different algorithms and defining a termination criterion against the minimal objective value it would be useful.

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.

4.2 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 scrdopi):

objfun='objdopi';
 bounds=feval(objfun,[],1);
 VLB=bounds(1,:); VUB=bounds(2,:);
 TSTART=0; TEND=2; XINIT=[0;-2]; XEND=[0;0];
 tbxmpga(objfun, [], VLB, VUB, 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];
       tbxmpga(objfun,[],VLB,VUB,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 advanced problem specific features, i.e.:

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 GOPTIONS parameters (GOPTIONS(29) and GOPTIONS(30)) and the corresponding name of the special function must be defined in global variables. When calling the evolutionary algorithm GOPTIONS(29) should be set to 1 and GOPTIONS(30) should be set to 1 or greater. 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 in global variables 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 start script scrdopi, the objective function objdopi, the initialization function initdopi, and the state plot function plotdopi.

4.3 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 author. If appropriate they will be included into the documentation of the toolbox and this tutorial as an example function.