The Fresa plant propagation algorithm for optimization

Table of Contents

Author Professor Eric S Fraga (e.fraga@ucl.ac.uk)
Last revision 2022-06-01 17:37:18
Latest news API change for solve method. See version information below.

1. Overview and introduction

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 quite interesting story of the strawberry plant can be found elsewhere.

The move to using Julia for the implementation of this algorithm 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 © 2022, Eric S Fraga. See the LICENSE file for licensing information. Please note that the code is made available with no warranty or guarantee of any kind. Use at own risk. Feedback, including bug reports, is most welcome.

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. The software itself may be cited using 10.5281/zenodo.5045812.

A tutorial example, the design of the operating temperature profile for a simple reactor, modelled as a dynamic optimization problem, is available. This example was the basis for a paper.1

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 at github as a Julia package. To use this code, you should add the package using the appropriate Julia tools. Fresa has been registered as a Julia package in the General registry. Therefore, all you should need to do to use it is 

add Fresa

in the Julia package manager.

Fresa makes use of the following Julia packages: Dates, Distributed, and Printf. The test problems below (and found in the test folder on github) may require one or more of the following extra packages: BenchmarkTools, PyPlot, Profile.

1.2. Version information

Major version log:

June 2021
v7.2, first public release via github and Zenodo: 10.5281/zenodo.5045812.
May 2021
v7.1, implemented multithreading in the evaluation of the population for each generation. This introduces a new option for the solve method: multithreading which can be set to either true or false with the latter being the default. Julia must be invoked with the --threads argument (or -t for short) with the number of threads to use or auto for automatic determination of the threads possible. Multithreading is useful when the evaluation of the objective function is computationally expensive. Otherwise, the overhead of multithreading is usually not worth it although it is not detrimental.
April 2021
v7,
  1. the domain for the search, which has to be bounded, is now defined by a Domain data structure which allows for different representations of solutions in the search space in a given population;
  2. allow setting the steepness of the fitness adjustment function. This is an outcome of the presentation by Wouter Vrielink at the PPA mini-zoomposium I organised in March 2021 to discuss the impact of PPA parameters on the effectiveness of the search procedure.
March 2021
v6, one of the required arguments to the solve function has been changed. Specifically, the initial guess must now be a population of Point objects and not a single decision variable. See examples below for how to create this initial population easily.
April 2020
moved all code to github. This should make it easier for others to use the code.
September 2019
v5, The objective function values, in the Point type, are now a generic Vector instead of an array of floating point numbers. This opens up Fresa to be used for objective functions which are not necessarily simple scalar values. The use case has been illustrated through a case study in stochastic optimization, specifically design under uncertainty. Details available from the author.
July 2019
v4, 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
v3, 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
v2, 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
v1, 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. 

