Fresa, a plant propagation algorithm

Author	Professor Eric S Fraga (e.fraga@ucl.ac.uk)
Last revision	2019-07-11 12:50:32

1 Overview and introduction

Update: Fresa now works with Julia version 1.0.0.

Fresa is an implementation of a plant propagation algorithm in Julia programming language. Our original version was called Strawberry, written in Octave (and usable in MATLAB). Please see the documentation of Strawberry for a longer description of the history of these codes, the conditions on which you may use these codes and which articles to cite. The move to Julia is motivated by the desire to use open source and free (as in beer) software only while also ensuring that the tools are computationally efficient. Octave has not kept up with developments in just in time compilation and, although has some features not available in MATLAB, it is not a modern language. Julia provides many of the capabilities of Octave while also being computationally efficient and providing many more modern features such as built-in parallelism and multi-dispatch, both of which are used by Fresa for efficiency and to make the search procedure generic with respect to data types used to represent points in the search space.

All code (either version) is copyright © 2019, Eric S Fraga, all rights reserved. Permission is given to use this code for evaluation purposes. The code is made available with no warranty or guarantee of any kind. Use at own risk.

Please let the author know if you download the code. Further, if either Fresa or Strawberry codes are used and this use leads to publications of any type, please do cite the publications listed in the first paragraph of the documentation for Strawberry. Feedback, including bug reports, is most welcome.

1.1 Installation

This document has been written using org mode in the Emacs text editor. Org allows for literate programming and uses tangling to generate the actual source files for the code. The code, tangled from this file, can be found here: Fresa.jl . To use this code, you could add this to your ~/.julia/config/startup.jl file (for Linux; for MS Windows, check the Julia documentation):

push!(LOAD_PATH,"/directory/where/Fresa/was/saved")

Alternatively, you can include the code directly before using it:

include("Fresa.jl")
using Fresa

The test examples that are described below all assume that the load path has been set as shown in the first code snippet above, i.e. in the configuration startup file.

Although Fresa has no external dependencies, the test problems below may require one or more of the following packages: BenchmarkTools, PyPlot, Profile.

1.2 Version information

Major version log:

July 2019: The returned values for Fresa.solve in the single objective case have changed. Instead of separate returned values for the decision variables, the objective function value, etc., a single Fresa.Point value is returned for the best point found, along with the full final population as an array of Fresa.point values.
June 2019: The calling interface for using the Fresa.solve method has changed. Specifically, when the search space is defined by data structures that are not a vector of Float64 values, the user must create a Fresa.neighbour function definition for the specific data structure type.
September 2017: moved to an object representation for points in the search space and allowed for parallel evaluation of the objective function when multiple processors are available.
November 2016: first Julia plant propagation algorithm implementation.

A list summary of recent change history is given below.

2 Fresa – The code and documentation

The Fresa method is a population based evolutionary algorithm which mimics the propagation used by plants. Throughout the module, the population object is an array of Point objects. Each point is a point in a search space, the objective function values for this point and a feasibility indication with g≤0 feasible and g>0 infeasible. See the documentation for the solve method below for more details on the data structures used and expected.

2.1 start of module and dependencies

Fresa depends on a number of packages that should be available in any Julia installation. These are packages for displaying output and using parallel computing capabilities when available on the actual hardware.

module Fresa
version = "[2019-07-11 12:50]"
using Dates
using Distributed
using Printf
function __init__()
    if myid() == 1
        println("# -*- mode: org; eval: (org-content 3); -*-")
        println("* Fresa PPA version $version")
    end
end

2.2 types

2.2.1 Point

Fresa uses one type, Point, which is a point in the search space. It includes the variable x, of indeterminate type to allow for a wide range of applications (e.g. integer versus real values), z, the value of the objective function, and g, the constraint violation (feasible with ≤0 and infeasible otherwise). An instance of a point is defined by the variable in the search space, the objective function used to evaluate the point, the ancestor of this point (see below), and optional parameters to pass to that function.

As Fresa is an evolutionary procedure, every point in the search space considered will be the descendent of a previously considered point. The sole exception is the initial starting point given by the procedure which invokes Fresa. The link between points is through a backward chain defined by the ancestor entry. This is not used by Fresa itself directly but provides extra meta-information that could be useful for post-optimization analysis, e.g. to see how effective the balance between exploration and exploitation may be for the given search parameter values.

