1# annealing.tcl -- 2# Package implementing simulated annealing for minimizing functions 3# of one or more parameters 4# 5# Copyright (c) 2007 by Arjen Markus <arjenmarkus@users.sourceforge.net> 6# 7# See the file "license.terms" for information on usage and redistribution 8# of this file, and for a DISCLAIMER OF ALL WARRANTIES. 9# 10# RCS: @(#) $Id: annealing.tcl,v 1.4 2008/02/22 13:34:07 arjenmarkus Exp $ 11#------------------------------------------------------------------------------ 12 13package require Tcl 8.4 14 15# ::simulation::annealing -- 16# Create the namespace 17# 18namespace eval ::simulation::annealing { 19} 20 21 22# getOption -- 23# Return the value of an option 24# 25# Arguments: 26# option Name of the option (without -) 27# Result: 28# The value or an error message if it does not exist 29# 30proc ::simulation::annealing::getOption {option} { 31 variable ann_option 32 33 if { [info exists ann_option(-$option)] } { 34 return $ann_option(-$option) 35 } else { 36 return -code error "No such option: $option" 37 } 38} 39 40 41# setOption -- 42# Set the value of an option 43# 44# Arguments: 45# option Name of the option (without -) 46# value Value of the option 47# Result: 48# None 49# 50proc ::simulation::annealing::setOption {option value} { 51 variable ann_option 52 53 set ann_option(-$option) $value 54} 55 56 57# hasOption -- 58# Return whether the given option exists or not 59# 60# Arguments: 61# option Name of the option (without -) 62# Result: 63# 1 if it exists, 0 if not 64# 65proc ::simulation::annealing::hasOption {option} { 66 variable ann_option 67 68 if { [info exists ann_option(-$option)] } { 69 return 1 70 } else { 71 return 0 72 } 73} 74 75 76# findMinimum -- 77# Find the (global) minimum of a function using simulated annealing 78# 79# Arguments: 80# args Option-value pairs: 81# -parameters list - triples defining parameters and ranges 82# -function expr - expression defining the function 83# -code body - body of code to define the function 84# (takes precedence over -function) 85# should set the variable "result" 86# -init code - code to be run at start up 87# -final code - code to be run at the end 88# -trials n - number of trials before reducing the temperature 89# -reduce factor - reduce the temperature by this factor 90# (between 0 and 1) 91# -initial-temp t - initial temperature 92# -scale s - scale of the function (order of 93# magnitude of the values) 94# -estimate-scale y/n - estimate the scale (only if -scale not present) 95# -verbose y/n - Turn verbose printing on (1) or off (0) 96# -reportfile file - Channel to write verbose output to 97# Any others can be used via the getOption procedure 98# in the body. 99# 100# Result: 101# Estimated minimum and the parameters involved: 102# function value param1 value param2 value ... 103# 104proc ::simulation::annealing::findMinimum {args} { 105 variable ann_option 106 107 # 108 # Handle the options 109 # 110 set ann_option(-parameters) {} 111 set ann_option(-function) {} 112 set ann_option(-code) {} 113 set ann_option(-init) {} 114 set ann_option(-final) {} 115 set ann_option(-trials) 300 116 set ann_option(-reduce) 0.95 117 set ann_option(-initial-temp) 1.0 118 set ann_option(-scale) {} 119 set ann_option(-estimate-scale) 0 120 set ann_option(-verbose) 0 121 set ann_option(-reportfile) stdout 122 123 foreach {option value} $args { 124 set ann_option($option) $value 125 } 126 127 if { $ann_option(-scale) == {} } { 128 if { ! $ann_option(-estimate-scale) } { 129 set ann_option(-scale) 1.0 130 } 131 } 132 133 if { $ann_option(-code) != {} } { 134 set ann_option(-function) {} 135 } 136 137 if { $ann_option(-code) == {} && $ann_option(-function) == {} } { 138 return -code error "Neither code nor function given! Nothing to optimize" 139 } 140 if { $ann_option(-parameters) == {} } { 141 return -code error "No parameters given! Nothing to optimize" 142 } 143 144 if { $ann_option(-function) != {} } { 145 set ann_option(-code) "set result \[expr {$ann_option(-function)}\]" 146 } 147 148 # 149 # Create the procedure 150 # 151 proc FindMin {} [string map \ 152 [list PARAMETERS $ann_option(-parameters) \ 153 CODE $ann_option(-code) \ 154 INIT $ann_option(-init) \ 155 FINAL $ann_option(-final)] { 156 # 157 # Give all parameters a value 158 # 159 foreach {_param_ _min_ _max_} {PARAMETERS} { 160 set $_param_ $_min_ 161 } 162 163 set _trials_ [getOption trials] 164 set _temperature_ [getOption initial-temp] 165 set _reduce_ [getOption reduce] 166 set _noparams_ [expr {[llength {PARAMETERS}]/3}] 167 set _verbose_ [getOption verbose] 168 set _reportfile_ [getOption reportfile] 169 170 INIT 171 172 # 173 # Estimate the scale 174 # 175 if { [getOption estimate-scale] == 1 } { 176 set _sum_ 0.0 177 for { set _trial_ 0 } { $_trial_ < $_trials_/3 } { incr _trial_ } { 178 set _randp_ [expr {3*int($_noparams_*rand())}] 179 set _param_ [lindex {PARAMETERS} $_randp_] 180 set _min_ [lindex {PARAMETERS} [expr {$_randp_+1}]] 181 set _max_ [lindex {PARAMETERS} [expr {$_randp_+2}]] 182 set _old_param_ [set $_param_] 183 set $_param_ [expr {$_min_ + rand()*($_max_-$_min_)}] 184 185 CODE 186 187 set _sum_ [expr {$_sum_ + abs($result)}] 188 } 189 set _scale_ [expr {3.0*$_sum_/$_trials_}] 190 } else { 191 set _scale_ [getOption scale] 192 } 193 if { $_verbose_ } { 194 puts $_reportfile_ "Scale value: $_scale_" 195 } 196 197 # 198 # Start the outer loop 199 # 200 set _changes_ 1 201 202 # 203 # Get the initial value of the function 204 # 205 CODE 206 set _old_result_ $result 207 208 if { $_verbose_ } { 209 puts $_reportfile_ "Result -- Mean of accepted values -- % accepted" 210 } 211 212 while {1} { 213 set _sum_ $_old_result_ 214 set _accepted_ 1 215 for { set _trial_ 0 } { $_trial_ < $_trials_} { incr _trial_ } { 216 set _randp_ [expr {3*int($_noparams_*rand())}] 217 set _param_ [lindex {PARAMETERS} $_randp_] 218 set _min_ [lindex {PARAMETERS} [expr {$_randp_+1}]] 219 set _max_ [lindex {PARAMETERS} [expr {$_randp_+2}]] 220 set _old_param_ [set $_param_] 221 set $_param_ [expr {$_min_ + rand()*($_max_-$_min_)}] 222 223 CODE 224 225 # 226 # Accept the new solution? 227 # 228 set _rand_ [expr {rand()}] 229 if { log($_rand_) < -($result-$_old_result_)/($_scale_*$_temperature_) } { 230 incr _changes_ 231 set _old_result_ $result 232 set _sum_ [expr {$_sum_ + $result}] 233 incr _accepted_ 234 } else { 235 set $_param_ $_old_param_ 236 } 237 } 238 239 if { $_verbose_ } { 240 puts $_reportfile_ \ 241 [format "%.5g -- %.5g -- %.2f %%" $_old_result_ \ 242 [expr {$_sum_/$_accepted_}] [expr {100.0*double($_changes_)/$_trials_}]] 243 } 244 245 set _temperature_ [expr {$_reduce_ * $_temperature_}] 246 if { $_changes_ == 0 } { 247 break 248 } else { 249 set _changes_ 0 250 } 251 } 252 253 set result [list result $_old_result_] ;# Note: we need the last accepted result! 254 foreach {_param_ _min_ _max_} {PARAMETERS} { 255 lappend result $_param_ [set $_param_] 256 } 257 258 FINAL 259 260 return $result 261 }] 262 263 # 264 # Do the actual computation and return the result 265 # 266 return [FindMin] 267} 268 269 270# findCombinatorialMinimum -- 271# Find the (global) minimum of a combinatorial function using simulated annealing 272# 273# Arguments: 274# args Option-value pairs: 275# -number-params n - number of (binary) parameters 276# -initial-values list - list of parameter values to start with 277# -function expr - expression defining the function 278# -code body - body of code to define the function 279# (takes precedence over -function) 280# should set the variable "result" 281# The values of the solutions 282# are stored as a list in the 283# variable params 284# -init code - code to be run at start up 285# -final code - code to be run at the end 286# -trials n - number of trials before reducing the temperature 287# -reduce factor - reduce the temperature by this factor 288# (between 0 and 1) 289# -initial-temp t - initial temperature 290# -scale s - scale of the function (order of 291# magnitude of the values) 292# -estimate-scale y/n - estimate the scale (only if -scale not present) 293# -verbose y/n - Turn verbose printing on (1) or off (0) 294# -reportfile file - Channel to write verbose output to 295# Any others can be used via the getOption procedure 296# in the body. 297# 298# Result: 299# Estimated minimum and the parameters involved: 300# function value, list of values 301# 302# Note: 303# The parameters have the values 0 or 1 304# 305# The stop criterion is that if the result value does not change in 306# sqrt(trials) then the iteration stops. Experiments with the 307# example below show that the function to be minimised can show 308# a very wide minimum due to the parameters being discrete. 309# sqrt(trials) is just an arbitrary value. 310# 311proc ::simulation::annealing::findCombinatorialMinimum {args} { 312 variable ann_option 313 314 # 315 # Handle the options 316 # 317 set ann_option(-number-params) {} 318 set ann_option(-initial-values) {} 319 set ann_option(-function) {} 320 set ann_option(-code) {} 321 set ann_option(-init) {} 322 set ann_option(-final) {} 323 set ann_option(-trials) 300 324 set ann_option(-reduce) 0.95 325 set ann_option(-initial-temp) 1.0 326 set ann_option(-scale) {} 327 set ann_option(-estimate-scale) 0 328 set ann_option(-verbose) 0 329 set ann_option(-reportfile) stdout 330 331 foreach {option value} $args { 332 set ann_option($option) $value 333 } 334 335 if { $ann_option(-scale) == {} } { 336 if { ! $ann_option(-estimate-scale) } { 337 set ann_option(-scale) 1.0 338 } 339 } 340 341 if { $ann_option(-code) != {} } { 342 set ann_option(-function) {} 343 } 344 345 if { $ann_option(-code) == {} && $ann_option(-function) == {} } { 346 return -code error "Neither code nor function given! Nothing to optimize" 347 } 348 if { $ann_option(-number-params) == {} } { 349 return -code error "Number of parameters not given! Nothing to optimize" 350 } 351 352 if { $ann_option(-initial-values) == {} } { 353 for { set i 0 } { $i < $ann_option(-number-params) } { incr i } { 354 lappend ann_option(-initial-values) 0 355 } 356 } 357 358 if { $ann_option(-function) != {} } { 359 set ann_option(-code) "set result \[expr {$ann_option(-function)}\]" 360 } 361 362 # 363 # Create the procedure 364 # 365 proc FindCombMin {params} [string map \ 366 [list CODE $ann_option(-code) \ 367 INIT $ann_option(-init) \ 368 FINAL $ann_option(-final)] { 369 370 set _trials_ [getOption trials] 371 set _temperature_ [getOption initial-temp] 372 set _reduce_ [getOption reduce] 373 set _noparams_ [llength $params] 374 set _verbose_ [getOption verbose] 375 set _reportfile_ [getOption reportfile] 376 377 INIT 378 379 # 380 # Estimate the scale 381 # 382 if { [getOption estimate-scale] == 1 } { 383 set _sum_ 0.0 384 set _old_params_ $params 385 for { set _trial_ 0 } { $_trial_ < $_trials_/3 } { incr _trial_ } { 386 set _randp_ [expr {int($_noparams_*rand())}] 387 lset params $_randp_ [expr {rand()>0.5? 0 : 1}] 388 389 CODE 390 391 set _sum_ [expr {$_sum_ + abs($result)}] 392 } 393 set _scale_ [expr {3.0*$_sum_/$_trials_}] 394 set params $_old_params_ 395 } else { 396 set _scale_ [getOption scale] 397 } 398 if { $_verbose_ } { 399 puts $_reportfile_ "Scale value: $_scale_" 400 } 401 402 # 403 # Start the outer loop 404 # 405 set _changes_ 1 406 407 # 408 # Get the initial value of the function 409 # 410 CODE 411 set _old_result_ $result 412 set _result_same_ 0 413 set _result_after_loop_ $result 414 415 if { $_verbose_ } { 416 puts $_reportfile_ "Result -- Mean of accepted values -- % accepted" 417 } 418 419 while {1} { 420 set _sum_ $_old_result_ 421 set _accepted_ 1 422 for { set _trial_ 0 } { $_trial_ < $_trials_} { incr _trial_ } { 423 set _old_params_ $params 424 set _randp_ [expr {int($_noparams_*rand())}] 425 lset params $_randp_ [expr {rand()>0.5? 0 : 1}] 426 427 CODE 428 429 # 430 # Accept the new solution? 431 # 432 set _rand_ [expr {rand()}] 433 if { log($_rand_) < -($result-$_old_result_)/($_scale_*$_temperature_) } { 434 incr _changes_ 435 set _old_result_ $result 436 set _sum_ [expr {$_sum_ + $result}] 437 incr _accepted_ 438 } else { 439 set params $_old_params_ 440 } 441 } 442 443 if { $_verbose_ } { 444 puts $_reportfile_ \ 445 [format "%.5g -- %.5g -- %.2f %%" $_old_result_ \ 446 [expr {$_sum_/$_accepted_}] [expr {100.0*double($_changes_)/$_trials_}]] 447 } 448 449 set _temperature_ [expr {$_reduce_ * $_temperature_}] 450 if { $_changes_ == 0 || $_result_same_ > sqrt($_trials_) } { 451 break 452 } else { 453 set _changes_ 0 454 } 455 456 if { $_result_after_loop_ == $_old_result_ } { 457 incr _result_same_ 458 } else { 459 set _result_after_loop_ $_old_result_ 460 } 461 } 462 463 set result [list result $_old_result_] ;# Note: we need the last accepted result! 464 lappend result solution $params 465 466 FINAL 467 468 return $result 469 }] 470 471 # 472 # Do the actual computation and return the result 473 # 474 return [FindCombMin $ann_option(-initial-values)] 475} 476 477# Announce the package 478# 479package provide simulation::annealing 0.2 480 481# main -- 482# Example 483# 484if { 0 } { 485puts [::simulation::annealing::findMinimum \ 486 -trials 300 \ 487 -verbose 1 \ 488 -parameters {x -5.0 5.0 y -5.0 5.0} \ 489 -function {$x*$x+$y*$y+sin(10.0*$x)+4.0*cos(20.0*$y)}] 490 491puts "Constrained:" 492puts [::simulation::annealing::findMinimum \ 493 -trials 3000 \ 494 -reduce 0.98 \ 495 -parameters {x -5.0 5.0 y -5.0 5.0} \ 496 -code { 497 if { hypot($x-5.0,$y-5.0) < 4.0 } { 498 set result [expr {$x*$x+$y*$y+sin(10.0*$x)+4.0*cos(20.0*$y)}] 499 } else { 500 set result 1.0e100 501 } 502 }] 503} 504 505# 506# A simple combinatorial problem: 507# We have 100 items and the function is optimal if the first 10 508# values are 1 and the result is 0. Can we find this solution? 509# 510# What if we have 1000 items? Or 10000 items? 511# 512# WARNING: 513# 10000 items take a very long time! 514# 515if { 0 } { 516 proc cost {params} { 517 set cost 0 518 foreach p [lrange $params 0 9] { 519 if { $p == 0 } { 520 incr cost 521 } 522 } 523 foreach p [lrange $params 10 end] { 524 if { $p == 1 } { 525 incr cost 526 } 527 } 528 return $cost 529 } 530 531 foreach n {100 1000 10000} { 532 break 533 puts "Problem size: $n" 534 puts [::simulation::annealing::findCombinatorialMinimum \ 535 -trials 300 \ 536 -verbose 0 \ 537 -number-params $n \ 538 -code {set result [cost $params]}] 539 } 540 541 # 542 # Second problem: 543 # Only the values of the first 10 items are important - 544 # they should be 1 545 # 546 proc cost2 {params} { 547 set cost 0 548 foreach p [lrange $params 0 9] { 549 if { $p == 0 } { 550 incr cost 551 } 552 } 553 return $cost 554 } 555 556 foreach n {100 1000 10000} { 557 puts "Problem size: $n" 558 puts [::simulation::annealing::findCombinatorialMinimum \ 559 -trials 300 \ 560 -verbose 0 \ 561 -number-params $n \ 562 -code {set result [cost2 $params]}] 563 } 564} 565