# All code copyright © Eric S Fraga. 
# Date of last change in version variable below.
module Fresa
version = "[2022-06-01 17:36]"
using Dates
using Distributed
using Printf
function __init__()
    if myid() == 1
        println("# -*- mode: org; eval: (org-content 3); -*-")
        println(": Fresa PPA last change $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 these entries:

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, as a vector, where the entries in the vector can be of any type that can be compared and sorted by sortperm 2 or, in the case of multiple criteria, where it can be determined whether one point dominates another,
g
the constraint violation (feasible with ≤0 and infeasible otherwise) always of type Float64 (for now), and
ancestor
another point in the search space, along with some extra information, that led to the creation of this point.

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.

An issue in Julia (as of 2021, at least) is that you cannot define two data structures that mutually refer to each other. Therefore, the type of the ancestor entry in the Point data structure has to be defined later (see Ancestor definition below). This is discussed in the issue for Julia on github. 

"""

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 :: Vector                 # objective function values
    g :: Float64                # constraint violation
    ancestor                    # the parent of this point
end

Customise how a Point is displayed. We display the objective function value(s) first and then the representation of the point. This allows for a population to have different representations without causing problems with any data analysis on the columns representing the objective function values. 

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, "|-")
        for i=1:nz
            print(io,"| z$(i) ")
        end
        println(io, "| g | x |")
        println(io,"|-")
        for i=1:length(p)
            for j=1:nz
                print(io,"| ", p[i].z[j], " ")
            end
            print(io, "| ", p[i].g, " ")
            print(io, "| ", p[i].x, " |\n")
        end
        println(io,"|-")
    else
        print(io,"empty")
    end
end

and also indicate that a Point is atomic in a sense: 

import Base.size
Base.size(p :: Point) = ()

2.2.2. Ancestor

The creation of any point in the search is based on one of the existing points in the population. This existing point is known as the ancestor of the new point. The Ancestor data structure is used to connect points to their ancestors and collect information about when and how the new point was created. 

struct Ancestor
    point :: Point        # the actual ancestor point
    fitness :: Float64    # the fitness of the ancestor
    generation :: Int32   # the generation when this point was created
end

Once the Ancestor data structure has been defined, we can now use an access constructor to define the type for the ancestor field in the Point object: 

ancestor(p :: Point) = p.ancestor :: Union{Ancestor,Nothing}

2.2.3. Domain

Fresa assumes a bounded domain for the search. Each design variable will have a lower and upper bound. To provide for domain specific design variable data structures, the Domain structure is used. In this structure, the lower and upper variables are functions which will be evaluated with a point in the search space and are expected to return appropriate data that the neighbour function (see below) will be able to use to ensure the domain bounds are respected in the creation of new search points. 

struct Domain
    lower                       # function which returns lower bound on search variable(s)
    upper                       # function which returns upper bound on search variable(s)
end

An example of a the use of this Domain structure is: 

d = Domain(x -> zeros(length(x)), x -> ones(length(x)))

which will define the domain as a unit square, \(x \in [0,1]^n\), as defined by the size of the x argument.

2.3. create a point

A trivial function that simply creates a new Point object. This exists for two reasons:

  1. It 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).
  2. The Point type is parametric. This makes defining a generic constructor difficult (at least, I was unable to find a working solution).

The optional parameters and ancestor arguments are passed through to their respective destinations: the objective function for the parameters and the point creation for the ancestor linking. 

function createpoint(x,f,parameters = nothing,ancestor = nothing)
    z = 0
    g = 0
    if ! ( parameters isa Nothing )
        (z, g) = f(x, parameters)
    else
        (z, g) = f(x)
    end
    if g isa Int
        g = float(g)
    end
    p = Nothing
    if rank(z) == 1
        p = Point(x, z, g, ancestor)
    elseif rank(z) == 0
        p = Point(x, [z], g, ancestor)
    else
        error("Fresa can only handle scalar and vector criteria, not $(typeof(z)).")
    end
    return p
end

(deprecated) and we provide two versions with simple calling sequences: 

function createpoint(x,f)
    return createpoint(x,f,nothing,nothing)
end
function createpoint(x,f,parameters)
    return createpoint(x,f,parameters,nothing)
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) with the same properties: best solutions have fitness values closer to 1 than less fit solutions.

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

function fitness(pop, fitnesstype, steepness, generation, ngen)
    l = length(pop)
    indexfeasible = (1:l)[map(p->p.g,pop) .<= 0]
    indexinfeasible = (1:l)[map(p->p.g,pop) .> 0]
    @debug "Feasible/infeasible breakdown" indexfeasible indexinfeasible maxlog=3
    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, steepness, generation, ngen)
        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,
                                             steepness,
                                             generation,
                                             ngen
                                             ) / 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 three alternative measures of
fitness ranking possible.  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 based on
each criterion ranking.  Finally, a measure based on dominance,
similar to that used by the popular NSGA-II genetic algorithm, is
available.

"""
function vectorfitness(v, fitnesstype, steepness, generation, ngen)
    # 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 
            fitness = [v[i][1] for i=1:l]
        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)
            elseif fitnesstype == :borda
                for i=1:m
                    rank[i,sortperm([v[j][i] for j=1:l])] = 1:l
                end
                # borda sum of ranks
                fitness = map(x->sum(x), rank[:,i] for i=1:l)
            elseif fitnesstype == :nondominated
                # similar to that used by NSGA-II (Deb 2000)
                fitness = zeros(l)
                maxl = assigndominancefitness!(fitness,v,1)
                # println("Resulting fitness: $fitness")
            else
                throw(ArgumentError("Type of fitness evaluation must be either :borda, :nondominated, or :hadamard, not $(repr(fitnesstype))."))
            end
        end
        # normalise (1=best, 0=worst) while avoiding
        # extreme 0,1 values using the hyperbolic tangent
        fit = adjustfitness(fitness, steepness, generation, ngen)
        # println(":  scaled fitness: $fit")
        @debug "Fitness calculations" v[1][1] v[2][1] v[l][1] fitness[1] fitness[2] fitness[l] fit[1] fit[2] fit[l] maxlog=3
    end
    fit
