[Arjen Markus] (2 october 2018) There exists a wide collection of algorithms for performing "global" optimisation of some (often) complicated function. (In my line of work they need not be mathematical functions of one or more variables that can be written down like: f(x,y) = x**2 + y**2 or something similar, but they are often the outcome of extensive computer programs where dozens of input parameters determine the exact output and the optimisation problem is to find the values for these input parameters that cause the outcome to be as close as possible to measurement values.) A few of those: genetic algorithms, simulated annealing, shuffled complex evolution and particle swarm optimisation. Many of these algorithms have variants and almost all depend on a fair number of options. The fun part is that there is no generally accepted method to determine a priori the best values of the options for a particular problem. Another thing they have in common is that there is a probabilistic component that is used to get to the ''global'' optimum. At least that is the theory. **Implementation** Okay, it is perhaps a bit of hobbyism, but I implemented two versions of the (a?) particle swarm optimisation method. No comprehensive documentation yet, except for the comments in the source code. ====== # pso.tcl -- # Straightforward implementation of particle swarm optimisation # # Options: # - The "best" is taken to be the globally best particle (global) # - The "best" is taken to be the best of the neighbourhood of # each particle # # TODO: # Implement a dictionary with intermediate results # # pso -- # Public interface to the PSO implementations # # Arguments: # func Function to be minimized # xmin Minimum values for the coordinates # xmax Maximum values for the coordinates # args Key-value pairs: # -swarmsize number Number of particles to consider (default: 50) # -vweight value Weight for the current "velocity" (0-1, default: 0.5) # -pweight value Weight for the individual particle's best position (0-1, default: 0.3) # -gweight value Weight for the "best" overall position as per particle (0-1, default: 0.3) # -type local/global Type of optimisation # -neighbours number Size of the neighbourhood (default: 5, used if "local") # -iterations number Maximum number of iterations # -tolerance value Absolute minimal improvement for minimum value # # Result: # Vector with the coordinates of the "best" position and the value # proc pso {func xmin xmax args} { set options [dict create -swarmsize 50 -vweight 0.5 -pweight 0.3 -gweight 0.3 \ -type global -neighbours 5 -iterations 50 -tolerance 0.0] foreach {key value} \$args { if { [dict key \$options \$key] != "" } { dict set options \$key \$value } else { return -code error "Unknown option: \$key" } } if { [llength \$xmin] != [llength \$xmax] } { return -code error "The lists of minimum and maximum values for the coordinates should have the same length" } set type [dict get \$options -type] if { \$type == "global" } { return [Pso_global \$func \$xmin \$xmax \$options] } if { \$type == "local" } { return [Pso_local \$func \$xmin \$xmax \$options] } return -code error "Unknown type: \$type" } # Pso_global -- # Use the "global" PSO algorithm # # Arguments: # func Function to be minimized # xmin Minimum values for the coordinates # xmax Maximum values for the coordinates # options Dictionary of options # proc Pso_global {func xmin xmax options} { # # Set up the initial positions # set swarmsize [dict get \$options -swarmsize] set maxiters [dict get \$options -iterations] set particle_bests {} set positions {} set velocities {} set global_best -1 set global_best_value {} for {set i 0} {\$i < \$swarmsize} {incr i} { set coords [Pso_position \$xmin \$xmax] set fvalue [\$func \$coords] lappend positions \$coords lappend velocities [lrepeat [llength \$coords] 0.0] lappend particle_bests [list \$fvalue \$coords] if { \$global_best == -1 || \$global_best_value > \$fvalue } { set global_best \$i set global_best_value \$fvalue } } set vweight [dict get \$options -vweight] set pweight [dict get \$options -pweight] set gweight [dict get \$options -gweight] set tolerance [dict get \$options -tolerance] for {set iteration 0} {\$iteration < \$maxiters} {incr iteration} { set new_positions {} set new_velocities {} set new_particle_bests {} # # Determine the new positions for all particles # for {set i 0} {\$i < \$swarmsize} {incr i} { set old_velocity [lindex \$velocities \$i] set old_position [lindex \$positions \$i] set old_bestvalue [lindex \$particle_bests \$i 0] set old_bestpos [lindex \$particle_bests \$i 1] set velocity [Pso_update_vel \$vweight \$pweight \$gweight \$old_velocity \$old_position \$old_bestpos [lindex \$positions \$global_best]] set position [Pso_new_position \$old_position \$velocity] set fvalue [\$func \$position] lappend new_positions \$position lappend new_velocities \$velocity if { \$fvalue < \$old_bestvalue } { lappend new_particle_bests [list \$fvalue \$position] } else { lappend new_particle_bests [list \$old_bestvalue \$old_position] } } set positions \$new_positions set velocities \$new_velocities set particle_bests \$new_particle_bests # # Determine the globally best position # for {set i 0} {\$i < \$swarmsize} {incr i} { set fvalue [lindex \$particle_bests \$i 0] if { \$fvalue < \$global_best_value } { set global_best_value \$fvalue set global_best \$i set global_best_pos [lindex \$particle_bests \$i 1] } } # # Have we reached the tolerance yet? # if { \$iteration > 0 } { if { \$prev_best_value - \$global_best_value < \$tolerance && \$prev_best_value != \$global_best_value } { break } } set prev_best_value \$global_best_value } return [list \$global_best_pos \$global_best_value] } # Pso_local -- # Use the "local" PSO algorithm, i.e. look only at the "neighbours" # # Arguments: # func Function to be minimized # xmin Minimum values for the coordinates # xmax Maximum values for the coordinates # options Dictionary of options # proc Pso_local {func xmin xmax options} { # # Set up the initial positions # set swarmsize [dict get \$options -swarmsize] set maxiters [dict get \$options -iterations] set neighbours [dict get \$options -neighbours] set particle_bests {} set positions {} set velocities {} set global_best -1 set global_best_value {} for {set i 0} {\$i < \$swarmsize+\$neighbours} {incr i} { lappend neighbours_idx [expr {\$i % \$swarmsize}] } for {set i 0} {\$i < \$swarmsize} {incr i} { set coords [Pso_position \$xmin \$xmax] set fvalue [\$func \$coords] lappend positions \$coords lappend velocities [lrepeat [llength \$coords] 0.0] lappend particle_bests [list \$fvalue \$coords] lappend local_best \$i lappend local_best_value \$fvalue lappend local_best_pos \$coords lappend prev_best_value \$fvalue } # # Initial estimates: # Examine the neighbouring particles and determine which holds the best result # for {set i 0} {\$i < \$swarmsize} {incr i} { set current_best [lindex \$local_best \$i] set current_best_value [lindex \$particle_bests \$current_best 0] set current_best_pos [lindex \$particle_bests \$current_best 1] for {set n 0} {\$n < \$neighbours} {incr n} { set nth [lindex \$neighbours_idx [expr {\$i+\$n}]] set fvalue [lindex \$particle_bests \$nth 0] if { \$current_best_value > \$fvalue } { set current_best \$nth set current_best_value \$fvalue set current_best_pos [lindex \$particle_bests \$nth 1] } } lset local_best \$i \$current_best lset local_best_value \$i \$current_best_value lset local_best_pos \$i \$current_best_pos } # # Actual loop # set vweight [dict get \$options -vweight] set pweight [dict get \$options -pweight] set gweight [dict get \$options -gweight] set tolerance [dict get \$options -tolerance] set stop 0 for {set iteration 0} {\$iteration < \$maxiters && \$stop == 0} {incr iteration} { set new_positions {} set new_velocities {} set new_particle_bests {} # # Determine the new positions for all particles # for {set i 0} {\$i < \$swarmsize} {incr i} { set old_velocity [lindex \$velocities \$i] set old_position [lindex \$positions \$i] set old_bestvalue [lindex \$particle_bests \$i 0] set old_bestpos [lindex \$particle_bests \$i 1] set idx [lindex \$local_best \$i] set velocity [Pso_update_vel \$vweight \$pweight \$gweight \$old_velocity \$old_position \$old_bestpos [lindex \$positions \$idx]] set position [Pso_new_position \$old_position \$velocity] set fvalue [\$func \$position] lappend new_positions \$position lappend new_velocities \$velocity if { \$fvalue < \$old_bestvalue } { lappend new_particle_bests [list \$fvalue \$position] } else { lappend new_particle_bests [list \$old_bestvalue \$old_position] } } set positions \$new_positions set velocities \$new_velocities set particle_bests \$new_particle_bests # # Examine the neighbouring particles and determine which holds the best result # for {set i 0} {\$i < \$swarmsize} {incr i} { set current_best [lindex \$local_best \$i] set current_best_value [lindex \$local_best_value \$i] set current_best_pos [lindex \$local_best_pos \$i] for {set n 0} {\$n < \$neighbours} {incr n} { set nth [lindex \$neighbours_idx [expr {\$i+\$n}]] set fvalue [lindex \$particle_bests \$nth 0] if { \$current_best_value > \$fvalue } { set current_best \$nth set current_best_value \$fvalue set current_best_pos [lindex \$particle_bests \$nth 1] } } lset local_best \$i \$current_best lset local_best_value \$i \$current_best_value lset local_best_pos \$i \$current_best_pos # # Have we reached the tolerance yet? # Note: local citerium - one group reaching a minimum? Then stop # if { \$iteration > 0 } { if { [lindex \$prev_best_value \$i] - \$current_best_value < \$tolerance && [lindex \$prev_best_value \$i] != \$current_best_value } { set stop 1 break } } lset prev_best_value \$i \$current_best_value } } # # Now determine the overall best position # set global_best_value [lindex \$local_best_value 0] set global_best_pos [lindex \$local_best_pos 0] for {set i 1} {\$i <\$swarmsize} {incr i} { set particle_best_value [lindex \$local_best_value \$i] if { \$global_best_value > \$particle_best_value } { set global_best_value \$particle_best_value set global_best_pos [lindex \$local_best_pos \$i] } } return [list \$global_best_pos \$global_best_value] } # Pso_position -- # Determine the initial position # # Arguments: # xmin Minimum values for the coordinates # xmax Maximum values for the coordinates # # Result: # Vector of coordinates # proc Pso_position {xmin xmax} { set new_position {} foreach min \$xmin max \$xmax { lappend new_position [expr {\$min + (\$max - \$min) * rand()}] } return \$new_position } # Pso_new_position -- # Update the position # # Arguments: # position Position vector # velocity Velocity vector # # Result: # Vector of new coordinates # proc Pso_new_position {position velocity} { set new_position {} foreach p \$position v \$velocity { lappend new_position [expr {\$p + \$v}] } return \$new_position } # Pso_update_vel -- # Update the velocity # # Arguments: # vweight Weight for the old vector # pweight Weight for the particle's best position # gweight Weight for the globally (locally) best position # old_velocity Old velocity vector # old_position Old position vector # old_particle_best Old best vector for the particle # old_global_best Old globally (locally) best vector # # Result: # Vector of new coordinates # proc Pso_update_vel {vweight pweight gweight old_velocity old_position old_particle_best old_global_best} { set new_velocity {} set pw [expr {\$pweight * rand()}] set gw [expr {\$gweight * rand()}] foreach v \$old_velocity p \$old_position b \$old_particle_best g \$old_global_best { lappend new_velocity [expr {\$vweight * \$v + \$pw * (\$b - \$p) + \$gw * (\$g - \$p)}] } return \$new_velocity } ====== **Some simple tests** To put the implementation to the test and get at least some experience with the method, I wrote a few simple test cases. For these cases it does its work nicely, but as you can see, the function is very simple - it has only one minimum, so it is not a real challenge. The challenge was in the implementation instead :) ====== # test -- # Simple test: f(x,y) = (x-1)**2 + (y-1)**2, x and y in [0,4] # proc f {coords} { lassign \$coords x y return [expr {(\$x-1.0)**2 + (\$y-1.0)**2}] } set coords_value [pso f {0.0 0.0} {4.0 4.0} -iterations 50] puts "\$coords_value -- [f [lindex \$coords_value 0]]" set coords_value [pso f {0.0 0.0} {4.0 4.0} -iterations 50 -tolerance 1.0e-4] puts "\$coords_value -- [f [lindex \$coords_value 0]]" set coords_value [pso f {0.0 0.0} {4.0 4.0} -iterations 50 -type local] puts "\$coords_value -- [f [lindex \$coords_value 0]]" set coords_value [pso f {0.0 0.0} {4.0 4.0} -iterations 50 -type local -tolerance 0.1e-4] puts "\$coords_value -- [f [lindex \$coords_value 0]]" ====== <>Category Mathematics