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-06-01 17:37:18 |
Latest news | API change for solve method. See version information below. |
1. Overview and introduction
Fresa is an implementation of a plant propagation algorithm in Julia programming language. Our original version was called Strawberry, written in Octave (and usable in MATLAB). Please see the documentation of Strawberry for a longer description of the history of these codes, the conditions on which you may use these codes and which articles to cite. The quite interesting story of the strawberry plant can be found elsewhere.
The move to using Julia for the implementation of this algorithm is motivated by the desire to use open source and free (as in beer) software only while also ensuring that the tools are computationally efficient. Octave has not kept up with developments in just in time compilation and, although has some features not available in MATLAB, it is not a modern language. Julia provides many of the capabilities of Octave while also being computationally efficient and providing many more modern features such as built-in parallelism and multi-dispatch, both of which are used by Fresa for efficiency and to make the search procedure generic with respect to data types used to represent points in the search space.
All code (either version) is copyright © 2022, Eric S Fraga. See the LICENSE file for licensing information. Please note that the code is made available with no warranty or guarantee of any kind. Use at own risk. Feedback, including bug reports, is most welcome.
Please let the author know if you download the code. Further, if either Fresa or Strawberry codes are used and this use leads to publications of any type, please do cite the publications listed in the first paragraph of the documentation for Strawberry. The software itself may be cited using 10.5281/zenodo.5045812.
A tutorial example, the design of the operating temperature profile for a simple reactor, modelled as a dynamic optimization problem, is available. This example was the basis for a paper.1
1.1. Installation
This document has been written using org mode in the Emacs text editor. Org allows for literate programming and uses tangling to generate the actual source files for the code. The code, tangled from this file, can be found at github as a Julia package. To use this code, you should add the package using the appropriate Julia tools. Fresa has been registered as a Julia package in the General registry. Therefore, all you should need to do to use it is
add Fresa
Fresa makes use of the following Julia packages: Dates
, Distributed
, and Printf
. The test problems below (and found in the test folder on github) may require one or more of the following extra packages: BenchmarkTools
, PyPlot
, Profile
.
1.2. Version information
Major version log:
- June 2021
- v7.2, first public release via github and Zenodo: 10.5281/zenodo.5045812.
- May 2021
- v7.1, implemented multithreading in the evaluation of the population for each generation. This introduces a new option for the
solve
method:multithreading
which can be set to eithertrue
orfalse
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 orauto
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 ofPoint
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 genericVector
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 singleFresa.Point
value is returned for the best point found, along with the full final population as an array ofFresa.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 ofFloat64
values, the user must create aFresa.neighbour
function definition for the specific data structure type. - September 2017
- v2, moved to an object representation for points in the search space and allowed for parallel evaluation of the objective function when multiple processors are available.
- November 2016
- v1, first Julia plant propagation algorithm implementation.
A list summary of recent change history is given below.
2. Fresa – The code and documentation
The Fresa method is a population based evolutionary algorithm which mimics the propagation used by plants. Throughout the module, the population object is an array of Point
objects. Each point is a point in a search space, the objective function values for this point and a feasibility indication with g≤0
feasible and g>0
infeasible. See the documentation for the solve
method below for more details on the data structures used and expected.
2.1. start of module and dependencies
Fresa depends on a number of packages that should be available in any Julia installation. These are packages for displaying output and using parallel computing capabilities when available on the actual hardware.
# All code copyright © Eric S Fraga. # Date of last change in version variable below. module Fresa version = "[2022-06-01 17:36]" using Dates using Distributed using Printf function __init__() if myid() == 1 println("# -*- mode: org; eval: (org-content 3); -*-") println(": Fresa PPA last change $version") end end
2.2. types
2.2.1. Point
Fresa uses one type, Point
, which is a point in the search space. It includes these entries:
x
- of indeterminate type to allow for a wide range of applications (e.g. integer versus real values),
z
- the value of the objective function, as a vector, where the entries in the vector can be of any type that can be compared and sorted by
sortperm
2 or, in the case of multiple criteria, where it can be determined whether one point dominates another, g
- the constraint violation (feasible with ≤0 and infeasible otherwise) always of type
Float64
(for now), and ancestor
- another point in the search space, along with some extra information, that led to the creation of this point.
An instance of a point is defined by the variable in the search space, the objective function used to evaluate the point, the ancestor of this point (see below), and optional parameters to pass to that function.
As Fresa is an evolutionary procedure, every point in the search space considered will be the descendent of a previously considered point. The sole exception is the initial starting point given by the procedure which invokes Fresa. The link between points is through a backward chain defined by the ancestor entry. This is not used by Fresa itself directly but provides extra meta-information that could be useful for post-optimization analysis, e.g. to see how effective the balance between exploration and exploitation may be for the given search parameter values.
An issue in Julia (as of 2021, at least) is that you cannot define two data structures that mutually refer to each other. Therefore, the type of the ancestor
entry in the Point
data structure has to be defined later (see Ancestor
definition below). This is discussed in the issue for Julia on github.
""" Point (`x`) in the search space along with objective function values (`z[]`) and feasbility indication (`g`). The type of `x` is problem specific. `z[]` and `g` hold `Float64` values. `g` should be of length 1. """ struct Point x :: Any # decision point z :: Vector # objective function values g :: Float64 # constraint violation ancestor # the parent of this point end
Customise how a Point is displayed. We display the objective function value(s) first and then the representation of the point. This allows for a population to have different representations without causing problems with any data analysis on the columns representing the objective function values.
import Base Base.show(io::IO, p::Fresa.Point) = print(io, "f(", p.x, ")=", p.z, " g=", p.g) # and also an array of points function Base.show(io::IO, p::Array{Point,1}) np = length(p) if np > 0 nz = length(p[1].z) println(io, "|-") for i=1:nz print(io,"| z$(i) ") end println(io, "| g | x |") println(io,"|-") for i=1:length(p) for j=1:nz print(io,"| ", p[i].z[j], " ") end print(io, "| ", p[i].g, " ") print(io, "| ", p[i].x, " |\n") end println(io,"|-") else print(io,"empty") end end
and also indicate that a Point
is atomic in a sense:
import Base.size Base.size(p :: Point) = ()
2.2.2. Ancestor
The creation of any point in the search is based on one of the existing points in the population. This existing point is known as the ancestor of the new point. The Ancestor
data structure is used to connect points to their ancestors and collect information about when and how the new point was created.
struct Ancestor point :: Point # the actual ancestor point fitness :: Float64 # the fitness of the ancestor generation :: Int32 # the generation when this point was created end
Once the Ancestor
data structure has been defined, we can now use an access constructor to define the type for the ancestor
field in the Point
object:
ancestor(p :: Point) = p.ancestor :: Union{Ancestor,Nothing}
2.2.3. Domain
Fresa assumes a bounded domain for the search. Each design variable will have a lower and upper bound. To provide for domain specific design variable data structures, the Domain
structure is used. In this structure, the lower
and upper
variables are functions which will be evaluated with a point in the search space and are expected to return appropriate data that the neighbour
function (see below) will be able to use to ensure the domain bounds are respected in the creation of new search points.
struct Domain lower # function which returns lower bound on search variable(s) upper # function which returns upper bound on search variable(s) end
An example of a the use of this Domain
structure is:
d = Domain(x -> zeros(length(x)), x -> ones(length(x)))
which will define the domain as a unit square, , 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 theremotecall
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. Some3 use ≼ (\preceq
); others4 use ≻ (\succ
). I have decide to use the latter as it gives the impression of dominating.
function dominates(a, b) all(a .<= b) && any(a .< b) end ≻(a,b) = dominates(a,b)
This operator will be extended by other packages that wish to make comparisons between non-scalar values of each objective function. The easiest way may often be to ensure that ≤ and < operators are defined for the individual entries in the vector of objective function values.
2.6.2. find Pareto set
The following code splits a population into those points that are non-dominated (i.e. would be considered an approximation to a Pareto frontier) and those that are dominated. The function returns indices into the population passed to it.
function paretoindices(z) n = length(z) dominance = [reduce(&, [!(z[i] ≻ z[j]) for i ∈ 1:n]) for j ∈ 1:n] paretoindices = filter(j -> dominance[j], 1:n) dominatedindices = filter(j -> !dominance[j], 1:n) (paretoindices, dominatedindices) end
Given a population of Point
objects, this function identifies those that are non-dominated (see above). If the population includes both feasible and infeasible points, only those that are feasible are considered.
# indices of non-dominated and dominated points from the population of # Point objects function pareto(pop :: Vector{Point}) l = length(pop) indexfeasible = (1:l)[map(p->p.g,pop) .<= 0] indexinfeasible = (1:l)[map(p->p.g,pop) .> 0] if length(indexfeasible) > 0 subset = view(pop,indexfeasible) indices = indexfeasible else #println(": Fresa.pareto warning: no feasible solutions. Pareto set meaningless?") subset = pop indices = 1:l end z = map(p->p.z, subset) # use function below to return indices of non-dominated and # dominated from objective function values alone in the subset of # feasible solutions (p, d) = paretoindices(z) (indices[p], indices[d]) end
2.7. printHistoryTrace - show history of a given solution
Each point encountered in the search, other than points in the initial population, is the result of propagating another point. When a new point is created, a link back to its parent point is created. This allows us to explore the history of all points in the search. This function prints out the historical trace of a given point, using an org
table for formatting.
function printHistoryTrace(p :: Point) a = p.ancestor while ! (a isa Nothing) println("| $(a.generation) | $(a.fitness) |") a = a.point.ancestor end end
2.8. prune - control population diversity
Due to the stochastic nature of the method and also the likely duplication of points when elitism is used, there is a need to prune the population. We wish to remove members that have objective function values (or alternatively decision variable values) that are too close to each other. The main difficulty is the definition of too close. We use a tolerance based on the range of values present in the population for the objective function pruning and the tolerance directly for decision variable pruning.
Previously, pruning was done on objective function values. In the case where that is not possible (e.g. cannot find difference of values), we consider the decision variables as well. The latter assume that we have a -(subtraction) operator for the decision variable type; if not, we do no pruning at all.
function prune(pop :: AbstractArray, tolerance) npruned = 0 z = map(p->p.z, pop) # @show z[1] # println("typeof(z)=$(typeof(z))") l = length(z) # println("typeof(z[1])=$(typeof(z[1]))") n = length(z[1]) # @show n zmin = zeros(n) zmax = zeros(n) try for i=1:n row = [z[j][i] for j=1:l] zmin[i] = minimum(row) zmax[i] = maximum(row) if zmax[i] - zmin[i] < 100*eps() zmax[i] = zmin[i]+100*eps() end end pruned = [pop[1]] for i=2:l similar = false for j=1:length(pruned) if all(abs.(z[i]-pruned[j].z) .< tolerance*(zmax-zmin)) similar = true; break; end end if !similar push!(pruned,pop[i]) else npruned += 1 end end (pruned, npruned) catch e # println("prune method error: $e") if e isa MethodError # probably (possibly) due to objective function type not # being a number. In this case, we try again but looking # at the decision variable values instead. x = map(p->p.x, pop) # println("typeof(z)=$(typeof(z))") l = length(x) # start building up the population that remains after # pruning. The first entry will always be there as any # similar solutions will not be included by the search # that follows. pruned = [pop[1]] try for i=2:l similar = false # now check this solution against all those already in # the list we are collating for j=1:length(pruned) if all(Float64.(abs.(x[i]-pruned[j].x)) .< tolerance) similar = true; break; end end if !similar push!(pruned,pop[i]) else npruned += 1 end end (pruned, npruned) catch e # println("prune method second error: $e") if e isa MethodError # this is now probably/possibly due to not being # to find the difference between two decision # points. In that case, return the whole # original population (pop, 0) else @error "Unexpected error in prune method for points" e end # if method error end # try pruning on points else @error "Unexpected error in prune method for objective function values" e end # if is a method error end # try pruning on objective function values end # function
2.9. randompopulation – for testing other methods
Create a random population of size n
evaluated using f
. A single point, x
, in the search domain must be given as the domain definition is function based and the lower and upper bounds are potentially a function of the location in the space. The randompoint
method below is suitable for domains defined by float valued vectors.
function randompopulation(n, f, parameters, p0, domain :: Domain) p = Point[] # population object for j in 1:n # l = domain.lower(p0.x) # @show l # u = domain.upper(p0.x) # @show u # x = randompoint(l,u) # push!(p, createpoint(x, f, parameters)) push!(p, createpoint(randompoint(domain.lower(p0.x), domain.upper(p0.x)), f, parameters)) end p end
By default, the following method generates a random point within the search domain. This does not attempt to find a feasible point, simply one within the box defined by lower, a
, and upper, b
, bounds.
function randompoint(a :: Float64, b :: Float64) x = a + rand()*(b-a) end function randompoint(a, b) x = a + rand(length(a)).*(b-a) end
2.10. select – choose a member of the population
Given a fitness, f
, choose two solutions randomly and select the one with the better fitness. This is known as a tournament selection procedure of size 2. Other options are possible but not currently implemented.
function select(f) l = length(f) ind1 = rand(1:l) if ind1 == 0 ind1 = 1 end ind2 = rand(1:l) # println("Comparing $ind1 to $ind2") if f[ind1] > f[ind2] return ind1 else return ind2 end end
2.11. solve – use the PPA to solve the optimisation problem
The function expects the objective function, f
, an initial population, p0
, with at least one point, and the Domain
for the search. It returns the optimum, the objective function value(s) at this point, the constraint at that point and the whole population at the end. The actual return values and data structures depends on the number of criteria:
- 1
- returns best point as a
Fresa.Point
object (which includes the decision variable values, the objective function value, and the constraint value) and also the full population; - >1
- returns the set of non-dominated points (as an array including objective function values and constraint value) and the full population.
The objective function, f
, should return two results: z
, the objective function value(s) which must be of type Float64
, single or array, and g
, the constraint violation. If g≤0
, the point is feasible; any value g>0
means an infeasible point. The value of g
for infeasible points will be used to rank the fitness of the infeasible solution, with lower values being fitter.
The calling sequence for f
is a point in the search space plus, optionally, the parameters
defined in the call to solve
.
p0
is a population with initial points in the search space; these points can be of any type. domain
is a valid Domain
object with appropriate functions for determining the lower and upper bounds of the search space in terms of the optimization variables. These should be consistent with the representations use for the individual points in the search space.
If the decision vector is not an array of Float64
, a type specific Fresa.neighbour
function will need to be defined. The calling sequence for Fresa.neighbour
is (x,a,b,fitness)
where x
, a
, and b
, should all be of the desired type and the function itself must also return an object of that type. The fitness
will always be a Float64
. See the mixed integer nonlinear example below for an example.
The fitnesstype
is used for ranking members of a population for multi-objective problems. The default is to use a Hadamard product of the rank each solution has for each objective individually. One alternative, specifying fitnesstype=:borda
uses a sum of the rank, i.e. a Borda count. The former tends to emphasise points near the extrema of the individual criteria while the latter is possibly better distributed but possibly at providing less emphasis on the Pareto points themselves. There is also the option fitnesstype=:nondominated
which bases the fitness on levels of dominance, as used by the NSGA-II
genetic algorithm.
The size of the population, npop
, may be a single integer value or a Tuple
of two integer values. The latter, which is only for multi-objective optimization problems, gives a range of possible values for the population size. This size will be chosen dynamically within this range depending on the size of the non-dominated set at the start of each generation. Specifically, the population will be set to 2 times that size. This allows for sufficient diversity in the population while minimizing computation time. It has been seen that Fresa is largely insensitive to the population size: there is an interesting video by Marleen de Jonge & Daan van den Berg discussing the robustness of the plant propagation algorithm with respect to the parameters for the algorithm, using a slightly different version of the algorithm which does not use tournament selection but instead selects the top npop
members of the population for propagation.
The output of the progress during the search is controlled by the output
optional argument. This should be an integer value that indicates how often a summary of the current population is generated and sent to standard output. It will be the initial value used. The value will go up in powers of 10 as the generations proceed to ensure that there is sufficient granularity without overwhelming the output file. The default is 1 to output every generation until the 10th, then 10 until the 100th, and so on. A value of 0 will eliminate all output from the solve method.
""" Solve an optimisation problem, defined as the minimization of the values returned by the objective function, `f`. `f` returns not only the objective function values, an array of `Float64` values, but also a measure of feasibility (≤0) or infeasibility (>0). The problem is solved using the Fresa algorithm. `p0` is the initial population which has to have at least one member, a `Point`, and `a` and `b` are *bounds* on the search space. The return values for the solution of a single criterion problem are the best point and the full population at the end of the search. For a multi-objective problem, the returned values are the set of indices for the points within the full population (the second returned value) approximating the *Pareto* front. The population will consist of an array of `Fresa.Point` objects, each of which will have the point in the search space, the objective function value and the feasibility measure. """ function solve(f, p0, domain; # required arguments parameters = nothing, # allow parameters for objective function archiveelite = false, # save thinned out elite members elite = true, # elitism by default fitnesstype = :hadamard, # how to rank solutions in multi-objective case multithreading = false, # use multiple threads for objective function evaluation ngen = 100, # number of generations npop = 10, # population size: fixed (single value) or dynamic (tuple) nrmax = 5, # number of runners maximum ns = 100, # number of stable solutions for stopping output = 1, # how often to output information plotvectors = false, # generate output file for search plot populationoutput = false, # output population every generation? steepness = 1.0, # show steep is the adjustment shape for fitness tolerance = 0.001, # tolerance for similarity detection usemultiproc = false) # parallel processing by Fresa itself? output > 0 && println("** solve $f $(orgtimestamp(now()))") tstart = time() nf = 1 # number of function evaluations npruned = 0 # number solutions pruned from population nz = length(p0[1].z) # number of criteria pop = copy(p0); # create/initialise the population object if archiveelite archive = Point[] end if output > 0 println("#+name: $(f)settings") println("| variable | value |") println("|-") println("| ngen | $ngen |") println("| npop | $npop |") println("| nrmax | $nrmax |") println("| ns | $ns |") println("| elite | $elite |") println("| archive | $archiveelite |") println("| fitness | $fitnesstype |") println("| steepness | $steepness |") println("|-") # output != 0 && println(": solving with ngen=$ngen npop=$npop nrmax=$nrmax ns=$ns") # output != 0 && println(": elite=$elite archive elite=$archiveelite fitness type=$fitnesstype") end if plotvectors plotvectorio = open("fresa-vectors-$(orgtimestamp(now())).data", create=true, write=true) output > 0 && println(": output of vectors for subsequent plotting") end # if npop was given as a tuple, we are to have a dynamic # population size. This only makes sense for multi-objective # optimization problems so a warning will be given otherwise. npopmin = npop npopmax = npop if npop isa Tuple if nz > 1 npopmin = npop[1] npopmax = npop[2] if npopmin > npopmax error("Dynamic population sizing requires min <= max; you specified $npop") end npop = npopmin # start with minimum possible else println("*Warning*: you have specified a tuple for population size: $npop") println("This only makes sense for multi-objective optimization problems.") println("npop will be set to $(npop[1]).") npop = npop[1] # be optimistic and use minimum given end end # we use multithreading if asked for *and* if we have more than # one thread available multithreading = multithreading && Threads.nthreads() > 1 # we use parallel computing if we have more than one processor parallel = usemultiproc && nprocs() > 1 # parallel = false if output > 0 println(": function evaluations performed ", parallel ? "in parallel with $(nprocs()) processors." : (multithreading ? "in parallel with $(Threads.nthreads()) threads." : "sequentially.")) println("*** initial population") println("#+name: $(f)initial") println(pop) end if output > 0 println("*** evolution") println("#+name: $(f)evolution") println("#+plot: ind:1 deps:(6) with:\"points pt 7\" set:\"logscale x\"") @printf("| %9s | %9s | %9s | %9s | %9s |", "gen", "npop", (elite && nz > 1) ? "pareto" : "nf", "pruned", "t (s)") for i in 1:nz @printf(" z%-8d |", i) end @printf(" %9s |", "g") @printf("\n|-\n") end # now evolve the population for a predetermined number of generations for gen in 1:ngen # evaluate fitness which is adjusted depending on value of # steepness, a value that may depend on the generation fit = fitness(pop, fitnesstype, steepness, gen, ngen) if gen == 1 @debug "Initial fitness" f=fit end # sort index = sortperm(fit) if populationoutput println("\nGeneration $gen full population is:") println(pop) println("Fitness vector: $fit") end # and remember best which really only makes sense in single # criterion problems but is best in multi-objective case in # the ranking measure used by Fresa best = pop[index[end]] # if elitism is used if elite if nz > 1 # elite set is whole pareto set unless it is too # big. Recall that the pareto function returns the set # of indices into the population wholepareto = pareto(pop)[1] # if using dynamic population sizing, adjust the population npop = 2 * length(wholepareto) if npop < npopmin npop = npopmin end if npop > npopmax npop = npopmax end # now check that the pareto is not too big. if it is, thin it out if length(wholepareto) > ceil(npop/2) newpop, removed = thinout(pop, fit, wholepareto, ceil(Int,npop/2)) if archiveelite archive = prune(append!(archive, removed), tolerance)[1] archive = archive[pareto(archive)[1]] end else newpop = pop[wholepareto] end else # elite set is single element only newpop = [best] end # if plotting vectors for the search, include elitism if plotvectors for p in newpop write(plotvectorio, "$(gen-1) $(p.x)\n$gen $(p.x)\n\n") end end else newpop = Point[] end if output >= 0 print(stderr, ": $gen np=$(length(newpop))/$npop", archiveelite ? " na=$(length(archive))" : "", " with most fit z=$(best.z) \r") # if output has been requested, check to see if output is # required now and then also check to see if the frequency # needs to be reduced. if output > 0 if gen%output == 0 @printf("| %9d | %9d | %9d | %9d | %9.2f |", gen, length(fit), (elite && nz > 1) ? length(newpop) : nf, npruned, time()-tstart) for i = 1:length(best.z) print(" $(best.z[i]) |") end print(" $(best.g) |") println() end if 10^(floor(log10(gen))) > output output = 10^(Int(floor(log10(gen)))) end end end # if we are using any form of multiprocessing, either threads # or multiple cores, create an array to store all new points # which we evaluate later in parallel. Ideally, also keep # track of the points from which new points are derived to # provide the backward link through the evolution but this is # currently disabled as the creation of the Ancestor object # requires more information than I am currently storing away. if multithreading || parallel x = Any[] # typeof(newpop[1].x)[] # points = Point[] end # now loop through population, applying selection and then # generating neighbours l = length(pop) for i in 1:min(l,npop) s = select(fit) # println(": selection $i is $s") # println(": size of pop is $(size(pop))") selected = pop[s] if !elite # if no elitism, we ensure selected members remain in population push!(newpop, selected) if plotvectors write(plotvectorio, "$(gen-1) $(selected.x)\n$gen $(selected.x)\n\n") end end # number of runners to generate, function of fitness nr = ceil(fit[s]*nrmax*rand()) if nr < 1 nr = 1 end # println(": generating $nr runners") for r in 1:nr # create a neighbour, also function of fitness, # passing through the lower and upper bounds # appropriate for the particular solution point newx = neighbour(pop[s].x, domain.lower(pop[s].x), domain.upper(pop[s].x), fit[s]) # for parallel evaluation, we store the neighbours and # evaluate them later; otherwise, we evaluate # immediately and save the resulting point if multithreading || parallel push!(x, newx) # push!(points, pop[s]) else push!(newpop, createpoint(newx, f, parameters, Ancestor(pop[s],fit[s],gen))) if plotvectors write(plotvectorio, "$(gen-1) $(pop[s].x)\n$gen $newx\n\n") end nf += 1 end end # remove selected member from the original population so # it is not selected again splice!(fit, s) splice!(pop, s) end # if we are making use of parallel computing, we evaluate all # points generated in previous loop. Parallel processing is # done either via multithreading or with multiple # processors. The former is easier as it's based on shared # memory. if multithreading # using threads and shared memory results = Array{Point}(undef,length(x)) Threads.@threads for i ∈ 1:length(x) results[i] = createpoint(x[i],f,parameters) end append!(newpop, results) elseif parallel # using multiple processors with remote calls # will be used to collect results from worker processors results = Array{Future,1}(undef, nprocs()) i = 0; while i < length(x) # issue remote evaluation call for j=1:nprocs() if i+j <= length(x) # TODO: the information about the ancestor is # not available; this needs to be stored above results[j] = @spawn createpoint(x[i+j],f,parameters) nf += 1 end end # now wait for results for j=1:nprocs() if i+j <= length(x) push!(newpop, fetch(results[j])) end end i += nprocs() end end # and finally, if we have elitism, remove any duplicate points # in the new population and make it the current population for # the next generation; otherwise, simply copy over if elite (pop, nn) = prune(newpop, tolerance) npruned += nn else pop = newpop end end output > 0 && println("*** Fresa run finished\n: nf=$nf npruned=$npruned", archiveelite ? " archived=$(length(archive))" : "") if plotvectors close(plotvectorio) end if nz == 1 fit = fitness(pop, fitnesstype, steepness, ngen, ngen) index = sortperm(fit) best = pop[index[end]] return best, pop else return pareto(archiveelite ? append!(pop,archive) : pop)[1], pop end end
2.12. thinout – make Pareto set smaller
If we use elitism, for multi-objective problems, we use the Pareto set as the elite set. However, this set may grow to be large, causing performance challenges as well as making the search less effective at exploration, essentially getting stuck in the local area defined by this elite set. Therefore, we need to sometimes thin out the Pareto set for its use as an elite set.
The arguments are the whole population, the fitness of the members, the indices in this population for the Pareto set and the number of elements to keep. We keep the most fit ones.
function thinout(pop, fit, pareto, n::Int) indices = sortperm(fit[pareto]) return pop[pareto[indices[end-n+1:end]]], pop[pareto[indices[1:end-n]]] end
2.13. utility functions
Some functions that are not necessary for Fresa but provide some useful features, especially output related.
2.13.1. org time stamp
function orgtimestamp(dt::DateTime) return @sprintf("[%d-%02d-%02d %02d:%02d]", Dates.year(dt), Dates.month(dt), Dates.day(dt), Dates.hour(dt), Dates.minute(dt)) end
2.13.2. rank – dimension of a variable
Sometimes, we need to determine whether a variable (e.g. the objective function value returned by the evaluation of the model) is a scalar or a vector.
rank(x :: Any) = length(size(x))
2.14. module end
end
3. Tests
The following are simple tests for either the Fresa optimiser or just individual functions in the module. You can cut and paste these codes into your own editor and run them.
3.1. simple objective function
This test uses a simple quadratic objective function, defined within. All points are feasible within the domain defined by the lower and upper bounds. All Fresa settings are the defaults.
# load in the Fresa optimization package using Fresa # specify the dimension of the search space nx = 2 # create an initial point in the search space x0 = 0.5*ones(nx) # specify the domain for the search, x ∈ [0,10]ⁿ domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x))) # the actual objective function f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0) # create the initial population consisting of this single point p0 = [Fresa.createpoint(x0,f)] # now invoke Fresa to solve the problem @time best, pop = Fresa.solve(f, p0, domain) # output the results println("Population at end:\n$pop") println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")
One of the features that Fresa provides is a trace of how each solution has been created. That is, each solution has a link back to the ancestor solution that led to its creation, along with information about when this happened (the generation) and how fit the ancestor solution was. There is a function defined in Fresa for outputting a history trace. The output is in form of an org mode
table but is simple text that can be imported into a spreadsheet program, for instance.
println("\nHistory trace, by generation number, of fitness value of solution selected for propagation which results in a new best solution:") println("#+plot: ind:1 deps:(2) with:\"linespoints pt 7 ps 0.25\" set:nokey set:\"yrange [0:1]\" set:\"xrange [0:*]\" set:\"xlabel 'Generation'\" set:\"ylabel 'fitness'\"") Fresa.printHistoryTrace(best)
3.2. rosenbrock
using Fresa nx = 2 x0 = 0.5*ones(nx) # specify the domain for the search, x ∈ [0,10]ⁿ domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x))) rosenbrock(x) = ([(1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2], 0) # f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0) # create the initial population consisting of this single point p0 = [Fresa.createpoint(x0,rosenbrock)] # now invoke Fresa to solve the problem best, pop = Fresa.solve(rosenbrock, p0, domain; ngen=1000, tolerance=1e-8) println("Population at end: $pop") println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")
3.3. rosenbrock higher dimensions
The generalised Rosenbrock function is
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-06-01 17:37 improved HTML export for code listings, using css highlighting 2021-12-03 13:11 updated random search procedures for case study on fitness types 2021-12-03 11:59 labelled the domain type definition src block 2021-12-03 11:58 generalised the randompopulation method to use domain types 2021-12-03 11:55 fixed version update code although not actually used currently 2021-12-02 16:34 added JuMP to do item: interface to JuMP for use with Fresa 2021-08-05 17:01 use float instead of Float64 for conversion 2021-07-29 17:25 replaced all isa() and typeof() with x isa t constructs 2021-07-29 15:25 added description text in output for simple test case 2021-07-29 15:12 commented out some debug statements in prune method 2021-07-27 16:37 action for use of isa infix operator 2021-07-27 15:56 added Fresa specific actions from our todo list 2021-07-09 14:28 removed old instructions for using the package 2021-07-09 14:27 fixed package name in instructions for adding the package 2021-07-03 13:59 moved to semantic versioning, as required by Julia Registry 2021-06-30 14:22 added release news 2021-06-30 14:21 added link to LICENSE file and the DOI for the software release 2021-06-30 13:07 added optional argument to control of output of population 2021-05-28 11:10 changed name of web page to match bibliography entry 2021-05-26 12:25 added news item on multithreading capability
A full history is at the github site.
Footnotes:
E S Fraga (2019), An example of multi-objective optimization for dynamic processes, Chemical Engineering Transactions 74:601-606, 10.3303/CET1974101.