end

The fitness should be a value ∈ (0,1), i.e. not including the bounds themselves as those values cause some silly behaviour in the definition of individual neighbouring solutions (i.e. the runners) and the number of runners. Therefore, we adjust the fitness values to ensure that the bounds are not included.

See below for a discussion about the second function argument, steepness, and how the value s is calculated if steepness is a tuple and not a single value. 

function adjustfitness(fitness, steepness, generation, ngen)
    if (maximum(fitness)-minimum(fitness)) > eps()
        s = steepness
        if steepness isa Tuple
            a = (2*steepness[1]-2*steepness[2])/3
            b = - (3*steepness[1] - 3*steepness[2])/ngen^2
            d = steepness[1]
            s = a*generation^3 + b*generation^2 + c*generation + d
            @debug "Steepness " s "at generation" g
        end  
        fit = 0.5*(tanh.(4*s*(maximum(fitness) .- fitness)
                         / (maximum(fitness)-minimum(fitness))
                         .- 2*s) .+ 1)
    else
        # only one solution (or all solutions the same) in population
        fit = 0.5*ones(length(fitness))
    end
    fit
end

This function takes, as an argument, the steepness of the transition from poor fitness to good fitness. Some plots are useful for comparison. This first plot shows the default fitness adjustment function which gives some emphasis to the extreme values but also ensures that the fitness values are quite some distance from the boundary of the fitness domain:

defaultsteepness.png

Making the fitness adjustment steeper, e.g. with a value of steepness of 2 instead of the default value of 1, the function has a more pronounced emphasis towards the boundaries and allows values closer to those boundaries:

sharpersteepness.png

The steepness may be specified as a tuple in which case it represents the initial value for the steepness and the final value. The evolution of the steepness is based on a cubic with 0 slope at the start and at the end. The following maxima code is the solution of the that cubic given the need to pass through the points \((0,s_1)\) and \((n_g,s_2)\) where \(s_1\) and \(s_2\) are the two values of the tuple and \(n_g\) is the number of generations: 

c(g) := a*g^3 + b*g^2 + c*g + d;
define(d(g), diff(c(g),g));
equations: [c(0) = s1,
d(0) = 0,
c(n) = s2,
d(n) = 0];
solution: solve(equations, [a, b, c, d]);
for i: 1 thru length(solution[1]) do print(solution[1][i])$

    2 s1 - 2 s2
a = ----------- 
         3
        n
      3 s1 - 3 s2
b = - ----------- 
           2
          n
c = 0 
d = s1

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) = paretoindices(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).

2.6.1. dominates: determine dominance

To cater for generic comparisons between points in the objective function space (e.g. distributions instead of single values for each objective function), we introduce an operator used to determine dominance. The community differs on the symbol to use for dominates. Some3 use ≼ (\preceq); others4 use ≻ (\succ). I have decide to use the latter as it gives the impression of dominating. 

function dominates(a, b)
    all(a .<= b) && any(a .< b)
end
≻(a,b) = dominates(a,b)

This operator will be extended by other packages that wish to make comparisons between non-scalar values of each objective function. The easiest way may often be to ensure that ≤ and < operators are defined for the individual entries in the vector of objective function values.

