Examples for simulated annealing functions added and/or improved
[scilab.git] / scilab / modules / simulated_annealing / help / en_US / algorithms / optim_sa.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>
5  * Copyright (C) 2010 - Michael Baudin
6  *
7  * This file must be used under the terms of the CeCILL.
8  * This source file is licensed as described in the file COPYING, which
9  * you should have received as part of this distribution.  The terms
10  * are also available at
11  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
12  *
13  -->
14 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="optim_sa" xml:lang="en">
15     <refnamediv>
16         <refname>optim_sa</refname>
17         <refpurpose>A Simulated Annealing optimization method</refpurpose>
18     </refnamediv>
19     <refsynopsisdiv>
20         <title>Calling Sequence</title>
21         <synopsis>
22             x_best = optim_sa(x0,f,ItExt,ItInt,T0,Log,temp_law,param_temp_law,neigh_func,param_neigh_func)
23             [x_best,f_best] = optim_sa(..)
24             [x_best,f_best,mean_list] = optim_sa(..)
25             [x_best,f_best,mean_list,var_list] = optim_sa(..)
26             [x_best,f_best,mean_list,var_list,f_history] = optim_sa(..)
27             [x_best,f_best,mean_list,var_list,f_history,temp_list] = optim_sa(..)
28             [x_best,f_best,mean_list,var_list,f_history,temp_list,x_history] = optim_sa(..)
29             [x_best,f_best,mean_list,var_list,f_history,temp_list,x_history,iter] = optim_sa(..)
30         </synopsis>
31     </refsynopsisdiv>
32     <refsection>
33         <title>Arguments</title>
34         <variablelist>
35             <varlistentry>
36                 <term>x0</term>
37                 <listitem>
38                     <para>the initial solution</para>
39                 </listitem>
40             </varlistentry>
41             <varlistentry>
42                 <term>f</term>
43                 <listitem>
44                     <para>
45                         the objective function to be optimized (the prototype if
46                         f(x))
47                     </para>
48                 </listitem>
49             </varlistentry>
50             <varlistentry>
51                 <term>ItExt</term>
52                 <listitem>
53                     <para>the number of temperature decrease</para>
54                 </listitem>
55             </varlistentry>
56             <varlistentry>
57                 <term>ItInt</term>
58                 <listitem>
59                     <para>the number of iterations during one temperature stage</para>
60                 </listitem>
61             </varlistentry>
62             <varlistentry>
63                 <term>T0</term>
64                 <listitem>
65                     <para>
66                         the initial temperature (see compute_initial_temp to compute
67                         easily this temperature)
68                     </para>
69                 </listitem>
70             </varlistentry>
71             <varlistentry>
72                 <term>Log</term>
73                 <listitem>
74                     <para>
75                         if %T, some information will be displayed during the run of
76                         the simulated annealing
77                     </para>
78                 </listitem>
79             </varlistentry>
80             <varlistentry>
81                 <term>temp_law</term>
82                 <listitem>
83                     <para>
84                         the temperature decrease law (see temp_law_default for an
85                         example of such a function)
86                     </para>
87                 </listitem>
88             </varlistentry>
89             <varlistentry>
90                 <term>param_temp_law</term>
91                 <listitem>
92                     <para>
93                         a structure (of any kind - it depends on the temperature law
94                         used) which is transmitted as a parameter to temp_law
95                     </para>
96                 </listitem>
97             </varlistentry>
98             <varlistentry>
99                 <term>neigh_func</term>
100                 <listitem>
101                     <para>
102                         a function which computes a neighbor of a given point (see
103                         neigh_func_default for an example of such a function)
104                     </para>
105                 </listitem>
106             </varlistentry>
107             <varlistentry>
108                 <term>param_neigh_func</term>
109                 <listitem>
110                     <para>
111                         a structure (of any kind like vector, list, it depends on the
112                         neighborhood function used) which is transmitted as a parameter to
113                         neigh_func
114                     </para>
115                 </listitem>
116             </varlistentry>
117             <varlistentry>
118                 <term>x_best</term>
119                 <listitem>
120                     <para>the best solution found so far</para>
121                 </listitem>
122             </varlistentry>
123             <varlistentry>
124                 <term>f_best</term>
125                 <listitem>
126                     <para>the objective function value corresponding to x_best</para>
127                 </listitem>
128             </varlistentry>
129             <varlistentry>
130                 <term>mean_list</term>
131                 <listitem>
132                     <para>
133                         the mean of the objective function value for each temperature
134                         stage. A vector of float (optional)
135                     </para>
136                 </listitem>
137             </varlistentry>
138             <varlistentry>
139                 <term>var_list</term>
140                 <listitem>
141                     <para>
142                         the variance of the objective function values for each
143                         temperature stage. A vector of float (optional)
144                     </para>
145                 </listitem>
146             </varlistentry>
147             <varlistentry>
148                 <term>f_history</term>
149                 <listitem>
150                     <para>
151                         the computed objective function values for each iteration.
152                         Each input of the list corresponds to a temperature stage. Each
153                         input of the list is a vector of float which gathers all the
154                         objective function values computed during the corresponding
155                         temperature stage - (optional)
156                     </para>
157                 </listitem>
158             </varlistentry>
159             <varlistentry>
160                 <term>temp_list</term>
161                 <listitem>
162                     <para>
163                         the list of temperature computed for each temperature stage. A
164                         vector of float (optional)
165                     </para>
166                 </listitem>
167             </varlistentry>
168             <varlistentry>
169                 <term>x_history</term>
170                 <listitem>
171                     <para>
172                         the parameter values computed for each iteration. Each input
173                         of the list corresponds to a temperature stage. Each input of the
174                         list is a vector of input variables which corresponds to all the
175                         variables computed during the corresponding temperature stage -
176                         (optional - can slow down a lot the execution of optim_sa)
177                     </para>
178                 </listitem>
179             </varlistentry>
180             <varlistentry>
181                 <term>iter</term>
182                 <listitem>
183                     <para>
184                         a double, the actual number of external iterations in the
185                         algorithm (optional).
186                     </para>
187                 </listitem>
188             </varlistentry>
189         </variablelist>
190     </refsection>
191     <refsection>
192         <title>Description</title>
193         <para>A Simulated Annealing optimization method.</para>
194         <para>
195             Simulated annealing (SA) is a generic probabilistic meta-algorithm for the global optimization
196             problem, namely locating a good approximation to the global optimum of a given function in a
197             large search space. It is often used when the search space is discrete (e.g., all tours that
198             visit a given set of cities).
199         </para>
200         <para>
201             The current solver can find the solution of an optimization problem without constraints or 
202             with bound constraints. The bound constraints can be customized with the neighbour
203             function. This algorithm does not use the derivatives of the objective function.
204         </para>
205         <para>
206             The solver is made of Scilab macros, which enables a high-level programming model for
207             this optimization solver. The SA macros are based on the <literal>parameters</literal>
208             Scilab module for the management of the (many) optional parameters.
209         </para>
210         <para>
211             To use the SA algorithm, one should perform the following steps :
212             <itemizedlist><listitem>
213                     configure the parameters with calls to <literal>init_param</literal> 
214                     and <literal>add_param</literal> especially the neighbor function, the 
215                     acceptance function, the temperature law,
216                 </listitem>
217                 <listitem>
218                     compute an initial temperature with a call to <literal>compute_initial_temp</literal>,
219                 </listitem>
220                 <listitem>
221                     find an optimum by using the <literal>optim_sa</literal> solver.
222                 </listitem>
223             </itemizedlist>
224         </para>
225         <para>
226             The algorithm is based on an iterative update of two points :
227             <itemizedlist><listitem>
228                     the current point is updated by taking into account the neighbour and the acceptance
229                     functions,
230                 </listitem>
231                 <listitem>
232                     the best point is the point which achieved the minimum of the objective function over the
233                     iterations.
234                 </listitem>
235             </itemizedlist>
236             While the current point is used internally to explore the domain, only the best point is returned
237             by the function.
238             The algorithm is based on an external loop and an internal loop. In the external loop,
239             the temperature is updated according to the temperature function. In the internal loop, the 
240             point is updated according to the neighbour function. A new point is accepted depending 
241             on its associated function value or the value of the acceptance function, which value 
242             depends on the current temperature and a uniform random number.
243         </para>
244         <para>
245             The acceptance of the new point depends on the output values produced
246             by the <literal>rand</literal> function. This implies that two consecutive 
247             calls to the <literal>optim_sa</literal> will not produce the same result.
248             In order to always get exactly the same results, please initialize the random number 
249             generator with a valid seed.
250         </para>
251         <para>
252             See the Demonstrations, in the "Optimization" section and "Simulated Annealing" subsection
253             for more examples.
254         </para>
255     </refsection>
256     <refsection>
257         <title>The objective function</title>
258         <para>
259             The objective function is expected to have the following header.
260         </para>
261         <programlisting role="no-scilab-exec"><![CDATA[ 
262 function y = f ( x )
263  ]]></programlisting>
264         <para>
265             In the case where the objective function needs additionnal parameters,
266             the objective function can be defined as a list, where the first 
267             argument is the cost function, and the second argument is the 
268             additionnal parameter. See below for an example.
269         </para>
270     </refsection>
271     <refsection>
272         <title>Examples</title>
273         <para>
274             In the following example, we search the minimum of the
275             Rastriging function. This function has many local minimas, but only
276             one single global minimum located at x = (0,0), where the function value is
277             f(x) = -2. We use the simulated annealing algorithm with default settings
278             and the default neighbour function neigh_func_default.
279         </para>
280         <programlisting role="example"><![CDATA[ 
281 function y = rastrigin(x)
282   y = x(1)^2+x(2)^2-cos(12*x(1))-cos(18*x(2));
283 endfunction
284     
285 x0          = [2 2];
286 Proba_start = 0.7;
287 It_Pre      = 100;
288 It_extern   = 100;
289 It_intern   = 1000;
290 x_test = neigh_func_default(x0);
291
292 T0 = compute_initial_temp(x0, rastrigin, Proba_start, It_Pre);
293
294 Log = %T;
295 [x_opt, f_opt, sa_mean_list, sa_var_list] = optim_sa(x0, rastrigin, It_extern, It_intern, T0, Log);
296
297 mprintf("optimal solution:\n"); disp(x_opt);
298 mprintf("value of the objective function = %f\n", f_opt);
299
300 t = 1:length(sa_mean_list);
301 plot(t,sa_mean_list,"r",t,sa_var_list,"g");
302  ]]></programlisting>
303         <scilab:image>
304 function y = rastrigin(x)
305   y = x(1)^2+x(2)^2-cos(12*x(1))-cos(18*x(2));
306 endfunction
307     
308 x0          = [2 2];
309 Proba_start = 0.7;
310 It_Pre      = 100;
311 It_extern   = 100;
312 It_intern   = 1000;
313 x_test = neigh_func_default(x0);
314
315 T0 = compute_initial_temp(x0, rastrigin, Proba_start, It_Pre);
316
317 Log = %T;
318 [x_opt, f_opt, sa_mean_list, sa_var_list] = optim_sa(x0, rastrigin, It_extern, It_intern, T0, Log);
319
320 t = 1:length(sa_mean_list);
321 plot(t,sa_mean_list,"r",t,sa_var_list,"g");
322         </scilab:image>
323     </refsection>
324     <refsection>
325         <title>Configuring a neighbour function</title>
326         <para>
327             In the following example, we customize the
328             neighbourhood function. In order to pass this function to the
329             <literal>optim_sa</literal> function, we setup a parameter where the
330             <literal>"neigh_func"</literal> key is associated with our particular neighbour function.
331             The neighbour function can be customized at will, provided that the
332             header of the function is the same. The particular implementation shown
333             below is the same, in spirit, as the <literal>neigh_func_default</literal>
334             function.
335         </para>
336         <programlisting role="example"><![CDATA[ 
337 function f = quad ( x )
338   p = [4 3];
339   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
340 endfunction
341
342 // We produce a neighbor by adding some noise to each component of a given vector
343 function x_neigh = myneigh_func ( x_current, T , param)
344   nxrow = size(x_current,"r")
345   nxcol = size(x_current,"c")
346   sa_min_delta = -0.1*ones(nxrow,nxcol);
347   sa_max_delta = 0.1*ones(nxrow,nxcol);
348   x_neigh = x_current + (sa_max_delta - sa_min_delta).*rand(nxrow,nxcol) + sa_min_delta;
349 endfunction
350
351 x0          = [2 2];
352 Proba_start = 0.7;
353 It_Pre      = 100;
354 It_extern   = 50;
355 It_intern   = 100;
356
357 saparams = init_param();
358 saparams = add_param(saparams,"neigh_func", myneigh_func);
359 // or: saparams = add_param(saparams,"neigh_func", neigh_func_default);
360 // or: saparams = add_param(saparams,"neigh_func", neigh_func_csa);
361 // or: saparams = add_param(saparams,"neigh_func", neigh_func_fsa);
362 // or: saparams = add_param(saparams,"neigh_func", neigh_func_vfsa);
363
364 T0 = compute_initial_temp(x0, quad, Proba_start, It_Pre, saparams);
365 Log = %f;
366 // This should produce x_opt = [4 3]
367 [x_opt, f_opt] = optim_sa(x0, quad, It_extern, It_intern, T0, Log,saparams)
368  ]]></programlisting>
369     </refsection>
370     <refsection>
371         <title>Passing extra parameters</title>
372         <para>
373             In the following example, we use an objective function which requires
374             an extra parameter <literal>p</literal>. This parameter is the second 
375             input argument of the <literal>quadp</literal> function. In order to 
376             pass this parameter to the objective function, we define the objective 
377             function as <literal>list(quadp,p)</literal>. In this case,
378             the solver makes so that the calling sequence includes a second argument.
379         </para>
380         <programlisting role="example"><![CDATA[ 
381   function f = quadp ( x , p )
382   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
383   endfunction
384
385   x0 = [-1 -1];
386   p = [4 3];
387   Proba_start = 0.7;
388   It_Pre      = 100;
389   T0 = compute_initial_temp(x0, list(quadp,p) , Proba_start, It_Pre);
390   [x_opt, f_opt] = optim_sa(x0, list(quadp,p) , 10, 1000, T0, %f)
391  ]]></programlisting>
392     </refsection>
393     <refsection>
394         <title>Configuring an output function</title>
395         <para>
396             In the following example, we define an output function, which also 
397             provide a stopping rule. We define the function <literal>outfun</literal>
398             which takes as input arguments the data of the algorithm at the current 
399             iteration and returns the boolean <literal>stop</literal>. This function
400             prints a message into the console to inform the user about the 
401             current state of the algorithm. It also computes the boolean <literal>stop</literal>,
402             depending on the value of the function.
403             The stop variable becomes true when the function value is near zero. In order to let <literal>optim_sa</literal>
404             know about our output function, we configure the <literal>"output_func"</literal>
405             key to our <literal>outfun</literal> function and call the solver.
406             Notice that the number of external iterations is <literal>%inf</literal>, so
407             that the external loop never stops.
408             This allows to check that the output function really allows to
409             stop the algorithm.
410         </para>
411         <programlisting role="example"><![CDATA[ 
412 function f = quad ( x )
413   p = [4 3];
414   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
415 endfunction
416
417 function stop = outfunc ( itExt , x_best , f_best , T , saparams )
418   [mythreshold,err] = get_param(saparams,"mythreshold",0);
419   mprintf ( "Iter = #%d, \t x_best=[%f %f], \t f_best = %e, \t T = %e\n", itExt , x_best(1),x_best(2) , f_best , T )
420   stop = ( abs(f_best) < mythreshold )
421 endfunction
422
423 x0 = [-1 -1];
424 saparams = init_param();
425 saparams = add_param(saparams,"output_func", outfunc );
426 saparams = add_param(saparams,"mythreshold", 1.e-2 );
427
428 rand("seed",0);
429
430 T0 = compute_initial_temp(x0, quad , 0.7, 100, saparams);
431 [x_best, f_best, mean_list, var_list, temp_list, f_history, x_history , It ] = optim_sa(x0, quad , %inf, 100, T0, %f, saparams);
432  ]]></programlisting>
433         <para>
434             The previous script produces the following output. Notice that the actual 
435             output of the algorithm depends on the state of the random number generator <literal>rand</literal>:
436             if we had not initialize the seed of the uniform random number generator, 
437             we would have produced a different result.
438         </para>
439         <programlisting role="no-scilab-exec">
440             Iter = #1,   x_best=[-1.000000 -1.000000],   f_best = 4.100000e+001,         T = 1.453537e+000
441             Iter = #2,   x_best=[-0.408041 -0.318262],   f_best = 3.044169e+001,         T = 1.308183e+000
442             Iter = #3,   x_best=[-0.231406 -0.481078],   f_best = 3.002270e+001,         T = 1.177365e+000
443             Iter = #4,   x_best=[0.661827 0.083743],     f_best = 1.964796e+001,         T = 1.059628e+000
444             Iter = #5,   x_best=[0.931415 0.820681],     f_best = 1.416565e+001,         T = 9.536654e-001
445             Iter = #6,   x_best=[1.849796 1.222178],     f_best = 7.784028e+000,         T = 8.582988e-001
446             Iter = #7,   x_best=[2.539775 1.414591],     f_best = 4.645780e+000,         T = 7.724690e-001
447             Iter = #8,   x_best=[3.206047 2.394497],     f_best = 9.969957e-001,         T = 6.952221e-001
448             Iter = #9,   x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 6.256999e-001
449             Iter = #10,          x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 5.631299e-001
450             Iter = #11,          x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 5.068169e-001
451             Iter = #12,          x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 4.561352e-001
452             Iter = #13,          x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 4.105217e-001
453             Iter = #14,          x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 3.694695e-001
454             Iter = #15,          x_best=[3.931929 3.003181],     f_best = 4.643767e-003,         T = 3.325226e-001
455         </programlisting>
456     </refsection>
457     <refsection role="see also">
458         <title>See Also</title>
459         <simplelist type="inline">
460             <member>
461                 <link linkend="compute_initial_temp">
462                     compute_initial_temp
463                 </link>
464             </member>
465             <member>
466                 <link linkend="neigh_func_default">
467                     neigh_func_default
468                 </link>
469             </member>
470             <member>
471                 <link linkend="temp_law_default">
472                     temp_law_default
473                 </link>
474             </member>
475         </simplelist>
476     </refsection>
477     <refsection>
478         <title>Bibliography</title>
479         <para>
480             "Simulated annealing : theory and applications", P.J.M. Laarhoven and E.H.L. Aarts, Mathematics and its applications, Dordrecht : D. Reidel, 1988
481         </para>
482         <para>
483             "Theoretical and computational aspects of simulated annealing", P.J.M. van Laarhoven, Amsterdam, Netherlands : Centrum voor Wiskunde en Informatica, 1988
484         </para>
485         <para>
486             "Genetic algorithms and simulated annealing", Lawrence Davis, London : Pitman Los Altos, Calif. Morgan Kaufmann Publishers, 1987
487         </para>
488     </refsection>
489 </refentry>