"""

Point (`x`) in the search space along with objective function values
(`z[]`) and feasbility indication (`g`).  The type of `x` is problem
specific.  `z[]` and `g` hold `Float64` values.  `g` should be of
length 1.

"""
struct Point
    x :: Any                    # decision point
    z :: Array{Float64}         # objective function values
    g :: Float64                # constraint violation
    ancestor :: Union{Point,Nothing} # the parent of this point
    function Point(x, f, parameters, ancestor)
        if typeof(parameters) == Nothing
            (z, g) = f(x)
        else
            (z, g) = f(x,parameters)
        end
        # println("typeof(g) = $(typeof(g)) with length $(length(g))")
        if typeof(g) != Float64
            # println("typeof(g)==Array{Float64,1} = $(typeof(g)==Array{Float64,1})")
            if typeof(g) == Array{Float64,1} && length(g) == 1
                g = g[1]
            elseif typeof(g) == Int
                g = Float64(g)
            else
                error("Constraint g must be a single Float64 value.  Not $(typeof(g))")
            end
        end
        if typeof(z) == Array{Float64} || typeof(z) == Array{Float64,1}
            new(x,z,g,ancestor)
        elseif typeof(z) == Float64
            new(x,[z],g,ancestor)
        else
            error("Objective function values should be Float64 values, not $(typeof(z))")
        end
    end
end

Customise how a Point is displayed:

import Base
Base.show(io::IO, p::Fresa.Point) = print(io, "f(", p.x, ")=", p.z, " g=", p.g)
# and also an array of points
function Base.show(io::IO, p::Array{Point,1})
    np = length(p)
    if np > 0
        nz = length(p[1].z)
        println(io, "|-")
        print(io,"| x |")
        for i=1:nz
            print(io," z$(i) |")
        end
        println(io, " g |")
        println(io,"|-")
        for i=1:length(p)
            print(io, "| ", p[i].x, " |")
            for j=1:nz
                print(io," ", p[i].z[j], " |")
            end
            print(io, " ", p[i].g, " |\n")
        end
        println(io,"|-")
    else
        print(io,"empty")
    end
end

2.3 create a point