2.6.2. find Pareto set

The following code splits a population into those points that are non-dominated (i.e. would be considered an approximation to a Pareto frontier) and those that are dominated. The function returns indices into the population passed to it. 

function paretoindices(z)
    n = length(z)
    dominance = [reduce(&, [!(z[i] ≻ z[j]) for i ∈ 1:n]) for j ∈ 1:n]
    paretoindices = filter(j -> dominance[j], 1:n)
    dominatedindices = filter(j -> !dominance[j], 1:n)
    (paretoindices, dominatedindices)
end

Given a population of Point objects, this function identifies those that are non-dominated (see above). If the population includes both feasible and infeasible points, only those that are feasible are considered. 

# indices of non-dominated and dominated points from the population of
# Point objects
function pareto(pop :: Vector{Point})
    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) = paretoindices(z)
    (indices[p], indices[d])
end

2.7. printHistoryTrace - show history of a given solution

Each point encountered in the search, other than points in the initial population, is the result of propagating another point. When a new point is created, a link back to its parent point is created. This allows us to explore the history of all points in the search. This function prints out the historical trace of a given point, using an org table for formatting. 

function printHistoryTrace(p :: Point)
    a = p.ancestor
    while ! (a isa Nothing)
        println("| $(a.generation) | $(a.fitness) |")
        a = a.point.ancestor
    end
end

2.8. 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 (or alternatively decision variable 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 for the objective function pruning and the tolerance directly for decision variable pruning.

Previously, pruning was done on objective function values. In the case where that is not possible (e.g. cannot find difference of values), we consider the decision variables as well. The latter assume that we have a -(subtraction) operator for the decision variable type; if not, we do no pruning at all. 

function prune(pop :: AbstractArray, tolerance)
    npruned = 0
    z = map(p->p.z, pop)
    # @show z[1]
    # println("typeof(z)=$(typeof(z))")
    l = length(z)
    # println("typeof(z[1])=$(typeof(z[1]))")
    n = length(z[1])
    # @show n
    zmin = zeros(n)
    zmax = zeros(n)
    try 
        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)
    catch e
        # println("prune method error: $e")
        if e isa MethodError
            # probably (possibly) due to objective function type not
            # being a number.  In this case, we try again but looking
            # at the decision variable values instead.
            x = map(p->p.x, pop)
            # println("typeof(z)=$(typeof(z))")
            l = length(x)
            # start building up the population that remains after
            # pruning.  The first entry will always be there as any
            # similar solutions will not be included by the search
            # that follows.
            pruned = [pop[1]]
            try
                for i=2:l
                    similar = false
                    # now check this solution against all those already in
                    # the list we are collating
                    for j=1:length(pruned)
                        if all(Float64.(abs.(x[i]-pruned[j].x)) .< tolerance)
                            similar = true;
                            break;
                        end
                    end
                    if !similar
                        push!(pruned,pop[i])
                    else
                        npruned += 1
                    end
                end
                (pruned, npruned)        
            catch e
                # println("prune method second error: $e")
                if e isa MethodError
                    # this is now probably/possibly due to not being
                    # to find the difference between two decision
                    # points.  In that case, return the whole
                    # original population
                    (pop, 0)
                else
                    @error "Unexpected error in prune method for points" e
                end             # if method error
            end                 # try pruning on points
        else
            @error "Unexpected error in prune method for objective function values" e
        end                     # if is a method error
    end                         # try pruning on objective function values
end                             # function

2.9. randompopulation – for testing other methods

Create a random population of size n evaluated using f. A single point, x, in the search domain must be given as the domain definition is function based and the lower and upper bounds are potentially a function of the location in the space. The randompoint method below is suitable for domains defined by float valued vectors. 

function randompopulation(n, f, parameters, p0, domain :: Domain)
    p = Point[]                 # population object
    for j in 1:n
        # l = domain.lower(p0.x)
        # @show l
        # u = domain.upper(p0.x)
        # @show u
        # x = randompoint(l,u)
        # push!(p, createpoint(x, f, parameters))
        push!(p, createpoint(randompoint(domain.lower(p0.x), domain.upper(p0.x)), f, parameters))
    end
    p
