# 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