Optimizing model output

This is a simple example of optimizing output of a Simile model using SimileScripting and the Tcllib optimization package math::optimize. The script illustrates how to tackle this kind of problem.

"The package math::optimize is fully implemented in Tcl and no particular attention has been paid to the accuracy of the calculations. Instead, the algorithms have been used in a straightforward manner."

The script uses the procedure ::math::optimize::min_bound_1d that miminizes a function of one variable in the given interval. The procedure uses Brent's method of parabolic interpolation, protected by golden-section subdivisions if the interpolation is not converging. No guarantee is made that a global minimum is found. In the script the proceedure func is defined for ::math::optimize::min_bound_1d to evaluate. Proc func wraps up some SimileScript code that sets the input parameter value, resets the model and then runs the models. Finally it returns the value of a Simile model value for procedure ::math::optimize::min_bound_1d to minimize.

The package math::optimize also contains the procedure ::math::optimize::nelderMead that minimizes, in unconstrained fashion, a function of several variable over all of space, which may be worth experimenting with.

The model

Model mathOptimize
			Enumerated types: []
		 
		Variable   x
			x = 		Variable parameter 
			Minimum = 0, Maximum = 10
		 
		Variable   y
			y = 		x^3* exp(-3.0*x)

The response function of the model is shown below. To produce the plot the x values were set to the current simulation time.

Plot of the model response function.

The script

This script optimises the model output at the end of a run varying one input parameter. The script depends on math::optimize and uses the Plotchart package for a simple plot of the function output for a regular series of inputs the easiest way to get these package is probably to install ActiveTcl.

http://www.activestate.com/products/activetcl/

# This works to optimise for the last value of models runs.
	# The model is run for a set period (SetExecuteFor) and a parameter is set 
	# to a given value at the start of the run.
	# 
	
	# this script depends on math::optimize and uses the Plotchart for a
	# simple plot of the function output for a regular series of inputs
	# the easiest way to get these package is probably to install ActiveTcl
	# http://www.activestate.com/products/activetcl/
	
	
# Where the extra libraries live MODIFY FOR YOUR CASE
lappend auto_path {c:/program files/tcl/lib}; # eg ActiveTcl directories
package require math::optimize package require Plotchart # Simulistics provided package package require SimileAutoObj; # not actually needed in Windows using SimileScript.exe # Create a Tcl object/command (modelWin) with which to control that instance of Simile similescript::ModelWindow modelWin # Tell Simile not to use the single window Model Run Environment # We are not going to use any helper anyway here. # We then load the model and build it using C++ modelWin UseMRE false # MODIFY FOR YOUR CASE MUST BE FORWARD SLASHES EVEN UNDER WINDOWS!
modelWin Open "H:/mathOptimize.sml"
modelWin Run # Create a runControl command/object with which to control # (as you might expect) the run control. similescript::RunControl runControl runControl SetExecuteFor 1 runControl SetTimeStep 1 1 # Create proc to run model and return value of variable to optimize proc func {x} { runControl SetValue /x $x runControl Reset runControl Start # multiply the value by -1 as we want the maximum value using a minimising optimiser set y [runControl GetValue /y] return [expr -1*$y] } package require Plotchart toplevel .t canvas .t.c -background white -width 400 -height 200 pack .t.c -fill both # # Create the plot with its x- and y-axes # set s [::Plotchart::createXYPlot .t.c {0.0 20.0 2.0} {0.0 0.05 0.01}] $s title "Shape of function we are evaluating" # evalute the model output for a range of input values () # by plotting the output value for a series of input values # # multiply the value by -1 to get the original sign puts "Function shape" for {set i 0} {$i<=20} {incr i } { set x [expr {$i/2.0}] set y [expr {-1*[func $x]}] puts "$x $y" $s plot series1 $x $y } #puts "Maximum at(negative y value): [::math::optimize::min_bound_1d func 0.0 10.0]" set results [::math::optimize::min_bound_1d func 0.0 10.0 -trace on] set x [lindex $results 0] set y [expr {-1*[lindex $results 1]}] puts ##################################################### puts {::math::optimize::min_bound_1d func 0.0 10.0} puts "Maximum at: $x = $y" puts ##################################################### puts "" set results [::math::optimize::min_unbound_1d func 0.2 0.7 -trace on] set x [lindex $results 0] set y [expr {-1*[lindex $results 1]}] puts {::math::optimize::min_unbound_1d func 0.2 0.7} puts ##################################################### puts "Maximum at: $x = $y" puts #####################################################

Results

You can see its search for the minimum if the trace is on:

(Run) 159 % ::math::optimize::min_bound_1d func 0.0 10.0 -trace on
	f(3.8196601125) = -0.000588169133818 (initialisation)
	f(6.1803398875)=-2.09302565156e-006 (golden section)
	f(2.360679775)=-0.0110514849929 (golden section)
	f(1.4589803375)=-0.039018178181 (golden section)
	f(0.901699437495)=-0.0490203182137 (golden section)
	f(0.178385003365)=-0.0033239934159 (parabolic interpolation)
	f(0.625417908211)=-0.0374683267348 (golden section)
	f(1.05516631149)=-0.0495683007526 (parabolic interpolation)
	f(1.01193600577)=-0.0497765137855 (parabolic interpolation)
	f(1.00189882668)=-0.0497867994449 (parabolic interpolation)
	f(0.99971913441)=-0.0497870624755 (parabolic interpolation)
	f(0.999993651474)=-0.0497870683649 (parabolic interpolation)
	f(1.00000018344)=-0.0497870683679 (parabolic interpolation)
	f(0.999999991611)=-0.0497870683679 (parabolic interpolation)
	f(0.999999891511)=-0.0497870683679 (parabolic interpolation)
	f(0.999997508029)=-0.0497870683674 (golden section)
	f(0.999998981102)=-0.0497870683678 (golden section)
	f(0.999999543765)=-0.0497870683678 (golden section)
	f(0.999999758684)=-0.0497870683679 (golden section)
	f(0.999999658584)=-0.0497870683679 (golden section)
	0.999999658584 -0.0497870683679

Maximum at: 0.999999658584 = 0.0497870683679

Files