end

By default, the following method generates a random point within the search domain. This does not attempt to find a feasible point, simply one within the box defined by lower, a, and upper, b, bounds. 

function randompoint(a :: Float64, b :: Float64)
    x = a + rand()*(b-a)
end
function randompoint(a, b)
    x = a + rand(length(a)).*(b-a)
end

2.10. 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.11. solve – use the PPA to solve the optimisation problem

The function expects the objective function, f, an initial population, p0, with at least one point, and the Domain for the search. 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.

p0 is a population with initial points in the search space; these points can be of any type. domain is a valid Domain object with appropriate functions for determining the lower and upper bounds of the search space in terms of the optimization variables. These should be consistent with the representations use for the individual points in the search space.

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. One 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. There is also the option fitnesstype=:nondominated which bases the fitness on levels of dominance, as used by the NSGA-II genetic algorithm.

The size of the population, npop, may be a single integer value or a Tuple of two integer values. The latter, which is only for multi-objective optimization problems, gives a range of possible values for the population size. This size will be chosen dynamically within this range depending on the size of the non-dominated set at the start of each generation. Specifically, the population will be set to 2 times that size. This allows for sufficient diversity in the population while minimizing computation time. It has been seen that Fresa is largely insensitive to the population size: there is an interesting video by Marleen de Jonge & Daan van den Berg discussing the robustness of the plant propagation algorithm with respect to the parameters for the algorithm, using a slightly different version of the algorithm which does not use tournament selection but instead selects the top npop members of the population for propagation.

The output of the progress during the search is controlled by the output optional argument. This should be an integer value that indicates how often a summary of the current population is generated and sent to standard output. It will be the initial value used. The value will go up in powers of 10 as the generations proceed to ensure that there is sufficient granularity without overwhelming the output file. The default is 1 to output every generation until the 10th, then 10 until the 100th, and so on. A value of 0 will eliminate all output from the solve method. 

