scilab/modules/optimization/help/en_US/optim.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 - INRIA
   5  *
   6  * This file must be used under the terms of the CeCILL.
   7  * This source file is licensed as described in the file COPYING, which
   8  * you should have received as part of this distribution.  The terms
   9  * are also available at
  10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  11  *
  12  -->
  13 <refentry version="5.0-subset Scilab" xml:id="optim" xml:lang="en"
  14           xmlns="http://docbook.org/ns/docbook"
  15           xmlns:xlink="http://www.w3.org/1999/xlink"
  16           xmlns:svg="http://www.w3.org/2000/svg"
  17           xmlns:ns3="http://www.w3.org/1999/xhtml"
  18           xmlns:mml="http://www.w3.org/1998/Math/MathML"
  19           xmlns:db="http://docbook.org/ns/docbook">
  20   <info>
  21     <pubdate>$LastChangedDate: 2008-04-28 09:36:26 +0200 (lun, 28 avr 2008)
  22     $</pubdate>
  23   </info>
  24
  25   <refnamediv>
  26     <refname>optim</refname>
  27
  28     <refpurpose>non-linear optimization routine</refpurpose>
  29   </refnamediv>
  30
  31   <refsynopsisdiv>
  32     <title>Calling Sequence</title>
  33
  34     <synopsis>[f,xopt]=optim(costf,x0)
  35 [f [,xopt [,gradopt [,work]]]]=optim(costf [,&lt;contr&gt;],x0 [,algo] [,df0 [,mem]] [,work] [,&lt;stop&gt;] [,&lt;params&gt;] [,imp=iflag])</synopsis>
  36   </refsynopsisdiv>
  37
  38   <refsection>
  39     <title>Parameters</title>
  40
  41     <variablelist>
  42       <varlistentry>
  43         <term>costf</term>
  44
  45         <listitem>
  46           <para>external, i.e Scilab function list or string
  47           (<literal>costf</literal> is the cost function, that is, a Scilab
  48           script, a Fortran 77 routine or a C function.</para>
  49         </listitem>
  50       </varlistentry>
  51
  52       <varlistentry>
  53         <term>x0</term>
  54
  55         <listitem>
  56           <para>real vector (initial value of variable to be
  57           minimized).</para>
  58         </listitem>
  59       </varlistentry>
  60
  61       <varlistentry>
  62         <term>f</term>
  63
  64         <listitem>
  65           <para>value of optimal cost
  66           (<literal>f=costf(xopt)</literal>)</para>
  67         </listitem>
  68       </varlistentry>
  69
  70       <varlistentry>
  71         <term>xopt</term>
  72
  73         <listitem>
  74           <para>best value of <literal>x</literal> found.</para>
  75         </listitem>
  76       </varlistentry>
  77
  78       <varlistentry>
  79         <term>&lt;contr&gt;</term>
  80
  81         <listitem>
  82           <para>keyword representing the following sequence of arguments:
  83           <literal>'b',binf,bsup</literal> with <literal>binf</literal> and
  84           <literal>bsup</literal> are real vectors with same dimension as
  85           <literal>x0</literal>. <literal>binf</literal> and
  86           <literal>bsup</literal> are lower and upper bounds on
  87           <literal>x</literal>.</para>
  88         </listitem>
  89       </varlistentry>
  90
  91       <varlistentry>
  92         <term>algo</term>
  93
  94         <listitem>
  95           <itemizedlist>
  96             <listitem>
  97               <para><literal>'qn'</literal> : quasi-Newton (this is the
  98               default solver)</para>
  99             </listitem>
 100
 101             <listitem>
 102               <para><literal>'gc'</literal> : conjugate gradient</para>
 103             </listitem>
 104
 105             <listitem>
 106               <para><literal>'nd'</literal> : non-differentiable.</para>
 107
 108               <para>Note that the conjugate gradient solver does not accept
 109               bounds on <literal>x</literal>.</para>
 110             </listitem>
 111           </itemizedlist>
 112         </listitem>
 113       </varlistentry>
 114
 115       <varlistentry>
 116         <term>df0</term>
 117
 118         <listitem>
 119           <para>real scalar. Guessed decreasing of <literal>f</literal> at
 120           first iteration. (<literal>df0=1</literal> is the default
 121           value).</para>
 122         </listitem>
 123       </varlistentry>
 124
 125       <varlistentry>
 126         <term>mem :</term>
 127
 128         <listitem>
 129           <para>integer, number of variables used to approximate the Hessian.
 130           Default value is 10. This feature is available for the
 131           Gradient-Conjugate algorithm "gc" without constraints and the
 132           non-smooth algorithm "nd" without constraints.</para>
 133         </listitem>
 134       </varlistentry>
 135
 136       <varlistentry>
 137         <term>&lt;stop&gt;</term>
 138
 139         <listitem>
 140           <para>keyword representing the sequence of optional parameters
 141           controlling the convergence of the algorithm. <literal>'ar',nap
 142           [,iter [,epsg [,epsf [,epsx]]]]</literal></para>
 143
 144           <variablelist>
 145             <varlistentry>
 146               <term>"ar"</term>
 147
 148               <listitem>
 149                 <para>reserved keyword for stopping rule selection defined as
 150                 follows:</para>
 151               </listitem>
 152             </varlistentry>
 153
 154             <varlistentry>
 155               <term>nap</term>
 156
 157               <listitem>
 158                 <para>maximum number of calls to <literal>costf</literal>
 159                 allowed (default is 100).</para>
 160               </listitem>
 161             </varlistentry>
 162
 163             <varlistentry>
 164               <term>iter</term>
 165
 166               <listitem>
 167                 <para>maximum number of iterations allowed (default is
 168                 100).</para>
 169               </listitem>
 170             </varlistentry>
 171
 172             <varlistentry>
 173               <term>epsg</term>
 174
 175               <listitem>
 176                 <para>threshold on gradient norm.</para>
 177               </listitem>
 178             </varlistentry>
 179
 180             <varlistentry>
 181               <term>epsf</term>
 182
 183               <listitem>
 184                 <para>threshold controlling decreasing of
 185                 <literal>f</literal></para>
 186               </listitem>
 187             </varlistentry>
 188
 189             <varlistentry>
 190               <term>epsx</term>
 191
 192               <listitem>
 193                 <para>threshold controlling variation of <literal>x</literal>.
 194                 This vector (possibly matrix) of same size as
 195                 <literal>x0</literal> can be used to scale
 196                 <literal>x</literal>.</para>
 197               </listitem>
 198             </varlistentry>
 199           </variablelist>
 200         </listitem>
 201       </varlistentry>
 202
 203       <varlistentry>
 204         <term>&lt;params&gt;</term>
 205
 206         <listitem>
 207           <para>keyword representing the method to initialize the arguments
 208           <literal>ti, td</literal> passed to the objective function, provided
 209           as a C or Fortran routine. This option has no meaning when the cost
 210           function is a Scilab script. &lt;params&gt; can be set to only one
 211           of the following values.</para>
 212
 213           <itemizedlist>
 214             <listitem>
 215               <para>"in"</para>
 216
 217               <para>That mode allows to allocate memory in the internal Scilab
 218               workspace so that the objective function can get arrays with the
 219               required size, but without directly allocating the memory. "in"
 220               stands for "initialization". In that mode, before the value and
 221               derivative of the objective function is to be computed, there is
 222               a dialog between the optim Scilab primitive and the objective
 223               function. In this dialog, the objective function is called two
 224               times, with particular values of the "ind" parameter. The first
 225               time, ind is set to 10 and the objective function is expected to
 226               set the nizs, nrzs and ndzs integer parameters of the "nird"
 227               common.</para>
 228
 229               <programlisting role = ""><![CDATA[
 230 common /nird/ nizs,nrzs,ndzs
 231  ]]></programlisting>
 232
 233               <para>This allows Scilab to allocate memory inside its internal
 234               workspace. The second time the objective function is called, ind
 235               is set to 11 and the objective function is expected to set the
 236               ti, tr and tz arrays. After this initialization phase, each time
 237               it is called, the objective function is ensured that the ti, tr
 238               and tz arrays which are passed to it have the values that have
 239               been previously initialized.</para>
 240             </listitem>
 241
 242             <listitem>
 243               <para>"ti",valti</para>
 244
 245               <para>In this mode, valti is expected to be a Scilab vector
 246               variable containing integers. Whenever the objective function is
 247               called, the ti array it receives contains the values of the
 248               Scilab variable.</para>
 249             </listitem>
 250
 251             <listitem>
 252               <para>"td", valtd</para>
 253
 254               <para>In this mode, valtd is expected to be a Scilab vector
 255               variable containing double values. Whenever the objective
 256               function is called, the td array it receives contains the values
 257               of the Scilab variable.</para>
 258             </listitem>
 259
 260             <listitem>
 261               <para>"ti",valti,"td",valtd</para>
 262
 263               <para>This mode combines the two previous.</para>
 264             </listitem>
 265           </itemizedlist>
 266
 267           <para>The <literal>ti, td</literal> arrays may be used so that the
 268           objective function can be computed. For example, if the objective
 269           function is a polynomial, the ti array may may be used to store the
 270           coefficients of that polynomial.</para>
 271
 272           <para>Users should choose carefully between the "in" mode and the
 273           "ti" and "td" mode, depending on the fact that the arrays are Scilab
 274           variables or not. If the data is available as Scilab variables, then
 275           the "ti", valti, "td", valtd mode should be chosen. If the data is
 276           available directly from the objective function, the "in" mode should
 277           be chosen. Notice that there is no "tr" mode, since, in Scilab, all
 278           real values are of "double" type.</para>
 279
 280           <para>If neither the "in" mode, nor the "ti", "td" mode is chosen,
 281           that is, if &lt;params&gt; is not present as an option of the optim
 282           primitive, the user may should not assume that the ti,tr and td
 283           arrays can be used : reading or writing the arrays may generate
 284           unpredictable results.</para>
 285         </listitem>
 286       </varlistentry>
 287
 288       <varlistentry>
 289         <term>"imp=iflag"</term>
 290
 291         <listitem>
 292           <para>named argument used to set the trace mode. The possible values
 293           for iflag are 0,1,2 and &gt;2. Use this option with caution : most
 294           of these reports are written on the Scilab standard output.</para>
 295
 296           <itemizedlist>
 297             <listitem>
 298               <para>iflag=0: nothing (except errors) is reported (this is the
 299               default),</para>
 300             </listitem>
 301
 302             <listitem>
 303               <para>iflag=1: initial and final reports,</para>
 304             </listitem>
 305
 306             <listitem>
 307               <para>iflag=2: adds a report per iteration,</para>
 308             </listitem>
 309
 310             <listitem>
 311               <para>iflag&gt;2: add reports on linear search.</para>
 312             </listitem>
 313
 314             <listitem>
 315               <para>iflag&lt;0: calls the cost function with ind=1 every -imp iterations.</para>
 316             </listitem>
 317           </itemizedlist>
 318         </listitem>
 319       </varlistentry>
 320
 321       <varlistentry>
 322         <term>gradopt</term>
 323
 324         <listitem>
 325           <para>gradient of <literal>costf</literal> at
 326           <literal>xopt</literal></para>
 327         </listitem>
 328       </varlistentry>
 329
 330       <varlistentry>
 331         <term>work</term>
 332
 333         <listitem>
 334           <para>working array for hot restart for quasi-Newton method. This
 335           array is automatically initialized by <literal>optim</literal> when
 336           <literal>optim</literal> is invoked. It can be used as input
 337           parameter to speed-up the calculations.</para>
 338         </listitem>
 339       </varlistentry>
 340     </variablelist>
 341   </refsection>
 342
 343   <refsection>
 344     <title>Description</title>
 345
 346     <para>Non-linear optimization routine for programs without constraints or
 347     with bound constraints:</para>
 348
 349     <programlisting role = ""><![CDATA[
 350 min costf(x) w.r.t x.
 351  ]]></programlisting>
 352
 353     <para><literal>costf</literal> is an "external" i.e a Scilab function, a
 354     list or a string giving the name of a C or Fortran routine (see
 355     "external"). This external must return the value <literal>f</literal> of
 356     the cost function at the point <literal>x</literal> and the gradient
 357     <literal>g</literal> of the cost function at the point
 358     <literal>x</literal>.</para>
 359
 360     <variablelist>
 361       <varlistentry>
 362         <term>- Scilab function case</term>
 363
 364         <listitem>
 365           <para>If <literal>costf</literal> is a Scilab function, the calling
 366           sequence for <literal>costf</literal> must be:</para>
 367
 368           <programlisting role = ""><![CDATA[
 369 [f,g,ind]=costf(x,ind)
 370  ]]></programlisting>
 371
 372           <para>Here, <literal>costf</literal> is a function which returns
 373           <literal>f</literal>, value (real number) of cost function at
 374           <literal>x</literal>, and <literal>g</literal>, gradient vector of
 375           cost function at <literal>x</literal>. The variable
 376           <literal>ind</literal> is described below.</para>
 377         </listitem>
 378       </varlistentry>
 379
 380       <varlistentry>
 381         <term>- List case</term>
 382
 383         <listitem>
 384           <para>If <literal>costf</literal> is a list, it should be of the
 385           form: <literal>list(real_costf, arg1,...,argn)</literal> with
 386           <literal>real_costf</literal> a Scilab function with calling
 387           sequence : <literal>[f,g,ind]=costf(x,ind,arg1,... argn)</literal>.
 388           The <literal>x</literal>, <literal>f</literal>,
 389           <literal>g</literal>, <literal>ind</literal> arguments have the same
 390           meaning that above. <literal>argi</literal> arguments can be used to
 391           pass function parameters.</para>
 392         </listitem>
 393       </varlistentry>
 394
 395       <varlistentry>
 396         <term>- String case</term>
 397
 398         <listitem>
 399           <para>If <literal>costf</literal> is a character string, it refers
 400           to the name of a C or Fortran routine which must be linked to
 401           Scilab</para>
 402
 403           <variablelist>
 404             <varlistentry>
 405               <term>* Fortran case</term>
 406
 407               <listitem>
 408                 <para>The interface of the Fortran subroutine computing the
 409                 objective must be :</para>
 410
 411                 <programlisting role = ""><![CDATA[
 412 subroutine costf(ind,n,x,f,g,ti,tr,td)
 413  ]]></programlisting>
 414
 415                 <para>with the following declarations:</para>
 416
 417                 <programlisting role = ""><![CDATA[
 418 integer ind,n ti(*)
 419 double precision x(n),f,g(n),td(*)
 420 real tr(*)
 421  ]]></programlisting>
 422
 423                 <para>The argument <literal>ind</literal> is described
 424                 below.</para>
 425
 426                 <para>If ind = 2, 3 or 4, the inputs of the routine are :
 427                 <literal>x, ind, n, ti, tr,td</literal>.</para>
 428
 429                 <para>If ind = 2, 3 or 4, the outputs of the routine are :
 430                 <literal>f</literal> and <literal>g</literal>.</para>
 431               </listitem>
 432             </varlistentry>
 433
 434             <varlistentry>
 435               <term>* C case</term>
 436
 437               <listitem>
 438                 <para>The interface of the C function computing the objective
 439                 must be :</para>
 440
 441                 <programlisting role = ""><![CDATA[
 442 void costf(int *ind, int *n, double *x, double *f, double *g, int *ti, float *tr, double *td)
 443  ]]></programlisting>
 444
 445                 <para>The argument <literal>ind</literal> is described
 446                 below.</para>
 447
 448                 <para>The inputs and outputs of the function are the same as
 449                 in the fortran case.</para>
 450               </listitem>
 451             </varlistentry>
 452           </variablelist>
 453         </listitem>
 454       </varlistentry>
 455     </variablelist>
 456
 457     <para>If <literal>ind=2</literal> (resp. <literal>3, 4</literal>),
 458     <literal>costf</literal> must provide <literal>f</literal> (resp.
 459     <literal>g, f</literal> and <literal>g</literal>).</para>
 460
 461     <para>If <literal>ind=1</literal> nothing is computed (used for display
 462     purposes only).</para>
 463
 464     <para>On output, <literal>ind&lt;0</literal> means that
 465     <literal>f</literal> cannot be evaluated at <literal>x</literal> and
 466     <literal>ind=0</literal> interrupts the optimization.</para>
 467   </refsection>
 468
 469   <refsection>
 470     <title>Example #1 : Scilab function</title>
 471
 472     <para>The following is an example with a Scilab function. Notice, for
 473     simplifications reasons, the Scilab function "cost" of the following
 474     example computes the objective function f and its derivative no matter of
 475     the value of ind. This allows to keep the example simple. In practical
 476     situations though, the computation of "f" and "g" may raise performances
 477     issues so that a direct optimization may be to use the value of "ind" to
 478     compute "f" and "g" only when needed.</para>
 479
 480     <programlisting role="example"><![CDATA[
 481 // External function written in Scilab
 482 xref=[1;2;3];
 483 x0=[1;-1;1];
 484 function [f,g,ind] = cost(x,ind)
 485   f=0.5*norm(x-xref)^2;
 486   g=x-xref;
 487 endfunction
 488
 489 // Simplest call
 490 [f,xopt]=optim(cost,x0)
 491
 492 // By conjugate gradient - you can use 'qn', 'gc' or 'nd'
 493 [f,xopt,gopt]=optim(cost,x0,'gc')
 494
 495 //Seen as non differentiable
 496 [f,xopt,gopt]=optim(cost,x0,'nd')
 497
 498 // Upper and lower bounds on x
 499 [f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0)
 500
 501 // Upper and lower bounds on x and setting up the algorithm to 'gc'
 502 [f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0,'gc')
 503
 504 // Bound on the number of call to the objective function
 505 [f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0,'gc','ar',3)
 506
 507 // Set max number of call to the objective function (3)
 508 // Set max number of iterations (100)
 509 // Set stopping threshold on the value of f (1e-6),
 510 // on the value of the norm of the gradient of the objective function (1e-6)
 511 // on the improvement on the parameters x_opt (1e-6;1e-6;1e-6)
 512 [f,xopt,gopt]=optim(cost,'b',[-1;0;2],[0.5;1;4],x0,'gc','ar',3,100,1e-6,1e-6,[1e-3;1e-3;1e-3])
 513
 514 // Print information messages while optimizing
 515 // Be careful, some messages are printed in a terminal. You must
 516 // Scilab from the command line to see these messages.
 517 [f,xopt]=optim(cost,x0,imp=3)
 518
 519 // Use the 'derivative' function to compute the partial
 520 // derivatives of the previous problem
 521 deff('y=my_f(x)','y=0.5*norm(x-xref)^2');
 522 deff('y=my_df(x)','y=derivative(my_f,x)');
 523 deff('[f,g,ind]=cost(x,ind)','f=my_f(x); ...
 524                               g=my_df(x)');
 525
 526 // Simplest call
 527 xref=[1;2;3];x0=[1;-1;1]
 528 [f,xopt]=optim(cost,x0)
 529  ]]></programlisting>
 530   </refsection>
 531
 532   <refsection>
 533     <title>Example #2 : C function</title>
 534
 535     <para>The following is an example with a C function, where a C source code
 536     is written into a file, dynamically compiled and loaded into Scilab, and
 537     then used by the "optim" solver. The interface of the "rosenc" function is
 538     fixed, even if the arguments are not really used in the cost function.
 539     This is because the underlying optimization solvers must assume that the
 540     objective function has a known, constant interface. In the following
 541     example, the arrays ti and tr are not used, only the array "td" is used,
 542     as a parameter of the Rosenbrock function. Notice that the content of the
 543     arrays ti and td are the same that the content of the Scilab variable, as
 544     expected.</para>
 545
 546     <programlisting role="example"><![CDATA[
 547 // External function written in C (C compiler required)
 548 // write down the C code (Rosenbrock problem)
 549 C=['#include <math.h>'
 550    'double sq(double x)'
 551    '{ return x*x;}'
 552    'void rosenc(int *ind, int *n, double *x, double *f, double *g, '
 553    '                                int *ti, float *tr, double *td)'
 554    '{'
 555    '  double p;'
 556    '  int i;'
 557    '  p=td[0];'
 558    '  if (*ind==2||*ind==4) {'
 559    '    *f=1.0;'
 560    '    for (i=1;i<*n;i++)'
 561    '      *f+=p*sq(x[i]-sq(x[i-1]))+sq(1.0-x[i]);'
 562    '  }'
 563    '  if (*ind==3||*ind==4) {'
 564    '    g[0]=-4.0*p*(x[1]-sq(x[0]))*x[0];'
 565    '    for (i=1;i<*n-1;i++)'
 566    '      g[i]=2.0*p*(x[i]-sq(x[i-1]))-4.0*p*(x[i+1]-sq(x[i]))*x[i]-2.0*(1.0-x[i]);'
 567    '    g[*n-1]=2.0*p*(x[*n-1]-sq(x[*n-2]))-2.0*(1.0-x[*n-1]);'
 568    '  }'
 569    '}'];
 570 mputl(C,TMPDIR+'/rosenc.c')
 571
 572 // compile the C code
 573 l=ilib_for_link('rosenc','rosenc.o',[],'c',TMPDIR+'/Makefile');
 574
 575 // incremental linking
 576 link(l,'rosenc','c')
 577
 578 //solve the problem
 579 x0=[40;10;50];
 580 p=100;
 581 [f,xo,go]=optim('rosenc',x0,'td',p)
 582  ]]></programlisting>
 583   </refsection>
 584
 585   <refsection>
 586     <title>Example #3 : Fortran function</title>
 587
 588     <para>The following is an example with a Fortran function.</para>
 589
 590     <programlisting role="example"><![CDATA[
 591 // External function written in Fortran (Fortran compiler required)
 592 // write down the Fortran  code (Rosenbrock problem)
 593 F=[ '      subroutine rosenf(ind, n, x, f, g, ti, tr, td)'
 594     '      integer ind,n,ti(*)'
 595     '      double precision x(n),f,g(n),td(*)'
 596     '      real tr(*)'
 597     'c'
 598     '      double precision y,p'
 599     '      p=td(1)'
 600     '      if (ind.eq.2.or.ind.eq.4) then'
 601     '        f=1.0d0'
 602     '        do i=2,n'
 603     '          f=f+p*(x(i)-x(i-1)**2)**2+(1.0d0-x(i))**2'
 604     '        enddo'
 605     '      endif'
 606     '      if (ind.eq.3.or.ind.eq.4) then'
 607     '        g(1)=-4.0d0*p*(x(2)-x(1)**2)*x(1)'
 608     '        if(n.gt.2) then'
 609     '          do i=2,n-1'
 610     '            g(i)=2.0d0*p*(x(i)-x(i-1)**2)-4.0d0*p*(x(i+1)-x(i)**2)*x(i)'
 611     '     &           -2.0d0*(1.0d0-x(i))'
 612     '          enddo'
 613     '        endif'
 614     '        g(n)=2.0d0*p*(x(n)-x(n-1)**2)-2.0d0*(1.0d0-x(n))'
 615     '      endif'
 616     '      return'
 617     '      end'];
 618
 619 mputl(F,TMPDIR+'/rosenf.f')
 620
 621 // compile the Fortran code
 622 l=ilib_for_link('rosenf','rosenf.o',[],'f',TMPDIR+'/Makefile');
 623
 624 // incremental linking
 625 link(l,'rosenf','f')
 626
 627 //solve the problem
 628 x0=[40;10;50];
 629 p=100;
 630 [f,xo,go]=optim('rosenf',x0,'td',p)
 631  ]]></programlisting>
 632   </refsection>
 633
 634   <refsection>
 635     <title>Example #4 : Fortran function with initialization</title>
 636
 637     <para>The following is an example with a Fortran function in which the
 638     "in" option is used to allocate memory inside the Scilab environment. In
 639     this mode, there is a dialog between Scilab and the objective function.
 640     The goal of this dialog is to initialize the parameters of the objective
 641     function. Each part of this dialog is based on a specific value of the
 642     "ind" parameter.</para>
 643
 644     <para>At the beginning, Scilab calls the objective function, with the ind
 645     parameter equals to 10. This tells the objective function to initialize
 646     the sizes of the arrays it needs by setting the nizs, nrzs and ndzs
 647     integer parameters of the "nird" common. Then the objective function
 648     returns. At this point, Scilab creates internal variables and allocate
 649     memory for the variable izs, rzs and dzs. Scilab calls the objective
 650     function back again, this time with ind equals to 11. This tells the
 651     objective function to initialize the arrays izs, rzs and dzs. When the
 652     objective function has done so, it returns. Then Scilab enters in the real
 653     optimization mode and calls the optimization solver the user requested.
 654     Whenever the objective function is called, the izs, rzs and dzs arrays
 655     have the values that have been previously initialized.</para>
 656
 657     <programlisting role="example"><![CDATA[
 658 //
 659 // Define a fortran source code and compile it (fortran compiler required)
 660 //
 661 fortransource=['      subroutine rosenf(ind,n,x,f,g,izs,rzs,dzs)'
 662                'C     -------------------------------------------'
 663                'c     Example of cost function given by a subroutine'
 664                'c     if n<=2 returns ind=0'
 665                'c     f.bonnans, oct 86'
 666                '      implicit double precision (a-h,o-z)'
 667                '      real rzs(1)'
 668                '      double precision dzs(*)'
 669                '      dimension x(n),g(n),izs(*)'
 670                '      common/nird/nizs,nrzs,ndzs'
 671                '      if (n.lt.3) then'
 672                '        ind=0'
 673                '        return'
 674                '      endif'
 675                '      if(ind.eq.10) then'
 676                '         nizs=2'
 677                '         nrzs=1'
 678                '         ndzs=2'
 679                '         return'
 680                '      endif'
 681                '      if(ind.eq.11) then'
 682                '         izs(1)=5'
 683                '         izs(2)=10'
 684                '         dzs(2)=100.0d+0'
 685                '         return'
 686                '      endif'
 687                '      if(ind.eq.2)go to 5'
 688                '      if(ind.eq.3)go to 20'
 689                '      if(ind.eq.4)go to 5'
 690                '      ind=-1'
 691                '      return'
 692                '5     f=1.0d+0'
 693                '      do 10 i=2,n'
 694                '        im1=i-1'
 695                '10      f=f + dzs(2)*(x(i)-x(im1)**2)**2 + (1.0d+0-x(i))**2'
 696                '      if(ind.eq.2)return'
 697                '20    g(1)=-4.0d+0*dzs(2)*(x(2)-x(1)**2)*x(1)'
 698                '      nm1=n-1'
 699                '      do 30 i=2,nm1'
 700                '        im1=i-1'
 701                '        ip1=i+1'
 702                '        g(i)=2.0d+0*dzs(2)*(x(i)-x(im1)**2)'
 703                '30      g(i)=g(i) -4.0d+0*dzs(2)*(x(ip1)-x(i)**2)*x(i) - '
 704                '     &        2.0d+0*(1.0d+0-x(i))'
 705                '      g(n)=2.0d+0*dzs(2)*(x(n)-x(nm1)**2) - 2.0d+0*(1.0d+0-x(n))'
 706                '      return'
 707                '      end'];
 708 mputl(fortransource,TMPDIR+'/rosenf.f')
 709
 710 // compile the C code
 711 libpath=ilib_for_link('rosenf','rosenf.o',[],'f',TMPDIR+'/Makefile');
 712
 713 // incremental linking
 714 linkid=link(libpath,'rosenf','f');
 715
 716 x0=1.2*ones(1,5);
 717 //
 718 // Solve the problem
 719 //
 720 [f,x,g]=optim('rosenf',x0,'in');
 721  ]]></programlisting>
 722   </refsection>
 723
 724   <refsection>
 725     <title>Example #5 : Fortran function with initialization on Windows with
 726     Intel Fortran Compiler</title>
 727
 728     <para>Under the Windows operating system with Intel Fortran Compiler, one
 729     must carefully design the fortran source code so that the dynamic link
 730     works properly. On Scilab's side, the optimization component is
 731     dynamically linked and the symbol "nird" is exported out of the
 732     optimization dll. On the cost function's side, which is also dynamically
 733     linked, the "nird" common must be imported in the cost function
 734     dll.</para>
 735
 736     <para>The following example is a re-writing of the previous example, with
 737     special attention for the Windows operating system with Intel Fortran
 738     compiler as example. In that case, we introduce additionnal compiling
 739     instructions, which allows the compiler to import the "nird"
 740     symbol.</para>
 741
 742     <programlisting role="example"><![CDATA[
 743 fortransource=['subroutine rosenf(ind,n,x,f,g,izs,rzs,dzs)'
 744                'cDEC$ IF DEFINED (FORDLL)'
 745                'cDEC$ ATTRIBUTES DLLIMPORT:: /nird/'
 746                'cDEC$ ENDIF'
 747                'C     -------------------------------------------'
 748                'c     Example of cost function given by a subroutine'
 749                'c     if n<=2 returns ind=0'
 750                'c     f.bonnans, oct 86'
 751                '      implicit double precision (a-h,o-z)'
 752 [etc...]
 753  ]]></programlisting>
 754   </refsection>
 755
 756   <refsection>
 757     <title>Example #6 : Logging features</title>
 758
 759     <para>The imp flag may take negative integer values, say k.
 760     In that case, the cost function is called once every -k iterations.
 761     This allows to draw the function value or write a log file.
 762     </para>
 763
 764     <para>In the following example, we solve the Rosenbrock test case.
 765     For each iteration of the algorithm, we print the value of x, f and g.
 766     </para>
 767
 768 <programlisting role="example">
 769 xref=[1;2;3];
 770 x0=[1;-1;1];
 771 function [f,g,ind] = cost(x,ind)
 772   f=0.5*norm(x-xref)^2;
 773   g=x-xref;
 774 endfunction
 775 function [f,g,ind] = cost(x,ind)
 776   if ind == 2 | ind == 4 then
 777     f=0.5*norm(x-xref)^2;
 778   end
 779   if ind == 3 | ind == 4 then
 780     g=x-xref;
 781   end
 782   if ind == 1 then
 783     mprintf("x = %s\n", strcat(string(x)," "))
 784     mprintf("f = %e\n", f)
 785     g=x-xref;
 786     mprintf("g = %s\n", strcat(string(g)," "))
 787   end
 788 endfunction
 789 [f,xopt]=optim(cost,x0,imp=-1)
 790   </programlisting>
 791
 792     <para>In the following example, we solve the Rosenbrock test case.
 793     For each iteration of the algorithm, we plot the current value of x
 794     into a 2D graph containing the contours of Rosenbrock's function.
 795     This allows to see the progress of the algorithm while the
 796     algorithm is performing. We could as well write the
 797     value of x, f and g into a log file if needed.
 798     </para>
 799
 800 <programlisting role="example">
 801 // 1. Define rosenbrock
 802 function [ f , g , ind ] = rosenbrock ( x , ind )
 803   if ((ind == 1) | (ind == 2) | (ind == 4)) then
 804     f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
 805   end
 806   if ((ind == 1) | (ind == 3) | (ind == 4)) then
 807     g(1) = - 400. * ( x(2) - x(1)**2 ) * x(1) -2. * ( 1. - x(1) )
 808     g(2) = 200. * ( x(2) - x(1)**2 )
 809   end
 810   if (ind == 1) then
 811     plot ( x(1) , x(2) , "g." )
 812   end
 813 endfunction
 814 x0 = [-1.2 1.0];
 815 xopt = [1.0 1.0];
 816 // 2. Draw the contour of Rosenbrock's function
 817 xmin = -2.0
 818 xmax = 2.0
 819 stepx = 0.1
 820 ymin = -1.0
 821 ymax = 2.0
 822 stepy = 0.1
 823 nx = 100
 824 ny = 100
 825 stepx = (xmax - xmin)/nx;
 826 xdata = xmin:stepx:xmax;
 827 stepy = (ymax - ymin)/ny;
 828 ydata = ymin:stepy:ymax;
 829 for ix = 1:length(xdata)
 830     for iy = 1:length(ydata)
 831       x = [xdata(ix) ydata(iy)];
 832       f = rosenbrock ( x , 2 );
 833       zdata ( ix , iy ) = f;
 834     end
 835 end
 836 contour ( xdata , ydata , zdata , [1 10 100 500 1000])
 837 plot(x0(1) , x0(2) , "b.")
 838 plot(xopt(1) , xopt(2) , "r*")
 839 // 3. Plot the optimization process, during optimization
 840 [ fopt , xopt ] = optim ( rosenbrock , x0 , imp = -1)
 841   </programlisting>
 842   </refsection>
 843
 844
 845   <refsection>
 846     <title>Example #7 : Optimizing without derivatives</title>
 847
 848     <para>It is possible to optimize a problem without an explicit
 849     knowledge of the derivative of the cost function.
 850     For this purpose, we can use the numdiff or derivative function
 851     to compute a numerical derivative of the cost function.
 852     </para>
 853
 854     <para>In the following example, we use the numdiff function to
 855     solve Rosenbrock's problem.
 856     </para>
 857
 858
 859 <programlisting role="example">
 860 function f = rosenbrock ( x )
 861   f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
 862 endfunction
 863 function [ f , g , ind ] = rosenbrockCost ( x , ind )
 864   if ((ind == 1) | (ind == 2) | (ind == 4)) then
 865     f = rosenbrock ( x );
 866   end
 867   if ((ind == 1) | (ind == 3) | (ind == 4)) then
 868     g= numdiff ( rosenbrock , x );
 869   end
 870 endfunction
 871 x0 = [-1.2 1.0];
 872 [ fopt , xopt ] = optim ( rosenbrockCost , x0 )
 873   </programlisting>
 874
 875     <para>In the following example, we use the derivative function to
 876     solve Rosenbrock's problem. Given that the step computation strategy
 877     is not the same in numdiff and derivative, this might lead to improved
 878     results.
 879     </para>
 880
 881
 882 <programlisting role="example">
 883 function f = rosenbrock ( x )
 884   f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
 885 endfunction
 886 function [ f , g , ind ] = rosenbrockCost2 ( x , ind )
 887   if ((ind == 1) | (ind == 2) | (ind == 4)) then
 888     f = rosenbrock ( x );
 889   end
 890   if ((ind == 1) | (ind == 3) | (ind == 4)) then
 891     g= derivative ( rosenbrock , x.' , order = 4 );
 892   end
 893 endfunction
 894 x0 = [-1.2 1.0];
 895 [ fopt , xopt ] = optim ( rosenbrockCost2 , x0 )
 896   </programlisting>
 897
 898   </refsection>
 899
 900
 901   <refsection>
 902     <title>See Also</title>
 903
 904     <simplelist type="inline">
 905       <member><link linkend="external">external</link></member>
 906
 907       <member><link linkend="qpsolve">qpsolve</link></member>
 908
 909       <member><link linkend="datafit">datafit</link></member>
 910
 911       <member><link linkend="leastsq">leastsq</link></member>
 912
 913       <member><link linkend="numdiff">numdiff</link></member>
 914
 915       <member><link linkend="derivative">derivative</link></member>
 916
 917       <member><link linkend="NDcost">NDcost</link></member>
 918     </simplelist>
 919   </refsection>
 920
 921   <refsection>
 922     <title>References</title>
 923
 924     <para>The following is a map from the various options to the underlying
 925     solvers, with some comments about the algorithm, when available.</para>
 926
 927     <variablelist>
 928       <varlistentry>
 929         <term>"qn" without constraints</term>
 930
 931         <listitem>
 932           <para>n1qn1 : a quasi-Newton method with a Wolfe-type line
 933           search</para>
 934         </listitem>
 935       </varlistentry>
 936
 937       <varlistentry>
 938         <term>"qn" with bounds constraints</term>
 939
 940         <listitem>
 941           <para>qnbd : a quasi-Newton method with projection</para>
 942
 943           <para>RR-0242 - A variant of a projected variable metric method for
 944           bound constrained optimization problems, Bonnans Frederic, Rapport
 945           de recherche de l'INRIA - Rocquencourt, Octobre 1983</para>
 946         </listitem>
 947       </varlistentry>
 948
 949       <varlistentry>
 950         <term>"gc" without constraints</term>
 951
 952         <listitem>
 953           <para>n1qn3 : a conjugate gradient method with BFGS.</para>
 954         </listitem>
 955       </varlistentry>
 956
 957       <varlistentry>
 958         <term>"gc" with bounds constraints</term>
 959
 960         <listitem>
 961           <para>gcbd : a BFGS-type method with limited memory and
 962           projection</para>
 963         </listitem>
 964       </varlistentry>
 965
 966       <varlistentry>
 967         <term>"nd" without constraints</term>
 968
 969         <listitem>
 970           <para>n1fc1 : a bundle method</para>
 971         </listitem>
 972       </varlistentry>
 973
 974       <varlistentry>
 975         <term>"nd" with bounds constraints</term>
 976
 977         <listitem>
 978           <para>not available</para>
 979         </listitem>
 980       </varlistentry>
 981     </variablelist>
 982   </refsection>
 983   <refsection>
 984     <title>Author</title>
 985     <para>The Modulopt library : J.Frederic Bonnans, Jean-Charles Gilbert, Claude Lemarechal</para>
 986     <para>The interfaces to the Modulopt library : J.Frederic Bonnans</para>
 987     <para>This help : Michael Baudin</para>
 988   </refsection>
 989 </refentry>