# The Fresa plant propagation algorithm for optimization

## Table of Contents

- 1. Overview and introduction
- 2. Fresa – The code and documentation
- 2.1. start of module and dependencies
- 2.2. types
- 2.3. create a point
- 2.4. fitness
- 2.5. neighbour – generate random point
- 2.6. pareto – set of non-dominated points
- 2.7. printHistoryTrace - show history of a given solution
- 2.8. prune - control population diversity
- 2.9. randompopulation – for testing other methods
- 2.10. select – choose a member of the population
- 2.11. solve – use the PPA to solve the optimisation problem
- 2.12. thinout – make Pareto set smaller
- 2.13. utility functions
- 2.14. module end

- 3. Tests
- 4. Recent change history

Author |
Professor Eric S Fraga (e.fraga@ucl.ac.uk) |

Last revision |
2022-08-07 13:16:13 |

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.

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

### 1.1. Licence and citing

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 relevant sources. These include the software itself (https://doi.org/10.5281/zenodo.5045812):

@misc{fresa, author={Eric S. Fraga}, title={ericsfraga/Fresa.jl: First public release (r2021.06.30)}, journal={Zenodo}, note={{doi:10.5281/zenodo.5045812}}, year={2021}, doi={ }, url={ } }

You may also wish to cite the original paper which led to the development of both *Fresa* and Strawberry, the first implementation of the plant propagation algorithm in MATLAB:

@inproceedings{salhi-fraga-2011a, author={Abdellah Salhi and Eric S. Fraga}, title={Nature-inspired optimisation approaches and the new plant propagation algorithm}, booktitle={Proceedings of ICeMATH 2011, The International Conference on Numerical Analysis and Optimization}, year={2011}, pages={K2:1--8}, address={Yogyakarta}, month={6} }

### 1.2. 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

*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.3. Quick example

Once Fresa is available in your Julia installation, using it requires you to define three aspects to define a problem:

- the objective function;
- the search domain; and,
- an initial population.

I will illustrate how to use Fresa for a simple nonlinear problem, adapted from chapter 13 of a textbook from MIT:^{1}

subject to the following constraints

and with and .

We start by telling Julia that we wish to use the Fresa package:

```
using Fresa
```

The next step is to define the objective function. This function has two responsibilities: it must calculate the value of the objective and it must indicate whether the given point in the search space is feasible or not. The function returns a *tuple* consisting of *z*, the objective function value, and *g*, the indication of feasibility. *g* should be ≤ 0 if the point is feasible and greater than 0 otherwise. For constraints as in the example given above, the most straightforward approach can be to rewrite the constraints in the form :

With this transformation, the objective function can be written:

function objective(x) # calculate the objective function value z = 5*x[1]^2 + 4*x[2]^2 - 60*x[1] - 80*x[2] # evaluate the constraints so that feasible points result in a # non-positive value, i.e. 0 or less, but infeasible points give a # positive value. We choose the maximum of both constraints as # the value to return as an indication of feasibility g = maximum( [ 6*x[1] + 5*x[2] - 60 10*x[1] + 12*x[2] - 150 ] ) # return the objective function value along with indication of # feasibility (z, g) end

The second requirement is the definition of the search domain. For flexibility, for instance to allow the use of problem specific data structures, Fresa expects the search domain to be a function of the search points. The domain therefore is defined by providing two functions, one which returns the lower bounds for the given point in the search space and the other returning the upper bounds. For problems involving a search in the real number domain, such as the example above, it is usually straightforward to define a box search domain.

In the example above, the second optimization variable is unbounded above. However, looking at the constraints and taking into account the domain for the first optimization variable, we can determine that must hold for feasible points. The search domain, for Fresa, is therefore defined as follows:

domain = Fresa.Domain( x -> [ 0.0, 0.0 ], # lower bounds x -> [ 8.0, 12.5 ] ) # upper bounds

This code says that for any point in the search space, `x`

, the lower bounds are given by the vector `[0.0, 0.0]`

and the upper bounds by `[8.0, 12.5]`

Finally, an initial population must be provided to Fresa. This population can be of any size so long as there is at least one member. Fresa usually works well even if only one initial point in the search domain is provided. We consider starting at the midpoint of the search domain defined above and create a point in the search domain:

initialpopulation = [ Fresa.createpoint( [4.0, 6.25 ], objective) ]

The `Fresa.createpoint`

function expects two arguments: the actual point in the search domain and the actual objective function.

Given the above code, Fresa can now be used to solve the problem:

best, population = Fresa.solve( objective, initialpopulation, domain )

The arguments given here for the `solve`

function are those that are required. There are also a number of optional arguments, as described in the code section below

If we execute all the above lines of code in Julia (see a Julia file with this code), the output will be similar to this:

# -*- mode: org; -*- : Fresa PPA last change [2022-08-06 16:04] objective (generic function with 1 method) Fresa.Domain(var"#1#3"(), var"#2#4"()) 1-element Vector{Fresa.Point}: f([4.0, 6.25])=[-503.75] g=-4.75 ** solve objective [2022-08-06 16:04] | variable | value | |- | ngen | 100 | | npop | 10 | | nrmax | 5 | | ns | 100 | | elite | true | | archive | false | | fitness | hadamard | | steepness | 1.0 | |- : function evaluations performed sequentially. *** initial population |- | z1 | g | x | |- | -503.75 | -4.75 | [4.0, 6.25] | |- *** evolution | gen | npop | nf | pruned | t (s) | z1 | g | |- | 1 | 1 | 1 | 0 | 1.62 | -503.75 | -4.75 | | 2 | 2 | 2 | 0 | 2.24 | -503.75 | -4.75 | | 3 | 3 | 4 | 0 | 2.24 | -508.98024613972393 | -3.963481299254326 | | 4 | 10 | 14 | 1 | 2.24 | -516.708162060443 | -2.25158356765359 | | 5 | 26 | 40 | 2 | 2.24 | -520.6840859905025 | -1.1924432568181942 | | 6 | 29 | 68 | 2 | 2.24 | -523.7646499331361 | -0.21175198131829376 | | 7 | 30 | 103 | 8 | 2.30 | -523.7646499331361 | -0.21175198131829376 | | 8 | 24 | 129 | 11 | 2.30 | -524.1229955815547 | -0.35088947745457233 | | 9 | 30 | 160 | 13 | 2.30 | -524.6890696989392 | -0.1186704914001453 | | 10 | 28 | 188 | 14 | 2.30 | -525.265180854414 | -0.7525770581036753 | | 20 | 18 | 458 | 67 | 2.30 | -529.2960864238084 | -0.07945091191650988 | | 30 | 23 | 764 | 125 | 2.31 | -529.6657551053829 | -0.017628945114793737 | | 40 | 22 | 1028 | 191 | 2.31 | -529.6657551053829 | -0.017628945114793737 | | 50 | 28 | 1289 | 249 | 2.31 | -529.6657551053829 | -0.017628945114793737 | | 60 | 15 | 1575 | 364 | 2.32 | -529.6657551053829 | -0.017628945114793737 | | 70 | 33 | 1894 | 416 | 2.32 | -529.6657551053829 | -0.017628945114793737 | | 80 | 21 | 2182 | 495 | 2.33 | -529.6657551053829 | -0.017628945114793737 | | 90 | 22 | 2451 | 545 | 2.34 | -529.6657551053829 | -0.017628945114793737 | | 100 | 19 | 2741 | 597 | 2.34 | -529.6657551053829 | -0.017628945114793737 | *** Fresa run finished : nf=2767 npruned=601 (f([3.7020260052569407, 7.554043004668713])=[-529.6657551053829] g=-0.017628945114793737, |- | z1 | g | x | |- | -529.6657551053829 | -0.017628945114793737 | [3.7020260052569407, 7.554043004668713] | | -366.54435704507523 | -27.21876533995116 | [1.7692502532766028, 4.433146628077845] | | -414.2529968479937 | -20.951491873449086 | [1.3966036701257092, 6.133777221159331] | | -528.8463299506853 | -0.22900672389539523 | [3.6795486014901346, 7.53874033343276] | | -530.384414390094 | 0.16910579472867937 | [3.7209424919011393, 7.568690168664369] | | -443.1260111837601 | -16.228285953791016 | [2.86399171858194, 5.317552746943469] | | -445.6210305308239 | -15.601715891265954 | [3.016294412254267, 5.260103527041689] | | -414.6744628574651 | -20.448352990168345 | [2.5321364109887003, 4.871765708779892] | | -429.18276623025895 | -18.390573964068352 | [2.672394747864535, 5.1150115097488875] | | -527.1298300049249 | -0.5625086104374191 | [3.825709902044894, 7.296646395458643] | | -524.7657659486813 | -1.1402676069200197 | [3.7873107530029255, 7.227173575012485] | | -524.5668961935723 | -1.1915319178131085 | [3.780900324334527, 7.224613227235946] | | -525.8468778251139 | -0.8906651408825468 | [3.7913260301461618, 7.272275735648097] | | -527.7529757793898 | -0.4269495828330747 | [3.8153534947185843, 7.336185889771085] | | -528.3084551642806 | -0.347630580465065 | [3.734607564840653, 7.448944806098202] | | -528.6635168151147 | -0.22202810180591115 | [3.801579393238273, 7.39369910775289] | | -392.9526971218319 | -22.862738531064636 | [2.5817576982271797, 4.329343055914457] | | -457.4148735934202 | -13.81074842840058 | [3.112459710647433, 5.5028986615429645] | | -424.5002108062361 | -12.353125520526376 | [4.648811676286212, 3.9508008843512705] | | -423.1306398623829 | -15.348331900055584 | [3.973866475752854, 4.161693849085458] | | -427.5637392336674 | -11.136583646876801 | [4.8401741154869224, 3.964474332040334] | | -359.63180208155165 | -15.950095856635826 | [5.141747031659222, 2.6398843906817673] | | -414.0489650018614 | -9.615662219485117 | [5.402713812562485, 3.593610981027995] | |- )

The output includes details on the settings of all tunable parameters for the method (all of which can be adjusted, as noted above), the best solution in the population as it evolves, and the best in the final population along with that full population at the end. Note that the output is formatted to be best viewed using `org`

mode^{2} in the Emacs^{3} text editor but the output should be readable as it is all just text.

A more complex 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.^{4} It illustrates the generic nature of Fresa, allowing its application to problems with domain specific data structures.

### 1.4. Version information

Major version log:

- June 2021
**v7.2**, first public release via github and Zenodo: https://doi.org/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**,- 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; - 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.

- the domain for the search, which has to be bounded, is now defined by a
- 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-08-07 13:14]" using Dates using Distributed using Printf function __init__() if myid() == 1 println("# -*- mode: org; -*-") println("#+startup: show3levels") 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`

^{5}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, , 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:

- 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). - 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:

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:

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 and where and are the two values of the tuple and 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*. Some^{6} use ≼ (`\preceq`

); others^{7} 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 for .

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-08-07 13:16 fix internal link 2022-08-07 13:10 rewrote introductory example text 2022-08-06 16:20 Merge branch 'development' 2022-08-06 16:12 upload example.jl and LaTeX images for maths in example 2022-08-06 16:06 Added illustrative example early on for solving NLP problems 2022-08-05 15:15 use show2levels startup directive 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-06 16:26 updated citation section 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

A full history is at the github site.

## Footnotes:

^{1}

Exercise 9, page 452: //web.mit.edu/15.053/www/AMP-Chapter-13.pdf

^{2}

^{4}

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