""" 

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.  `p0` is the initial population
which has to have at least one member, a `Point`, 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, p0, domain;        # 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
               multithreading = false, # use multiple threads for objective function evaluation
               ngen = 100,           # number of generations
               npop = 10,            # population size: fixed (single value) or dynamic (tuple)
               nrmax = 5,            # number of runners maximum
               ns = 100,             # number of stable solutions for stopping
               output = 1,           # how often to output information
               plotvectors = false,  # generate output file for search plot
               populationoutput = false, # output population every generation?
               steepness = 1.0,      # show steep is the adjustment shape for fitness
               tolerance = 0.001,    # tolerance for similarity detection
               usemultiproc = false) # parallel processing by Fresa itself?
    output > 0 && println("** solve $f $(orgtimestamp(now()))")
    tstart = time()
    nf = 1                   # number of function evaluations
    npruned = 0              # number solutions pruned from population
    nz = length(p0[1].z)     # number of criteria
    pop = copy(p0);          # create/initialise the population object
    if archiveelite
        archive = Point[]
    end
    if output > 0
        println("#+name: $(f)settings")
        println("| variable | value |")
        println("|-")
        println("| ngen | $ngen |")
        println("| npop | $npop |")
        println("| nrmax | $nrmax |")
        println("| ns | $ns |")
        println("| elite | $elite |")
        println("| archive | $archiveelite |")
        println("| fitness | $fitnesstype |")
        println("| steepness | $steepness |")
        println("|-")
        # output != 0 && println(": solving with ngen=$ngen npop=$npop nrmax=$nrmax ns=$ns")
        # output != 0 && println(": elite=$elite archive elite=$archiveelite fitness type=$fitnesstype")
    end
    if plotvectors
        plotvectorio = open("fresa-vectors-$(orgtimestamp(now())).data", create=true, write=true)
        output > 0 && println(": output of vectors for subsequent plotting")
    end
    # if npop was given as a tuple, we are to have a dynamic
    # population size.  This only makes sense for multi-objective
    # optimization problems so a warning will be given otherwise.
    npopmin = npop
    npopmax = npop
    if npop isa Tuple
        if nz > 1
            npopmin = npop[1]
            npopmax = npop[2]
            if npopmin > npopmax
                error("Dynamic population sizing requires min <= max; you specified $npop")
            end
            npop = npopmin      # start with minimum possible
        else
            println("*Warning*: you have specified a tuple for population size: $npop")
            println("This only makes sense for multi-objective optimization problems.")
            println("npop will be set to $(npop[1]).")
            npop = npop[1]      # be optimistic and use minimum given
        end
    end
    # we use multithreading if asked for *and* if we have more than
    # one thread available
    multithreading = multithreading && Threads.nthreads() > 1 
    # we use parallel computing if we have more than one processor
    parallel = usemultiproc && nprocs() > 1
    # parallel = false
    if output > 0
        println(": function evaluations performed ",
                parallel
                ? "in parallel with $(nprocs()) processors."
                : (multithreading
                   ? "in parallel with $(Threads.nthreads()) threads."
                   : "sequentially."))
        println("*** initial population")
        println("#+name: $(f)initial")
        println(pop)
    end
    if output > 0
        println("*** evolution")
        println("#+name: $(f)evolution")
        println("#+plot: ind:1 deps:(6) with:\"points pt 7\" set:\"logscale x\"")
        @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(" %9s |", "g")
        @printf("\n|-\n")
    end
    # now evolve the population for a predetermined number of generations
    for gen in 1:ngen
        # evaluate fitness which is adjusted depending on value of
        # steepness, a value that may depend on the generation
        fit = fitness(pop, fitnesstype, steepness, gen, ngen)
        if gen == 1
            @debug "Initial fitness" f=fit
        end
        # sort
        index = sortperm(fit)
        if populationoutput
            println("\nGeneration $gen full population is:")
            println(pop)
            println("Fitness vector: $fit")
        end
        # 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 using dynamic population sizing, adjust the population
                npop = 2 * length(wholepareto)
                if npop < npopmin
                    npop = npopmin
                end
                if npop > npopmax
                    npop = npopmax
                end
                # now check that the pareto is not too big.  if it is, thin it out
                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
        if output >= 0
            print(stderr, ": $gen np=$(length(newpop))/$npop",
                  archiveelite ? " na=$(length(archive))" : "",
              " with most fit z=$(best.z)           \r")
            # if output has been requested, check to see if output is
            # required now and then also check to see if the frequency
            # needs to be reduced.
            if output > 0
                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)
                        print(" $(best.z[i]) |")
                    end
                    print(" $(best.g) |")
                    println()
                end
                if 10^(floor(log10(gen))) > output
                    output = 10^(Int(floor(log10(gen))))
                end
            end
        end
        # if we are using any form of multiprocessing, either threads
        # or multiple cores, create an array to store all new points
        # which we evaluate later in parallel.  Ideally, also keep
        # track of the points from which new points are derived to
        # provide the backward link through the evolution but this is
        # currently disabled as the creation of the Ancestor object
        # requires more information than I am currently storing away.
        if multithreading || parallel
            x = Any[] # 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,
                # passing through the lower and upper bounds
                # appropriate for the particular solution point
                newx = neighbour(pop[s].x, domain.lower(pop[s].x), domain.upper(pop[s].x), fit[s])
                # for parallel evaluation, we store the neighbours and
                # evaluate them later; otherwise, we evaluate
                # immediately and save the resulting point
                if multithreading || parallel
                    push!(x, newx)
                    # push!(points, pop[s])
                else
                    push!(newpop, createpoint(newx, f, parameters, Ancestor(pop[s],fit[s],gen)))
                    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.  Parallel processing is
        # done either via multithreading or with multiple
        # processors.  The former is easier as it's based on shared
        # memory.
        if multithreading       # using threads and shared memory
            results = Array{Point}(undef,length(x))
            Threads.@threads for i ∈ 1:length(x)
                results[i] = createpoint(x[i],f,parameters)
            end
            append!(newpop, results)
        elseif parallel        # using multiple processors with remote calls
            # will be used to collect results from worker processors
            results = Array{Future,1}(undef, nprocs())
            i = 0;
            while i < length(x)
                # issue remote evaluation call
                for j=1:nprocs()
                    if i+j <= length(x) 
                        # TODO: the information about the ancestor is
                        # not available; this needs to be stored above
                        results[j] = @spawn createpoint(x[i+j],f,parameters)
                        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, if we have elitism, remove any duplicate points
        # in the new population and make it the current population for
        # the next generation; otherwise, simply copy over
        if elite
            (pop, nn) = prune(newpop, tolerance)
            npruned += nn
        else
            pop = newpop
        end
    end
    output > 0 && 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, steepness, ngen, ngen)
        index = sortperm(fit)
        best = pop[index[end]]
        return best, pop
    else
        return pareto(archiveelite ? append!(pop,archive) : pop)[1], pop
    end
end

2.12. 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.13. utility functions

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

2.13.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.13.2. rank – dimension of a variable

Sometimes, we need to determine whether a variable (e.g. the objective function value returned by the evaluation of the model) is a scalar or a vector. 

rank(x :: Any) = length(size(x))

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

# load in the Fresa optimization package
using Fresa
# specify the dimension of the search space
nx = 2
# create an initial point in the search space
x0 = 0.5*ones(nx)
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x)))
# the actual objective function
f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,f)]
# now invoke Fresa to solve the problem
@time best, pop = Fresa.solve(f, p0, domain)
# output the results
println("Population at end:\n$pop")
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")

One of the features that Fresa provides is a trace of how each solution has been created. That is, each solution has a link back to the ancestor solution that led to its creation, along with information about when this happened (the generation) and how fit the ancestor solution was. There is a function defined in Fresa for outputting a history trace. The output is in form of an org mode table but is simple text that can be imported into a spreadsheet program, for instance. 

println("\nHistory trace, by generation number, of fitness value of solution selected for propagation which results in a new best solution:")
println("#+plot: ind:1 deps:(2) with:\"linespoints pt 7 ps 0.25\" set:nokey set:\"yrange [0:1]\" set:\"xrange [0:*]\" set:\"xlabel 'Generation'\" set:\"ylabel 'fitness'\"")
Fresa.printHistoryTrace(best)

3.2. rosenbrock

using Fresa
nx = 2
x0 = 0.5*ones(nx)
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x)))
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)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,rosenbrock)]
# now invoke Fresa to solve the problem
best, pop = Fresa.solve(rosenbrock, p0, domain; ngen=1000, tolerance=1e-8)
println("Population at end: $pop")
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")

3.3. rosenbrock higher dimensions

The generalised Rosenbrock function is \[ f(x_1,\ldots,x_n) = \sum_{i=1}^{n-1} 100 \left ( x_{i+1}-x_i^2 \right)^2 + (1-x_i)^2 \] for \(n&gt;1\). 

using Fresa
nx = 20
x0 = 0.5*ones(nx)
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x)))
rosenbrock(x) = (sum([100 * (x[i+1]-x[i]^2)^2 + (1-x[i])^2 for i ∈ 1:length(x)-1]), 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,rosenbrock)]
# now invoke Fresa to solve the problem
best, pop = Fresa.solve(rosenbrock, p0, domain; npop=100, ngen=1000, tolerance=1e-8, multithreading=true)
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")

3.4. multi-objective test

using Fresa
nx = 2
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> ones(length(x)))
# initial point in domain
x = rand(nx)
# objective function 
f = x -> ( [sin(x[1]-x[2]); cos(x[1]+x[2])], 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x,f)]
# now invoke Fresa to solve the problem
pareto, population = Fresa.solve(f, p0, domain;
                                 #fitnesstype = :hadamard,
                                 #fitnesstype = :borda,
                                 fitnesstype = :nondominated,
                                 ngen=200,
                                 npop=(20,40),
                                 plotvectors=true,
                                 tolerance=0.01)

println("**** Pareto front:")
println("#+plot: ind:1 deps:(2) with:points")
println(population[pareto])
#using BenchmarkTools
#@benchmark

3.5. multi-objective test with 3 objectives

using Fresa
using Profile
nx = 5
# specify the domain for the search, x ∈ [0,1]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> ones(length(x)))
x = zeros(nx)
f = x -> ([ sum((x.-0.5).^2 .+ 1)
            sum(cos.(x))
            sum(sin.(x))],
          0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x,f)]
# now invoke Fresa to solve the problem
@profile for i=1
    pareto, population = Fresa.solve(f, p0, domain;
                                     archiveelite = false,
                                     npop=20, ngen=300,
                                     #output=100,
                                     tolerance=0.01)

    println("*** Pareto front:")
    println(population[pareto])
end
println("*** profile data")
println(": this may take some time so please wait")
Profile.print(format=:flat, sortedby=:count)

3.6. 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
domain = Fresa.Domain(x -> MI(1.0, 1),
                      x -> 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
# create the initial population consisting of a single MI point
p0 = [Fresa.createpoint(MI(1.0, 1),f)]
# now invoke Fresa to solve the problem
best, pop = Fresa.solve(f, p0, domain; 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: 

println("#+plot: ind:3 deps:(2) with:\"linespoints pt 7\" set:nokey set:\"yrange [0:1]\"")
ancestor = best.ancestor;
while ancestor != Some(nothing) && ! (ancestor isa Nothing)
    global ancestor
    println("| $(ancestor.point.z) | $(ancestor.fitness) | $(ancestor.generation) |")
    ancestor = ancestor.point.ancestor
end

3.7. 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.7.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.7.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.7.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.7.4. solve the multi-objective problem using Fresa

using Fresa
domain = Fresa.Domain(x -> [0.0;0.0;0.0], x -> [5.0;3.0;3.0])
x0 = [4.0;2.0;2.0]
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,fmo)]
# now invoke Fresa to solve the problem
pareto, population = Fresa.solve(fmo, p0, domain;
                                 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.7.5. solve the single objective version

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

4. Recent change history

Summary of the most recent changes to the package: 

2022-06-01 17:37 improved HTML export for code listings, using css highlighting
2021-12-03 13:11 updated random search procedures for case study on fitness types
2021-12-03 11:59 labelled the domain type definition src block
2021-12-03 11:58 generalised the randompopulation method to use domain types
2021-12-03 11:55 fixed version update code although not actually used currently
2021-12-02 16:34 added JuMP to do item: interface to JuMP for use with Fresa
2021-08-05 17:01 use float instead of Float64 for conversion
2021-07-29 17:25 replaced all isa() and typeof() with x isa t constructs
2021-07-29 15:25 added description text in output for simple test case
2021-07-29 15:12 commented out some debug statements in prune method
2021-07-27 16:37 action for use of isa infix operator
2021-07-27 15:56 added Fresa specific actions from our todo list
2021-07-09 14:28 removed old instructions for using the package
2021-07-09 14:27 fixed package name in instructions for adding the package
2021-07-03 13:59 moved to semantic versioning, as required by Julia Registry
2021-06-30 14:22 added release news
2021-06-30 14:21 added link to LICENSE file and the DOI for the software release
2021-06-30 13:07 added optional argument to control of output of population
2021-05-28 11:10 changed name of web page to match bibliography entry
2021-05-26 12:25 added news item on multithreading capability

A full history is at the github site.

Footnotes:

1

E S Fraga (2019), An example of multi-objective optimization for dynamic processes, Chemical Engineering Transactions 74:601-606, 10.3303/CET1974101.

2

https://docs.julialang.org/en/v1/base/sort/index.html

3

https://oklahomaanalytics.com/data-science-techniques/nsga-ii-explained/

4

https://stackoverflow.com/questions/28406878/about-correct-notation-of-pareto-dominance-in-optimization

Author: Professor Eric S Fraga

Created: 2022-06-01 Wed 17:37

Validate