83d92d72c4c2dd2f4d335170225e6d820e8664f2
[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="example"><![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>A simple example</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     </refsection>
304     <refsection>
305         <title>Configuring a neighbour function</title>
306         <para>
307             In the following example, we customize the
308             neighbourhood function. In order to pass this function to the
309             <literal>optim_sa</literal> function, we setup a parameter where the
310             <literal>"neigh_func"</literal> key is associated with our particular neighbour function.
311             The neighbour function can be customized at will, provided that the
312             header of the function is the same. The particular implementation shown
313             below is the same, in spirit, as the <literal>neigh_func_default</literal>
314             function.
315         </para>
316         <programlisting role="example"><![CDATA[ 
317 function f = quad ( x )
318   p = [4 3];
319   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
320 endfunction
321
322 // We produce a neighbor by adding some noise to each component of a given vector
323 function x_neigh = myneigh_func ( x_current, T , param)
324   nxrow = size(x_current,"r")
325   nxcol = size(x_current,"c")
326   sa_min_delta = -0.1*ones(nxrow,nxcol);
327   sa_max_delta = 0.1*ones(nxrow,nxcol);
328   x_neigh = x_current + (sa_max_delta - sa_min_delta).*rand(nxrow,nxcol) + sa_min_delta;
329 endfunction
330
331 x0          = [2 2];
332 Proba_start = 0.7;
333 It_Pre      = 100;
334 It_extern   = 50;
335 It_intern   = 100;
336
337 saparams = init_param();
338 saparams = add_param(saparams,"neigh_func", myneigh_func);
339 // or: saparams = add_param(saparams,"neigh_func", neigh_func_default);
340 // or: saparams = add_param(saparams,"neigh_func", neigh_func_csa);
341 // or: saparams = add_param(saparams,"neigh_func", neigh_func_fsa);
342 // or: saparams = add_param(saparams,"neigh_func", neigh_func_vfsa);
343
344 T0 = compute_initial_temp(x0, quad, Proba_start, It_Pre, saparams);
345 Log = %f;
346 // This should produce x_opt = [4 3]
347 [x_opt, f_opt] = optim_sa(x0, quad, It_extern, It_intern, T0, Log,saparams);
348  ]]></programlisting>
349     </refsection>
350     <refsection>
351         <title>Passing extra parameters</title>
352         <para>
353             In the following example, we use an objective function which requires
354             an extra parameter <literal>p</literal>. This parameter is the second 
355             input argument of the <literal>quadp</literal> function. In order to 
356             pass this parameter to the objective function, we define the objective 
357             function as <literal>list(quadp,p)</literal>. In this case,
358             the solver makes so that the calling sequence includes a second argument.
359         </para>
360         <programlisting role="example"><![CDATA[ 
361   function f = quadp ( x , p )
362   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
363   endfunction
364
365   x0 = [-1 -1];
366   p = [4 3];
367   T0 = compute_initial_temp(x0, list(quadp,p) , Proba_start, It_Pre);
368   [x_opt, f_opt] = optim_sa(x0, list(quadp,p) , 10, 1000, T0, %f);
369  ]]></programlisting>
370     </refsection>
371     <refsection>
372         <title>Configuring an output function</title>
373         <para>
374             In the following example, we define an output function, which also 
375             provide a stopping rule. We define the function <literal>outfun</literal>
376             which takes as input arguments the data of the algorithm at the current 
377             iteration and returns the boolean <literal>stop</literal>. This function
378             prints a message into the console to inform the user about the 
379             current state of the algorithm. It also computes the boolean <literal>stop</literal>,
380             depending on the value of the function.
381             The stop variable becomes true when the function value is near zero. In order to let <literal>optim_sa</literal>
382             know about our output function, we configure the <literal>"output_func"</literal>
383             key to our <literal>outfun</literal> function and call the solver.
384             Notice that the number of external iterations is <literal>%inf</literal>, so
385             that the external loop never stops.
386             This allows to check that the output function really allows to
387             stop the algorithm.
388         </para>
389         <programlisting role="example"><![CDATA[ 
390 function f = quad ( x )
391   p = [4 3];
392   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
393 endfunction
394
395 function stop = outfunc ( itExt , x_best , f_best , T , saparams )
396   [mythreshold,err] = get_param(saparams,"mythreshold",0);
397   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 )
398   stop = ( abs(f_best) < mythreshold )
399 endfunction
400
401 x0 = [-1 -1];
402 saparams = init_param();
403 saparams = add_param(saparams,"output_func", outfunc );
404 saparams = add_param(saparams,"mythreshold", 1.e-2 );
405
406 rand("seed",0);
407
408 T0 = compute_initial_temp(x0, quad , 0.7, 100, saparams);
409 [x_best, f_best, mean_list, var_list, temp_list, f_history, x_history , It ] = optim_sa(x0, quad , %inf, 100, T0, %f, saparams);
410  ]]></programlisting>
411         <para>
412             The previous script produces the following output. Notice that the actual 
413             output of the algorithm depends on the state of the random number generator <literal>rand</literal>:
414             if we had not initialize the seed of the uniform random number generator, 
415             we would have produced a different result.
416         </para>
417         <programlisting role="example">
418             Iter = #1,   x_best=[-1.000000 -1.000000],   f_best = 4.100000e+001,         T = 1.453537e+000
419             Iter = #2,   x_best=[-0.408041 -0.318262],   f_best = 3.044169e+001,         T = 1.308183e+000
420             Iter = #3,   x_best=[-0.231406 -0.481078],   f_best = 3.002270e+001,         T = 1.177365e+000
421             Iter = #4,   x_best=[0.661827 0.083743],     f_best = 1.964796e+001,         T = 1.059628e+000
422             Iter = #5,   x_best=[0.931415 0.820681],     f_best = 1.416565e+001,         T = 9.536654e-001
423             Iter = #6,   x_best=[1.849796 1.222178],     f_best = 7.784028e+000,         T = 8.582988e-001
424             Iter = #7,   x_best=[2.539775 1.414591],     f_best = 4.645780e+000,         T = 7.724690e-001
425             Iter = #8,   x_best=[3.206047 2.394497],     f_best = 9.969957e-001,         T = 6.952221e-001
426             Iter = #9,   x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 6.256999e-001
427             Iter = #10,          x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 5.631299e-001
428             Iter = #11,          x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 5.068169e-001
429             Iter = #12,          x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 4.561352e-001
430             Iter = #13,          x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 4.105217e-001
431             Iter = #14,          x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 3.694695e-001
432             Iter = #15,          x_best=[3.931929 3.003181],     f_best = 4.643767e-003,         T = 3.325226e-001
433         </programlisting>
434     </refsection>
435     <refsection role="see also">
436         <title>See Also</title>
437         <simplelist type="inline">
438             <member>
439                 <link linkend="compute_initial_temp">
440                     compute_initial_temp
441                 </link>
442             </member>
443             <member>
444                 <link linkend="neigh_func_default">
445                     neigh_func_default
446                 </link>
447             </member>
448             <member>
449                 <link linkend="temp_law_default">
450                     temp_law_default
451                 </link>
452             </member>
453         </simplelist>
454     </refsection>
455     <refsection>
456         <title>Bibliography</title>
457         <para>
458             "Simulated annealing : theory and applications", P.J.M. Laarhoven and E.H.L. Aarts, Mathematics and its applications, Dordrecht : D. Reidel, 1988
459         </para>
460         <para>
461             "Theoretical and computational aspects of simulated annealing", P.J.M. van Laarhoven, Amsterdam, Netherlands : Centrum voor Wiskunde en Informatica, 1988
462         </para>
463         <para>
464             "Genetic algorithms and simulated annealing", Lawrence Davis, London : Pitman Los Altos, Calif. Morgan Kaufmann Publishers, 1987
465         </para>
466     </refsection>
467 </refentry>