A trivial function that simply creates a new Point object. This is needed for the remotecall functionality when using parallel computing because the remotecall function has to be given a function and not just a constructor (for some obscure reason that means that a constructor is transformed to a conversion operation… don't ask me). The optional parameters argument is passed through: it is the Point creation method that checks to see if these are defined or not.

function createpoint(x,f,parameters,ancestor)
    return Point(x,f,parameters,ancestor)
end

2.4 fitness

The fitness function used depends on the number of objectives. For single criterion problems, the fitness is the objective function values normalised and reversed so that the minimum, i.e. the best solution, has a fitness of close to 1 and the worst a fitness close to 0. For multi-criteria problems, a Hadamard product of individual criteria rankings is used to create a fitness value Fraga & Amusat, 2016.

This function uses a helper function, defined below, to assign a fitness to a vector of objective function values.

function fitness(pop, fitnesstype)
    l = length(pop)
    indexfeasible = (1:l)[map(p->p.g,pop) .<= 0]
    indexinfeasible = (1:l)[map(p->p.g,pop) .> 0]
    # println("There are $(length(indexfeasible)) feasible entries and $(length(indexinfeasible)) infeasible")
    fit = zeros(l)
    factor = 1              # for placement in fitness interval (0,1)
    if length(indexfeasible) > 0
        feasible = view(pop,indexfeasible)
        # use objective function value(s) for ranking
        feasiblefit = vectorfitness(map(p->p.z,feasible), fitnesstype)
        if length(indexinfeasible) > 0
            feasiblefit = feasiblefit./2 .+ 0.5 # upper half of fitness interval
            factor = 2                        # have both feasible & infeasible
        end
        fit[indexfeasible] = (feasiblefit.+factor.-1)./factor
    end
    if length(indexinfeasible) > 0
        # squeeze infeasible fitness values into (0,0.5) or (0,1) depending
        # on factor, i.e. whether there are any feasible solutions as well or not
        infeasible = view(pop,indexinfeasible)
        # use constraint violation for ranking as objective function values
        # may not make any sense given that points are infeasible
        fit[indexinfeasible] = vectorfitness(map(p->p.g, infeasible), fitnesstype) / factor;
    end
    fit
end

The helper function works with a single vector of objective function values which may consist of single or multiple objectives.

"""
For single objective problems, the fitness is simply the normalised
objective function value.

For multi-objective cases, there are two types of fitness ranking that
we are interested in.  The first is based on the Hadamard product of
the rank of each member of population accoring to each criterion.  The
second is based on a weighted Borda ranking.
"""
function vectorfitness(v,fitnesstype)
    # determine number of objectives (or pseudo-objectives) to consider in
    # ranking
    l = length(v)
    if l == 1
        # no point in doing much as there is only one solution
        fit = [0.5]
    else
        m = length(v[1])
        # println("VF: v=$v")
        # println("  : of size $(size(v))")
        if m == 1                   # single objective 
            v = [v[i][1] for i=1:l]
            s = sortperm(v)
            zmin = v[s[1]]
            zmax = v[s[l]]
            if abs(zmax-zmin) < eps()
                fit = 0.5*ones(l)
            else
                # avoid extreme 0,1 values
                fit = tanh.((zmax .- v) ./ (zmax .- zmin) .- 0.5).+0.5
            end
        else                  # multi-objective
            rank = ones(m,l); #rank of each solution for each objective function 
            if fitnesstype == :hadamard
                for i=1:m
                    rank[i,sortperm([v[j][i] for j=1:l])] = 1:l;
                end
                # hadamard product of ranks
                fitness = map(x->prod(x), rank[:,i] for i=1:l)
                # normalise and reverse meaning (1=best, 0=worst) while avoiding
                # extreme 0,1 values using the hyperbolic tangent
                fit = tanh.(0.5 .- fitness ./ maximum(fitness)) .+ 0.5
            elseif fitnesstype == :borda
                for i=1:m
                    rank[i,sortperm([v[j][i] for j=1:l])] = l:-1:1;
                end
                # hadamard product of ranks
                fitness = map(x->sum(x), rank[:,i] for i=1:l)
                # normalise (1=best, 0=worst) while avoiding
                # extreme 0,1 values using the hyperbolic tangent
                if (maximum(fitness)-minimum(fitness)) > eps()
                    fit = tanh.((fitness .- minimum(fitness)) / (maximum(fitness)-minimum(fitness)) .- 0.5) .+ 0.5
                else
                    fit = 0.5*ones(l)
                end
            elseif fitnesstype == :nondominated
                # similar to that used by NSGA-II (Deb 2000)
                fitness = zeros(l)
                maxl = assigndominancefitness!(fitness,v,1)
                # println("Resulting fitness: $fitness")
                fit = tanh.((maxl.-fitness)./maxl .- 0.5) .+ 0.5
                # println(":  scaled fitness: $fit")
            else
                throw(ArgumentError("Type of fitness evaluation must be either :borda, :nondominated, or :hadamard, not $(repr(fitnesstype))."))
            end
        end
    end
    # println("VF: fit=$fit")
    fit
end

The following function is used by the vector fitness evaluation to recurse through the levels of non-dominance to assign fitness based on those levels.

function assigndominancefitness!(f,v,l)
    # assign value l to all members of v which dominate rest and then
    # recurse on those which are dominated
    (p, d) = pareto(v)
    # println("Assigning fitness $l to $p")
    f[p] .= l
    if !isempty(d)
        assigndominancefitness!(view(f,d),v[d],l+1)
    else
        l
    end
end

2.5 neighbour – generate random point

A random solution is generated with a distance from the original point being inversely proportional, in a stochastic sense, to the fitness of the solution. The new point is possibly adjusted to ensure it lies within the domain defined by the lower and upper bounds. The final argument is the fitness vector with values between 0 and 1, 1 being the most fit and 0 the least fit.

Fresa comes with two default methods for generating neighbouring solutions. The first is for a search space defined by vectors of Float64 values:

function neighbour(x :: Array{Float64,1},
                   a :: Array{Float64,1},
                   b :: Array{Float64,1},
                   f :: Float64
                   ) :: Array{Float64,1}
    xnew = x .+ (1.0 .- f) .* 2(rand(length(x)).-0.5) .* (b.-a)
    xnew[xnew.<a] = a[xnew.<a];
    xnew[xnew.>b] = b[xnew.>b];
    return xnew
end

There is also a version that expects single valued Float64 arguments.

function neighbour(x :: Float64,
                   a :: Float64,
                   b :: Float64,
                   f :: Float64
                   ) :: Float64
    # allow movements both up and down
    # in the domain for this variable
    newx = x + (b-a)*(2*rand()-1)/2.0 * (1-f)
    if newx < a
        newx = a
    elseif newx > b
        newx = b
    end
    newx
end

Should other decision point types be required, e.g. mixed-integer or domain specific data structures, the Fresa.neighbour function with parameters of the specific type will need to be defined. See the mixed integer nonlinear example below for an example of a simple mixed-integer case.

2.6 pareto – set of non-dominated points

Select a set consisting of those solutions in a population that are not dominated. This only applies to multi-objective optimisation; for a single criterion problem, the solution with minimum objective function value would be selected. This function is used only for returning the set of non-dominated solutions at the end of the solution procedure for multi-objective problems. It could be used for an alternative fitness function, a la Srinivas et al. (N Srinivas & K Deb (1995), Evolutionary Computation 2:221-248)

# indices of non-dominated and dominated points from the population of
# Point objects
function pareto(pop::Array{Point,1})
    l = length(pop)
    indexfeasible = (1:l)[map(p->p.g,pop) .<= 0]
    indexinfeasible = (1:l)[map(p->p.g,pop) .> 0]
    if length(indexfeasible) > 0
        subset = view(pop,indexfeasible)
        indices = indexfeasible
    else
        println(": Fresa.pareto warning: no feasible solutions.  Pareto set meaningless?")
        subset = pop
        indices = 1:l
    end
    z = map(p->p.z, subset)
    # use function below to return indices of non-dominated and
    # dominated from objective function values alone in the subset of
    # feasible solutions
    (p, d) = pareto(z)
    (indices[p], indices[d])
end

# set of non-dominated (and dominated) points from array of objective
# function values alone.
function pareto(z::Array{Array{Float64,1},1})
    l = length(z)
    p = Int[]                 # indices of pareto members in full population
    d = Int[]                 # indices for dominated members
    for i in 1:l
        dominated = false
        for j in 1:l
            if i != j
                if all(z[i] .>= z[j]) && any(z[i] .> z[j])
                    # println("$i dominated by $j")
                    # println("$(z[:,i]) >= $(z[:,j])")
                    dominated = true;
                    break;
                end
            end
        end
        # println("member $(pop[i]) is dominated: $dominated")
        if dominated
            push!(d,i)          # dominated
        else
            push!(p,i)          # pareto, i.e. non-dominated
        end
    end
    (p, d)
end

2.7 prune - control population diversity

Due to the stochastic nature of the method and also the likely duplication of points when elitism is used, there is a need to prune the population. We wish to remove members that have objective function values that are too close to each other. The main difficulty is the definition of too close. We use a tolerance based on the range of values present in the population.

function prune(pop::Array{Point}, tolerance)
    npruned = 0
    z = map(p->p.z, pop)
    # println("typeof(z)=$(typeof(z))")
    l = length(z)
    # println("typeof(z[1])=$(typeof(z[1]))")
    n = length(z[1])
    zmin = zeros(n)
    zmax = zeros(n)
    for i=1:n
        row = [z[j][i] for j=1:l]
        zmin[i] = minimum(row)
        zmax[i] = maximum(row)
        if zmax[i] - zmin[i] < 100*eps()
            zmax[i] = zmin[i]+100*eps()
        end
    end
    pruned = [pop[1]]
    for i=2:l
        similar = false
        for j=1:length(pruned)
            if all(abs.(z[i]-pruned[j].z) .< tolerance*(zmax-zmin))
                similar = true;
                break;
            end
        end
        if !similar
            push!(pruned,pop[i])
        else
            npruned += 1
        end
    end
    (pruned, npruned)
end

2.8 select – choose a member of the population

Given a fitness, f, choose two solutions randomly and select the one with the better fitness. This is known as a tournament selection procedure of size 2. Other options are possible but not currently implemented.

function select(f)
    l = length(f)
    ind1 = rand(1:l)
    if ind1 == 0
        ind1 = 1
    end
    ind2 = rand(1:l)
    # println("Comparing $ind1 to $ind2")
    if f[ind1] > f[ind2]
        return ind1
    else
        return ind2
    end
end

2.9 solve – use the PPA to solve the optimisation problem

The function expects the objective function, f, an initial guess, x0, and lower, a, and upper, b, bounds. It returns the optimum, the objective function value(s) at this point, the constraint at that point and the whole population at the end. The actual return values and data structures depends on the number of criteria:

1: returns best point as a Fresa.Point object (which includes the decision variable values, the objective function value, and the constraint value) and also the full population;
>1: returns the set of non-dominated points (as an array including objective function values and constraint value) and the full population.

The objective function, f, should return two results: z, the objective function value(s) which must be of type Float64, single or array, and g, the constraint violation. If g≤0, the point is feasible; any value g>0 means an infeasible point. The value of g for infeasible points will be used to rank the fitness of the infeasible solution, with lower values being fitter.

The calling sequence for f is a point in the search space plus, optionally, the parameters defined in the call to solve.

x0 is the initial guess and can be of any type. a and b are lower and upper bounds and should be of types consistent with each other and x0.

If the decision vector is not an array of Float64, a type specific Fresa.neighbour function will need to be defined. The calling sequence for Fresa.neighbour is (x,a,b,fitness) where x, a, and b, should all be of the desired type and the function itself must also return an object of that type. The fitness will always be a Float64. See the mixed integer nonlinear example below for an example.

The fitnesstype is used for ranking members of a population for multi-objective problems. The default is to use a Hadamard product of the rank each solution has for each objective individually. The alternative, specifying fitnesstype=:borda uses a sum of the rank, i.e. a Borda count. The former tends to emphasise points near the extrema of the individual criteria while the latter is possibly better distributed but possibly at providing less emphasis on the Pareto points themselves.

"""
Solve an optimisation problem, defined as the minimization of the
values returned by the objective function, `f`.  `f` returns not only
the objective function values, an array of `Float64` values, but also
a measure of feasibility (≤0) or infeasibility (>0).  The problem is
solved using the Fresa algorithm.  `x0` is the initial guess and `a`
and `b` are *bounds* on the search space.

The return values for the solution of a single criterion problem are
the best point and the full population at the end of the search. 

For a multi-objective problem, the returned values are the set of
indices for the points within the full population (the second returned
value) approximating the *Pareto* front.

The population will consist of an array of `Fresa.Point` objects, each
of which will have the point in the search space, the objective
function value and the feasibility measure.

"""
function solve(f, x0, a, b;     # required arguments
               parameters = nothing, # allow parameters for objective function 
               archiveelite = false,  # save thinned out elite members
               elite = true,    # elitism by default
               fitnesstype = :hadamard, # how to rank solutions in multi-objective case
               ngen = 100,      # number of generations
               npop = 10,       # population size
               nrmax = 5,       # number of runners maximum
               ns = 100,        # number of stable solutions for stopping
               output = 5,      # how often to output information
               plotvectors = false, # generate output file for search plot
               tolerance = 0.001) # tolerance for similarity detection
    println("** solve $f $(orgtimestamp(now()))")
    tstart = time()
    p0 = Point(x0, f, parameters, nothing)
    nf = 1                      # number of function evaluations
    npruned = 0                 # number solutions pruned from population
    nz = length(p0.z)           # number of criteria
    pop = [p0];                 # create/initialise the population object
    if archiveelite
        archive = Point[]
    end
    println(": solving with ngen=$ngen npop=$npop nrmax=$nrmax ns=$ns")
    println(": elite=$elite archive elite=$archiveelite fitness type=$fitnesstype")
    if plotvectors
        plotvectorio = open("fresa-vectors-$(orgtimestamp(now())).data", create=true, write=true)
        println(": output of vectors for subsequent plotting")
    end
    # we use parallel computing if we have more than one processor
    parallel = nprocs() > 1
    println(": function evaluations performed ",
            parallel ? "in parallel with $(nprocs()) processors." : "sequentially.")
    println("*** initial population")
    println(pop)
    if parallel
        # will be used to collect results from worker processors
        results = Array{Future,1}(nprocs())
    end
    println("*** evolution")
    @printf("| %9s | %9s | %9s | %9s | %9s |", "gen", "npop",
            (elite && nz > 1) ? "pareto" : "nf", "pruned", "t (s)")
    for i in 1:nz
        @printf(" z%-8d |", i)
    end
    @printf("\n|-\n")
    # now evolve the population for a predetermined number of generations
    for gen in 1:ngen
        # evaluate fitness
        fit = fitness(pop, fitnesstype)
        # sort
        index = sortperm(fit)
        # and remember best which really only makes sense in single
        # criterion problems but is best in multi-objective case in
        # the ranking measure used by Fresa
        best = pop[index[end]]
        # if elitism is used
        if elite
            if nz > 1
                # elite set is whole pareto set unless it is too
                # big. Recall that the pareto function returns the set
                # of indices into the population
                wholepareto = pareto(pop)[1]
                if length(wholepareto) > ceil(npop/2)
                    newpop, removed = thinout(pop, fit, wholepareto, ceil(Int,npop/2))
                    if archiveelite
                        archive = prune(append!(archive, removed), tolerance)[1]
                        archive = archive[pareto(archive)[1]]
                    end
                else
                    newpop = pop[wholepareto]
                end
            else
                # elite set is single element only
                newpop = [best]
            end
            # if plotting vectors for the search, include elitism
            if plotvectors
                for p in newpop
                    write(plotvectorio, "$(gen-1) $(p.x)\n$gen $(p.x)\n\n")
                end
            end
        else
            newpop = Point[]
        end
        print(stderr, ": $gen np=$(length(newpop))",
              archiveelite ? " na=$(length(archive))" : "",
              " with most fit z=$(best.z)           \r")
        if gen%output == 0
            @printf("| %9d | %9d | %9d | %9d | %9.2f |", gen, length(fit),
                    (elite && nz > 1) ? length(newpop) : nf, npruned, time()-tstart)
            for i = 1:length(best.z)
                @printf(" %9g |", best.z[i])
            end
            println()
        end
        if parallel
            # create array to store all new points; we evaluate them
            # later hopefully in parallel.  Also keep track of the
            # points from which new points are derived to provide the
            # backward link through the evolution
            x = typeof(newpop[1].x)[]
            points = Point[]
        end
        # now loop through population, applying selection and then
        # generating neighbours
        l = length(pop)
        for i in 1:min(l,npop)
            s = select(fit)
            # println(": selection $i is $s")
            # println(": size of pop is $(size(pop))")
            selected = pop[s]
            if !elite
                # if no elitism, we ensure selected members remain in population
                push!(newpop, selected)
                if plotvectors
                    write(plotvectorio, "$(gen-1) $(selected.x)\n$gen $(selected.x)\n\n")
                end
            end
            # number of runners to generate, function of fitness
            nr = ceil(fit[s]*nrmax*rand())
            if nr < 1
                nr = 1
            end
            # println(": generating $nr runners")
            for r in 1:nr
                # create a neighbour, also function of fitness
                newx = neighbour(pop[s].x,a,b,fit[s])
                # for parallel evaluation, we store the neighbours and
                # evaluate them later; otherwise, we evaluate
                # immediately and save the resulting point
                if parallel
                    push!(x, newx)
                    push!(points, pop[s])
                else
                    push!(newpop, Point(newx, f, parameters, pop[s]))
                    if plotvectors
                        write(plotvectorio, "$(gen-1) $(pop[s].x)\n$gen $newx\n\n")
                    end
                    nf += 1
                end
            end
            # remove selected member from the original population so
            # it is not selected again
            splice!(fit, s)
            splice!(pop, s)
        end
        # if we are making use of parallel computing, we evaluate all
        # points generated in previous loop.  
        if parallel
            i = 0;
            while i < length(x)
                # issue remote evaluation call
                for j=1:nprocs()
                    if i+j <= length(x) 
                        results[j] = @spawn createpoint(x[i+j],f,parameters,
                                                        points[i+j])
                        nf += 1
                    end
                end
                # now wait for results
                for j=1:nprocs()
                    if i+j <= length(x)
                        push!(newpop, fetch(results[j]))
                    end
                end
                i += nprocs()
            end
        end
        # and finally remove any duplicate points in the new
        # population and make it the current population for the next
        # generation
        (pop, nn) = prune(newpop, tolerance)
        npruned += nn
    end
    println("*** Fresa run finished\n: nf=$nf npruned=$npruned", archiveelite ? " archived=$(length(archive))" : "")
    if plotvectors
        close(plotvectorio)
    end
    if nz == 1
        fit = fitness(pop, fitnesstype)
        index = sortperm(fit)
        best = pop[index[end]]
        return best, pop
    else
        return pareto(archiveelite ? append!(pop,archive) : pop)[1], pop
    end
end

2.10 thinout – make Pareto set smaller

If we use elitism, for multi-objective problems, we use the Pareto set as the elite set. However, this set may grow to be large, causing performance challenges as well as making the search less effective at exploration, essentially getting stuck in the local area defined by this elite set. Therefore, we need to sometimes thin out the Pareto set for its use as an elite set.

The arguments are the whole population, the fitness of the members, the indices in this population for the Pareto set and the number of elements to keep. We keep the most fit ones.

function thinout(pop, fit, pareto, n::Int)
    indices = sortperm(fit[pareto])
    return pop[pareto[indices[end-n+1:end]]], pop[pareto[indices[1:end-n]]]
end

2.11 utility functions

Some functions that are not necessary for Fresa but provide some useful features, especially output related.

2.11.1 org time stamp

function orgtimestamp(dt::DateTime)
    return @sprintf("[%d-%02d-%02d %02d:%02d]",
                    Dates.year(dt),
                    Dates.month(dt),
                    Dates.day(dt),
                    Dates.hour(dt),
                    Dates.minute(dt))
end

2.12 module end

end

3 Tests

The following are simple tests for either the Fresa optimiser or just individual functions in the module. You can cut and paste these codes into your own editor and run them.

3.1 a GAMS interface example

The GAMS modelling system is used by many to write and solve optimization problems and many different solvers are available, including both local and global optimizers. However, there are some problems for which the solvers may not be able to find good solutions. Fresa may provide a suitable alternative solver for such problems. However, one of the best features of GAMS is that the model can be represented purely by the equations without the need to determine an evaluation sequence for these equations given a decision vector. It is therefore desirable to consider using Fresa with GAMS models.

This example implements an objective function which invokes GAMS to solve the model given values for some decision variables. This interface to GAMS requires writing and reading from files so will not be appropriate for small models due to the overheads in file access.

The files for this example can be downloaded: Julia code and GAMS model.

3.1.1 the GAMS model

We use, as an example, problem 8.26 in "Engineering Optimization" by Reklaitis, Ravindran and Ragsdell (1983). This problem seeks to minimise the square of the decision variables while minimising a residual value res. We treat this as a multi-objective problem which cannot be done directly in GAMS. By looking at it as a multi-objective problem, we can gain insight into the trade-offs between the residual and the primary objective function.

$TITLE Test Problem 
$OFFDIGIT
$OFFSYMXREF 
$OFFSYMLIST 

VARIABLES X1, X2, X3, Z, res ; 
POSITIVE VARIABLES X1, X2, X3 ; 

EQUATIONS CON1, CON2, CON3, OBJ ;

CON1..  X2 - X3 =G= 0 ; 
CON2..  X1 - X3 =G= 0 ; 
CON3..  X1 - X2**2 + X1*X2 - 4 =E= res ;
OBJ..   Z =E= SQR(X1) + SQR(X2) + SQR(X3) ; 

* Upper bounds 
X1.UP = 5 ; 
X2.UP = 3 ; 
X3.UP = 3 ; 

* Initial point 
X1.L = 4 ; 
X2.L = 2 ; 
X3.L = 2 ; 

MODEL TEST / ALL / ; 

OPTION LIMROW = 0; 
OPTION LIMCOL = 0;

3.1.2 a multi-objective function with interface to GAMS

The objective function for Fresa takes the decision variables, x, and uses these to set the GAMS model variables X1, X2, and X3. After solving the GAMS model, the results, consisting of the objective function value Z and the residual, res, are output to a file for subsequent reading into the Julia code. The absolute value of the residual is used as a second criterion.

function fmo(x::Array{Float64,1})
    open("gamsexample.gms", "w") do f
        write(f, "\$include gamsdeclarations.gms\n")
        write(f, "X1.fx = $(x[1]); \n")
        write(f, "X2.fx = $(x[2]); \n")
        write(f, "X3.fx = $(x[3]); \n")
        write(f, "solve TEST using NLP minimizing Z; \n")
        write(f, "file fresa /'gamsoutput.txt'/ ;\n")
        write(f, "put fresa ;\n")
        write(f, "put z.l /;\n")
        write(f, "put res.l /;\n")
        write(f, "put TEST.modelstat /;\n")
    end
    # execute GAMS
    run( `/opt/gams/latest/gams gamsexample.gms` )
    # read in results
    z = [0.0; 0.0]
    g = 0.0;
    open("gamsoutput.txt", "r") do f
        lines = readlines(f)
        z[1] = parse(Float64, lines[1])
        z[2] = abs(parse(Float64, lines[2]))
        modelstat = parse(Float64, lines[3])
        if modelstat != 1 && modelstat != 5
            g = 1
        end
    end
    # return results
    ( z, g )
end

3.1.3 a single objective function with interface to GAMS

In this case, the value of the residual, in absolute sense, is a measure of feasibility. We have a single criterion, the value of the GAMS objective function.

function fsingle(x::Array{Float64,1})
    open("gamsexample.gms", "w") do f
        write(f, "\$include gamsdeclarations.gms\n")
        write(f, "X1.fx = $(x[1]); \n")
        write(f, "X2.fx = $(x[2]); \n")
        write(f, "X3.fx = $(x[3]); \n")
        write(f, "solve TEST using NLP minimizing Z; \n")
        write(f, "file fresa /'gamsoutput.txt'/ ;\n")
        write(f, "put fresa ;\n")
        write(f, "put z.l /;\n")
        write(f, "put res.l /;\n")
        write(f, "put TEST.modelstat /;\n")
    end
    # execute GAMS
    run( `/opt/gams/latest/gams gamsexample.gms` )
    # read in results
    z = 0.0
    g = 0.0
    open("gamsoutput.txt", "r") do f
        lines = readlines(f)
        z = parse(Float64, lines[1])
        g = abs(parse(Float64, lines[2]))
        modelstat = parse(Float64, lines[3])
        if modelstat != 1 && modelstat != 5
            g = 10 # penalty function
        end
    end
    # return results
    ( z, g )
end

3.1.4 solve the multi-objective problem using Fresa

using Fresa
a = [0.0;0.0;0.0]
b = [5.0;3.0;3.0]
x0 = [4.0;2.0;2.0]
pareto, population = Fresa.solve(fmo, x0, a, b;
                                 fitnesstype = :borda,
                                 ngen = 100)
println("Pareto front:")
println(population[pareto])

and plot out the resulting Pareto set in objective function space:

using PyPlot
z = [population[pareto[i]].z for i in 1:length(pareto)];
PyPlot.plot([z[i][1] for i=1:length(z)],
            [z[i][2] for i=1:length(z)],
            "ro")
PyPlot.savefig("gamsmo.pdf")

3.1.5 solve the single objective version

best, pop = Fresa.solve(fsingle, x0, a, b; ngen = 100)
println("Population: $pop")
println("Best: f($(best.x)) = $(best.z), $( best.g )")

3.2 mixed integer nonlinear example

The MINLP example comes from: Tapio Westerlund & Joakim Westerlund, GGPECP – An algorithm for solving non-convex MINLP problems by cutting plane and transformation techniques, Proceedings of ICHEAP-6, Pisa, June 2003. It has one real variable and one integer variable. The search region is non-convex, consisting of two disjoint domains.

The aims of this example are to test the use of a non-default neighbour function and the use of a problem-specific type for solutions, a mixed-integer type in this case.

This example is also used, for the moment, to test out the parallel implementation of Fresa. The important aspects are that Fresa as well as the MI type be available on all worker processes. This is not a good example in that the parallel version takes longer than the sequential version.

using Distributed
using Printf
@everywhere using Fresa
# define new type for mixed integer problems
# in general, this would be vectors of real and integer values
@everywhere struct MI
    x :: Float64
    y :: Int32
end
import Base
Base.show(io::IO, m::MI) = print(io, m.x, " ", m.y)
f = s -> (3s.y - 5s.x,
          max(2s.y + 3s.x - 24,
              3s.x - 2s.y - 8,
              2s.y^2 - 2*√s.y + 11s.y + 8s.x - 39 - 2*√s.x*s.y^2))
# bounds
a = MI(1.0, 1)
b = MI(6.0, 6)
# function to find a neighbouring solution for MI type decision points
function Fresa.neighbour(s :: MI,
                         a :: MI,
                         b :: MI,
                         f :: Float64) :: MI
    x = s.x + (b.x-a.x)*(1-f)*2*(rand()-0.5)
    x = x < a.x ? a.x : (x > b.x ? b.x : x)
    # for the integer variable, we move in one direction or the other
    # a random number of places depending on fitness
    positive = rand(Bool)
    r = rand()
    # @printf(": neighbour: f=%g r=%g\n", f, r)
    inc = ceil(f*r*(b.y-a.y)/2)
    # @printf(": neighbour: positive=%s inc=%d\n", positive, inc)
    y = s.y + (positive ? inc : -inc)
    y = y < a.y ? a.y : (y > b.y ? b.y : y)
    return MI(x,y)
end
best, pop = Fresa.solve(f, MI(1.0, 1), a, b; ngen=100)
println("Population: $pop")
println("Best: f($(best.x)) = $(best.z), $(best.g)")

Using the results obtained above, we use the linking information for all the points in the search space encountered to do some simple analysis of the search. The simple analysis consists of printing out the path of evolution that led to the final best solution:

point = best;
while point != Some(nothing) && typeof(point) != Nothing
    global point
    println("$(point.x)")
    point = point.ancestor
end

3.3 multi-objective test

using Fresa
nx = 2
a = zeros(nx)
b = ones(nx)
x = rand(nx)
f = x -> ( [sin(x[1]-x[2]); cos(x[1]+x[2])], 0)
pareto, population = Fresa.solve(f, x, a, b;
                                 #fitnesstype = :hadamard,
                                 #fitnesstype = :borda,
                                 fitnesstype = :nondominated,
                                 ngen=100,
                                 npop=10,
                                 plotvectors=true,
                                 tolerance=0.01)

println("Pareto front:")
println(population[pareto])
#using BenchmarkTools
#@benchmark

using PyPlot
z = [population[pareto[i]].z for i in 1:length(pareto)];
PyPlot.plot([z[i][1] for i=1:length(z)],
            [z[i][2] for i=1:length(z)],
            "ro")
PyPlot.savefig("x.pdf")

3.4 multi-objective test with 3 objectives

using Fresa
using Profile
nx = 5
a = zeros(nx)
b = ones(nx)
x = zeros(nx)
f = x -> ([ sum((x.-0.5).^2 .+ 1)
            sum(cos.(x))
            sum(sin.(x))],
          0)
@profile for i=1
    pareto, population = Fresa.solve(f, x, a, b;
                                     archiveelite = false,
                                     npop=20, ngen=300,
                                     #output=100,
                                     tolerance=0.01)
    println("*** profile data")
    Profile.print(format=:flat, sortedby=:count)

    println("*** Pareto front:")
    println(population[pareto])

    using PyPlot
    z = [population[pareto[i]].z for i in 1:length(pareto)];
    PyPlot.plot3D([z[i][1] for i=1:length(z)],
                  [z[i][2] for i=1:length(z)],
                  [z[i][3] for i=1:length(z)],
                  "ro")
    PyPlot.savefig("x3.pdf")
end

3.5 rosenbrock

using Fresa
nx = 2
x0 = 0.5*ones(nx)
a = zeros(nx)
b = 10*ones(nx)
rosenbrock(x) = ([(1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2], 0)
# f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0)
best, pop = Fresa.solve(rosenbrock, x0, a, b)
println("Population at end: $pop")
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")

3.6 simple objective function

This test uses a simple quadratic objective function, defined within. All points are feasible within the domain defined by the lower and upper bounds. All Fresa settings are the defaults.

using Fresa
nx = 2
x0 = 0.5*ones(nx)
a = zeros(nx)
b = 10*ones(nx)
f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0.0)
@time best, pop = Fresa.solve(f, x0, a, b)
println("Population at end: $pop")
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")

4 Recent change history

2019-07-11 updated single objective test cases for new calling sequence
2019-07-11 minor documentation changes to reflect change in return values
2019-07-11 added a major change note for new return values
2019-07-11 added example of post-search analysis of evolution for MINLP example
2019-07-11 including ancestor field passes simple tests
2019-07-11 return the best Point object for single objective case
2019-07-11 added ancestor entry to Point type
2019-07-10 add elite members to the record of solutions generated
2019-06-11 added generation number to vector plot data
2019-06-11 added feature for plotting solutions and their neighbours